计算机仿真教案 39页

研究生课程计算机仿真讲稿计算机仿真教材：《控制系统计算机仿真与CAD—MATLAB语言应用》，天津大学出版社参考书：1．《控制系统计算机辅助设计》，清华大学出版社；2．《掌握和精通MATLAB》，北京航空航天大学出版社课程学习的目的：1.了解和掌握自动控制系统计算机仿真的基本原理、方法；2.采用MATLAB工具对控制系统进行仿真；3.采用MATLAB工具对控制系统进行辅助设计。前期课程基础：自动控制原理、自动控制系统、C语言基础授课学时与分配：总学时30学时1．绪论、控制系统数学模型的建立、转换及其连接（3学时）2．控制系统的数字仿真1（3学时）3．控制系统的数字仿真2（3学时）4．最优化的计算机辅助设计概述、线性规划（3学时）5．单变量寻优（6学时）7．多变量寻优（6学时）8．PID参数寻优（3学时）9．总结（2学时）授课时间、地点：时间为；地点为科技楼3楼多媒体教室。39\n研究生课程计算机仿真讲稿第1章绪论一、仿真的目的与方法仿真的目的：采用物理模拟和计算（数值计算、模拟计算）等仿真方法，分析系统的动、静态性能。仿真是对于复杂控制系统进行分析的重要手段。仿真的方法采用物理仿真、数学仿真的前题是数学相似——即仿真模型与对象的传递函数相同。二、数字仿真工具MATLAB的发展历史。创始人：CleveMoler教授，MathWorks公司的首席科学家。70年代中CleveMoler教授在美国国家科学基金会的资助下开发了两个FORTRAN子程序库。LINPACK—求解线性方程；EISPACK—解特征值。70年代末CleveMoler教授为了便于学生使用该程序库，设计了调用接口程序。取名为MATLAB（MATrixLABoratory）。84年MathWorks公司将MATLAB正式推向市场。MATLAB版本升级历程日期版本平台1987MATLAB3.0DOS1991MATLAB3.5DOS1993MATLAB3.0kWINDOWS3.01993MATLAB4.0WINDOWS3.x1993MATLAB4.1WINDOWS3.x1994MATLAB4.2WINDOWS3.x1995MATLAB4.2CWINDOWS3.x1997MATLAB5.0WINDOWS95三、MATLAB简介1．变量与表达式；?A=3*16+10-27/3A=49?3*16+10-27/3ans=39\n研究生课程计算机仿真讲稿491．基本运算符；运算包括算术运算和逻辑运算，参见书上表1-23．M文件例P1_004.M39研究生课程计算机仿真讲稿%2s+4%-------------------------%s^4+5s^3+8s^2+6sclcnum=2*[1,2];den=conv(conv([1,0],[1,3]),[1,2,2]);g1=tf(num,den);g=ss(g1);[a,b,c,d]=ssdata(g);ab=a-b*c;bb=b;cb=c;db=d;step(ab,bb,cb,db);39研究生课程计算机仿真讲稿4．矩阵、数组与数值运算功能；1）矩阵的输入?A=[123;4,5,6;789]A=123456789另外还可以用MATLAB的函数构造一些特殊的矩阵。2）矩阵与数组的运算A+B、A-B、A*B、A.*B、A\B、A.\B、inv(A)、A^k、矩阵转置、矩阵翻转等。5．多项式运算1）求根r=root(p)已知；p=[1–1141–6130](P1_014)%X^4-11X^3+41X^2-61X+30=0的根clcp=[1-1141-6130]r=roots(p)2)poly函数p=poly(r)39\n研究生课程计算机仿真讲稿(p1_015)clcA=[12;34]p=poly(A)r=roots(p)4．MATLAB符号运算功能(P1_FH)f1='(41*x^2+x+1)^2*(2*x-1)';f2='((3*x+5)*(2*x-1))';f=symdiv(f1,f2);solve(f,'x');例如：求f1/f2的表达式f。并求f=0的解。5．MATLAB语句流1）循环语句例如：求的值，可以用以下两种方法实现39研究生课程计算机仿真讲稿（1）sumi=0；fori=1:100:1sumi=sumi+i;endsumi运行结果是sumi=5050（2）sumi=0;i=1while(i<1000)sumi=sumi+i;i=i+1endsumi运行结果是sumi=505039研究生课程计算机仿真讲稿2)条件转移语句例如：当键入Y时给x赋值1，当键入N时给x赋值0，键入其它字符退出程序。ikey=0;while(ikey==0)s1=input(‘请键入字符’,’s’);if(s1==’Y’)x=1elseif(si==’N’)x=0elsebreak39\n研究生课程计算机仿真讲稿endend4．MATLAB的典型时域响应函数1）阶跃响应函数：step(sys);step(sys,T);step(num,den,T);step(A,B,C,D,T)T是一个时间或时间向量。如T=20或T=0:0.01:20例：求的单位阶跃响应（P1_020）39研究生课程计算机仿真讲稿clcholdonnum=[16];den=[14361817];T=20;sys=tf(num,den);step(sys,2,T);gridinputa;T=0:2:20;step(sys,T);gridholdoff39研究生课程计算机仿真讲稿2）任意输入下的响应函数：lsim(sys,U,T);lsim(sys,U,T,X0);lsim(num,den,U,T);lsim(A,B,C,D,U,T)例：，，，。（P1_021）39研究生课程计算机仿真讲稿clcA=[011;001;-4-3-2];B=[005]';C=[100];D=[0];sys=ss(A,B,C,D);t=0:0.001:10;u=20+0.4*sin(t);lsim(sys,u,t);grid39研究生课程计算机仿真讲稿5．MATLAB的图形绘制1）plot绘图函数。和传统绘图方法的区别（p1_022）39研究生课程计算机仿真讲稿clcholdony0=0;t0=0;fort=1:0.01:4*piy1=sin(t);y=[y0,y1];x=[t0,t];plot(x,y);t0=t;y0=y1;endgridoninputa;t=0:0.01:4*pi;y=sin(t+pi/6);39\n研究生课程计算机仿真讲稿plot(t,y,'m');holdoff39研究生课程计算机仿真讲稿2）同窗口下多幅曲线的绘制subplot(m,n,i)(p1_023)39研究生课程计算机仿真讲稿clcholdony0=0;t0=0;t=0:0.01:4*pi;y=sin(t);subplot(2,1,1);plot(t,y);gridont=0:0.01:4*pi;y=sin(t+pi/6);subplot(2,1,2);plot(t,y,'m');gridonholdoff39研究生课程计算机仿真讲稿3)坐标调整与图形标注四Simulink（P1_026）ScopeStepSumTransferFcn20.5s+s2439\n研究生课程计算机仿真讲稿39\n研究生课程计算机仿真讲稿第2章控制系统数学模型、模型之间的转换与连接一、控制系统的数学描述1．微分方程模型2．连续系统的传递函数模型传递函数模型的分子、分母定义：，建立模型3．离散系统脉冲传递函数模型传递函数模型的分子、分母定义：，建立模型4．连续系统零极点形式模型，5．连续系统的状态空间表达式A=[a11，a12，…a1n；a21，a22，…a2n；…an1，an2，…ann]；B=[b1，b2，…bn]；C=[c1，c2，…cn]；D=d；模型的基本格式：sys=ss（A，B，C，D）6．离散系统的状态空间表达式A=[a11，a12，…a1n；a21，a22，…a2n；…an1，an2，…ann]；B=[b1，b2，…bn]；C=[c1，c2，…cn]；D=d；模型的基本格式：sys=ss（A，B，C，D，T）二、系统模型的相互转换1．状态空间方程与传递函数之间的相互转换1）状态空间方程TO传递函数39\n研究生课程计算机仿真讲稿，转换方法1：[num，den]=ss2tf（A，B，C，D，iu）转换方法2：G1=ss（A，B，C，D）；G2=tf（G1）；2）传递函数TO状态空间方程转换方法：[A，B，C，D]=tf2ss（num，den）(p2_038)clc39研究生课程计算机仿真讲稿a=[200;041;004];b=[101]';c=[110];d=0;[num,den]=ss2tf(a,b,c,d);g=tf(num,den)inputl;g1=ss(a,b,c,d);g=tf(g1)end[a,b,c,d]=tf2ss(num,den)end39研究生课程计算机仿真讲稿2．状态空间方程与零极点模型之间的相互转换1）状态空间方程TO零极点模型转换方法1：[z，p，k]=ss2zp（A，B，C，D，IU）转换方法2：sys=ss（A，B，C，D）；Gzp=zpk（sys）2）零极点模型TO状态空间方程[A，B，C，D]=zp2ss（z，p，k）(p2_039)39研究生课程计算机仿真讲稿clca=[200;041;004];b=[101]';c=[110];d=0;[z,p,k]=ss2zp(a,b,c,d)inputl;g1=ss(a,b,c,d);g=zpk(g1)%[a,b,c,d]=zp2ss(z,p,k)end39研究生课程计算机仿真讲稿3．传递函数与零极点模型之间的相互转换1）传递函数TO零极点模型[z，p，k]=tf2zp（num，den）2）零极点模型TO传递函数[num，den]=zp2tf（z，p，k）(p2_040)39研究生课程计算机仿真讲稿clcnum=[2,1];den=[132];g1=tf(num,den)[z,p,k]=tf2zp(num,den)inputlg=zpk(g1)[num,den]=zp2tf(z,p,k);tf(num,den)39研究生课程计算机仿真讲稿4．连续系统和离散系统之间的转换1）连续系统TO离散系统先转换成为状态方程，再离散。[A，B，C，D]=tf2ss（num，den）G=ss（A，B，C，D，T）；sys=tf（G）39\n研究生课程计算机仿真讲稿如果要求离散时采用具体的方法：sysd=c2d（sysc，T，method）。method的定义为：‘zoh’零阶保持器；‘foh’：三角型保持器等。如果修改采样周期，格式为：sysd2=d2d（sysd1，T，method）2）离散系统TO连续系统sysc=d2c（sysd，method）例：要求：1将其离散化处理；2采用零阶保持器离散化，采样周期T=0.2；3改变T=0.4再次离散；4连续化.(p2-042)39研究生课程计算机仿真讲稿clcnum=[12];den=[.58.5115];T=.2;[A,B,C,D]=tf2ss(num,den);G=ss(A,B,C,D,T);sys=tf(G)inputl;G=ss(A,B,C,D);sysd=c2d(G,T,'zoh');sys=tf(sysd)inputl;t=.4;sys1=d2d(sys,t)inputl;sysc=d2c(sys,'zoh')inputl;sysc=d2c(sys1,'zoh')39研究生课程计算机仿真讲稿三、系统状态方程的变换与实现1．状态方程的相似变换已知，设，则。MATLAB实现方法GT=ss2ss（G，T）2．规范型状态方程的实现；3.系统的均衡实现；4.降阶四、控制系统的模型建立与典型连接1．基本模型1）二阶模型的建立[num，den]=ord2（Wn，z）或[A，B，C，D]=ord2（Wn，z）2）随机n阶系统的模型建立[num，den]=rmodel（n[，p，[m]]）[num，den]=drmodel（n[，p，[m]]）[A，B，C，D]=rmodel（n[，p，[m]]）[A，B，C，D]=drmodel（n[，p，[m]]）rmodel产生n阶p出m入的稳定的传递函数（状态方程）模型；drmodel产生相应的离散模型。3）模型结构的重构39\n研究生课程计算机仿真讲稿子系统的选取：[Ae，Be，Ce，De]=ssselect（A，B，C，D，inputs，outputs[，states]）子系统的删除：[Ae，Be，Ce，De]=ssdelect（A，B，C，D，inputs，outputs[，states]）状态增广：将变换为2．系统组合与连接模型串联：sys=series（sys1，sys2）模型并联：sys=parallel（sys1，sys2）反馈连接：sys=feedback（sys1，sys2）39\n研究生课程计算机仿真讲稿第3章控制系统数字仿真一、基本数值积分方法解决的问题求1．欧拉（Euler）法1）基本原理，所以。令，，则。描述成一般形式：2）物理意义与误差分析，其中，结论：欧拉法的实质为，用巨型面积替代曲边型的面积。2．梯形法原理：用体形面积替代曲边型的面积。。可得一般形式。此式中，求时，用到。解决方法是：，。3．龙格—库塔方法原理：将在附近展开，得到。将，代入，可得。设上式可以表示为：，式中，。将在附近展开，可得39\n研究生课程计算机仿真讲稿。比较可得；。因此写成一般形式。此式被称为二阶龙格—库塔方法。四阶龙格—库塔方法为4．四阶龙格—库塔方法的向量表示法状态方程一般表示方法是。则；；；；。例。求阶跃响应（P3_064）39研究生课程计算机仿真讲稿clca=[-5-2-1-.5;4000;0200;0010];b=[1000]';c=[00.25.5];d=0;u=10;x=[0000]';y=0;t=0;ts=15;h=0.01;n=ts/h;fori=1:nk1=a*x+b*u;k2=a*(x+h/2*k1)+b*u;k3=a*(x+h/2*k2)+b*u;k4=a*(x+h*k3)+b*u;x=x+h*(k1+2*k2+2*k3+k4)/6;y=[y,c*x];t=[t,t(i)+h];endplot(t,y);gridon;39研究生课程计算机仿真讲稿39\n研究生课程计算机仿真讲稿5．亚当斯法二、数值积分方法的计算稳定性1．基本概念一个稳定的系统当采用某种数值积分方法进行计算时，若选择的参数不当，则会产生意想不到的结果。例如：，。用欧拉法进行仿真计算。解：根据传递函数，可得微分方程。。讨论的取值对仿真结果的影响。（P3_065）39研究生课程计算机仿真讲稿clcholdony=dsolve('Dy=1-y','y(0)=0','t');t=0:.1:20;yy=subs(y,t,'t');yy1=numeric(yy);plot(t,yy1)gridonh=2;m=20/h;y0=0;t=[0];y=[0];fori=2:m;y(i)=y0+h*(1-y0);t(i)=t(i-1)+h;y0=y(i);endplot(t,y,'m')gridon39研究生课程计算机仿真讲稿2．各种方法稳定性分析三、连续系统仿真1．数值积分的MATLAB实现1）欧拉法的MATLAB实现例：已知用欧拉法求时的阶跃响应（25秒），。（P3_067）39研究生课程计算机仿真讲稿clch=0.01;x10=0;x20=0;u=25;t=[0];y=[0];m=25/h;fori=1:mx1=x10+h*(-.5572*x10-.7814*x20+u);x2=x20+h*(.7814*x10);y(i+1)=1.9691*x1+6.4493*x2;x10=x1;x20=x2;t(i+1)=t(i)+h;endplot(t,y)gridon39研究生课程计算机仿真讲稿不对39\n研究生课程计算机仿真讲稿2）梯形法的MAYTLAB实现解决上述问题。（P3_068）39研究生课程计算机仿真讲稿u=25;t=[0];y=[0];m=25/h;fori=1:mx1yb=x10+h*(-.5572*x10-.7814*x20+u);x2yb=x20+h*(.7814*x10);x1=x10+h/2*((-.5572*x10-.7814*x20+u)+(-.5572*x1yb-.7814*x2yb+u));x2=x20+h/2*(.7814*x10+.7814*x1yb);y(i+1)=1.9691*x1+6.4493*x2;x10=x1;x20=x2;t(i+1)=t(i)+h;endplot(t,y)gridon39研究生课程计算机仿真讲稿3）龙格—库塔法的MAYTLAB实现解决上述问题。（P3_069）39研究生课程计算机仿真讲稿clch=0.65;x=[0;0];u=25;t=[0];y=[0];m=25/h;a=[-.5572-.7814;.78140];b=[1;0];c=[1.96916.4493];d=0;fori=1:mk1=a*x+b*u;k2=a*(x+h*k1/2)+b*u;k3=a*(x+h*k2/2)+b*u;k4=a*(x+h*k3)+b*u;x=x+h/6*(k1+2*k2+2*k3+k4);y(i+1)=c*x+d*u;t(i+1)=t(i)+h;endplot(t,y)gridon39研究生课程计算机仿真讲稿2．面向微分方程的系统仿真与MATLAB实现1）基于ode函数的面向微分方程的系统仿真ode45函数[t,y]=ode45(‘F’,tspan,y0)；tspan：时间启始和终止时刻。y0：初试条件。例：，，，。求解。首先构造函数p3_071f，存放在p3_071f.m文件中。在编写p3_071.m程序计算。39研究生课程计算机仿真讲稿p3_071f.m内容functionxd=p3_071f(t,x)39\n研究生课程计算机仿真讲稿u=100;xd=[-14*x(1)+9*x(2)+10*x(3)+2.8*u;12*x(1)-9*x(2)+10*x(3)+4*u;24*x(1)-24*x(2)-18*x(3)+12*u];p3_071.m内容clcx0=[0;0;0];[t,x]=ode45('p3_071f',[0,3],x0);plot(t,x)gridon39研究生课程计算机仿真讲稿1）基于M函数的面向微分方程的系统仿真例：同上。首先构造函数p3_072f，存放在p3_072f.m文件中。在编写p3_072.m程序计算。39研究生课程计算机仿真讲稿p3_072f.m内容：functionxd=p3_072f(x,u)a=[-14910;12-910;24-24-18];b=[2.8;4;12];xd=a*x+b*u;p3_072.m内容：clch=0.01;t=0;m=3/h;x=[0;0;0];x1=[0,0,0];t=[0];u=100;fori=1:mk1=p3_072f(x,u);k2=p3_072f(x+h/2*k1,u);k3=p3_072f(x+h/2*k2,u);k4=p3_072f(x+h*k3,u);x=x+h/6*(k1+2*k2+2*k3+k4);x1=[x1;x'];t(i+1)=t(i)+h;endplot(t,x1);gridon39研究生课程计算机仿真讲稿3．面向传递函数的系统仿真与MATLAB实现系统的开环传递函数为，它对应的状态空间表达式为。当，代入得到，式中。例：单位反馈系统的开环传递函数为，求系统的阶跃响应，阶跃幅值为25。(p3_074)39研究生课程计算机仿真讲稿clcnum=2*[1,2];den=conv(conv([1,0],[1,3]),[1,2,2]);G=ss(tf(num,den));[aa,bb,cc,dd]=ssdata(G);h=0.01;x=[0;0;0;0];u=25;t=[0];y=[0];m=15/h;a=aa-bb*cc;b=bb;c=cc;d=dd;fori=1:mk1=a*x+b*u;k2=a*(x+h*k1/2)+b*u;39\n研究生课程计算机仿真讲稿k3=a*(x+h*k2/2)+b*u;k4=a*(x+h*k3)+b*u;x=x+h/6*(k1+2*k2+2*k3+k4);y(i+1)=c*x+d*u;t(i+1)=t(i)+h;endplot(t,y)gridon39研究生课程计算机仿真讲稿4．面向结构图的系统仿真1）基于典型环节的系统仿真a．系统典型化处理+-对于积分、惯性、一阶超前（滞后）环节可以描述为。对于二阶震荡环节可以处理为所有环节传递函数的一般形式为。将其写成矩阵形式为。其中，，b．连接矩阵，式中和为连接矩阵c．闭环系统的状态方程，式中；。例3-5(p3_079)39\n研究生课程计算机仿真讲稿--39研究生课程计算机仿真讲稿clcaa=diag([16.511]);bb=diag([4132]);cc=diag([412.45]);dd=diag([0001]);w=[00-40;100-1;0100;0010];w0=[1;0;0;0];Q1=bb-dd*w;Q=inv(Q1);R=cc*w-aa;a=Q*R;b=Q*cc*w0;c=[0010];d=0;h=0.01;x=[0;0;0;0];u=25;t=[0];y=[0];m=15/h;fori=1:mk1=a*x+b*u;k2=a*(x+h*k1/2)+b*u;k3=a*(x+h*k2/2)+b*u;k4=a*(x+h*k3)+b*u;x=x+h/6*(k1+2*k2+2*k3+k4);y(i+1)=c*x+d*u;t(i+1)=t(i)+h;endplot(t,y)gridon39研究生课程计算机仿真讲稿2）基于Connect连接函数的系统仿真例3-6(p3_080)39\n研究生课程计算机仿真讲稿81234567910--39研究生课程计算机仿真讲稿clcn=10;nn1=1;d1=[.011];nn2=[.171];d2=[.0850];nn3=1;d3=[.011];nn4=[.151];d4=[.0510];nn5=70;d5=[.00671];nn6=.21;d6=[.151];nn7=130;d7=[10];nn8=-0.212;d8=1;nn9=-.1;d9=[.011];nn10=-.0044;d10=[.011];blkbuild;Q=[100;2110;320;439;540;658;760;870;960;1070];[A,B,C,D]=connect(a,b,c,d,Q,1,7);G=ss(A,B,C,D);step(G,2);gridon39研究生课程计算机仿真讲稿5．离散相似法的系统仿真1）概念保持器u(t)u*(t)uh(t)y(t)y*(t)T139\n研究生课程计算机仿真讲稿2）面向闭环系统的方法已知，其中状态方程的解为。将上式写成一般形式，可得。采用零阶保持器时，代入可得其中：采用一阶保持器时。代入上式可以得到例：已知用离散相似法求时的阶跃响应（采用零阶保持器，）。(P3_084)39研究生课程计算机仿真讲稿clca=[-.5572-.7814;.78140];b=[1;0];c=[1.96916.4493];d=0;T=0.01;F=expm(a*T);G=inv(a)*(F-eye(size(F)))*b;t=0;y=0;x=[0;0];u=25;fori=2:25/Tx=F*x+G*u;y(i)=c*x+d*u;t(i)=t(i-1)+T;end39\n研究生课程计算机仿真讲稿plot(t,y);gridon39研究生课程计算机仿真讲稿3）非线性环节的处理例3—12（书上）(p3_090)39研究生课程计算机仿真讲稿clcT=0.05;F=[1,0;.5*(1-exp(-2*T)),exp(-2*T)];G=[T;T/2-1+exp(-2*T)];c=[01];t=0;y=0;x=[0;0];r=25;k=10;fori=2:30/Tu=r-y(i-1);ifu>ku=k;endifu<-ku=-k;endx=F*x+G*u;y(i)=c*x+d*u;t(i)=t(i-1)+T;endplot(t,y);gridon39研究生课程计算机仿真讲稿四、采样系统仿真1．采样系统组成控制器保持器r(t)e(t)e*(t)u(t)u*(t)uh(t)y(t)2．控制算法u(k)=u(k-1)+Kp*[e(k)-e(k-1)]+Ki*e(k)+Kd*[e(k)+e(k-2)-2*e(k-1)]3．控制对象的差分方程如，则，因为，可得。将其进行Z反变换得（P3_091）39研究生课程计算机仿真讲稿clck0=1;T0=1;T=.1;a=exp(-T/T0);r=10;y=0;t=0;u=0;e=0;e1=0;e2=0;kp=2;ki=.4;kd=1;fori=2:4/Te2=e1;e1=e;e=r-y(i-1);39\n研究生课程计算机仿真讲稿u1=u;u=u1+kp*(e-e1)+ki*e+kd*(e+e2-2*e1);y(i)=a*y(i-1)+k0*(1-a)*u1;t(i)=t(i-1)+T;endplot(t,y);gridon39研究生课程计算机仿真讲稿39\n研究生课程计算机仿真讲稿第六章控制系统最优化的计算机辅助设计最优控制是指：为了完成某一个目标，在一定的约束条件下，使得某项指标或某几项指标的综合达到极限（最大或最小）的控制方法和控制策略。一、概述最优控制系统的分类：最小误差控制系统；时间最优控制系统；能量最优控制系统；状态最优控制系统；最佳调度、最佳线路系统。1．目标函数积分型：；终值型：；综合型：2．约束条件等式约束和不等式约束。用数学描述为：1．参数（静态）最优化和函数（动态）最优化1）参数（静态）最优化控制器的结构已经确定，有个可以调整的控制参数，外加激励为。现在要求对控制器的参数进行优化，即找到一组最佳的具有固定数值的，使得达到最小值。2）和函数（动态）最优化如果上述问题变成寻找控制器输出，使得达到最小值。由于该问题要求确定的是最佳控制函数，它是虽时间变化的函数。所以这类问题称为函数的最优化问题，数学上称为范函极值问题。3）动态最优化转换为静态最优化在仿真研究中通常采用将离散，在一个时间段内，将一个时间的函数寻优问题，变成一个对非时间函数寻优问题，既把一个函数最优化的问题转化为参数最优化的问题。具体方法如下：设其中：；39\n研究生课程计算机仿真讲稿图6-3用矩形近似4．直接寻优与间接寻优1）间接寻优间接寻优的的方法是，根据给出的目标函数，利用数学解析的手段求取其极值。如已知的充分条件是剃度为零，即：（6-7）并且Hession矩阵为正定阵，即：为正定阵。按照这种方法，将对待寻优的参数分别求一次偏导并令其为零，这样接触一组，然后代入以检验其是否正定，如果正定则是一组能够使的最优解，这就是间接寻优。2）直接寻优直接寻优的思想是：从给定的初始状态起步，按照某种规则找到新的一点，然后计算是否达到要求的最小值，若是，则计算停止，否则，继续按照规定的规则寻找，直到满足要求为止。计算方法见6-8式。39\n研究生课程计算机仿真讲稿其中：为一个实数，被称为寻优步长；为一个矢量，被称为寻优方向。寻优的步骤见6-5所示的程序框图。寻优步长和寻优方向根据不同的直接寻优方法的不同取值的方法不同。二、线性规划方法例6-3设两个砖厂年产量分别为620万块，340万块。有三个工地需要用砖，年需求量分别为270万块，380万块，310万块。自产地运到销地的费用见表6-1。表6-1运输费用清单50元/万块60元/万块70元/万块60元/万块110元/万块95元/万块解设从砖厂调运到工地的砖数量为，（），总的运费为。以上问题用数学描述如下：+++++约束条件（）（6-8）例6-4某工厂在一个月内要求生产甲、乙两种产品。共需A、B、C、D四种设备加工这两种产品。两种产品在四种设备上加工的工时，以及设备在一个月内所能提供的工时数见表6-2。生产甲种产品可以获利200元，生产乙种产品可以获利300元。问如何安排生产才能够获得最大利润。表6-2工时供、需表ABCD甲2001004000乙2002000400最大工时120080016001200解甲、乙的计划产量分别为，总利润为。则上述问题的数学描述为39\n研究生课程计算机仿真讲稿（6-9）例6-3某农场要购买、、、四种型号的拖拉机，表6-3给出了各种型号拖拉机的价格和生产能力。农场现有农活为：春播5万亩；夏管3.5万亩；秋收8.5万亩。怎样作到最佳资的源配置。表6-3各种型号的拖拉机的价格和生产能力型号单价春播夏管秋收5000元450亩260亩620亩4500元440亩210亩640亩4400元480亩240亩630亩5200元460亩270亩650亩解设购置各种型号的拖拉机的台数分别为，总支出为。则上述问题可以描述为：（6-10）从以上三个例子看，线性规划解决的问题可以分为两大类：最大产值问题和最佳资源配置或最佳行动计划问题。前者为求目标函数的1．问题分类：线性规划解决的问题可以分为两大类：最大产值问题和最佳资源配置或最佳行动计划问题。前者为求目标函数的极大值，后者为求目标函数的极小值。约束条件可以分成等式约束条件和不等式约束条件。2．线性规划的标准型欲求极小值或极大值的函数称为目标函数；要求自变量应满足的条件称为约束条件。当目标函数和约束条件为线性表达式时，这一类问题被称为线性规划问题。线性规划的标准型式为：（6-11）39\n研究生课程计算机仿真讲稿其中；；；3．线性规划的图解法及一些基本概念6543210123456789（1）可行解集：满足约束条件所有点的集合称为可行解集，或称为凸集。（2）基本可行解集：约束条件所形成多边形，其顶点处的可行解被称为基本可行解。使得目标函数取到极值的点，必定是基本可行解中的一个。（3）最优基本可行解：使得目标函数取极值的基本可行解，被称为最优基本可行解，简称最优解。4．线性规划的单纯形解法定义与定理（1）为了使不等式约束条件转变成等式约束条件，需要引进松弛变量（附加变量）。原来问题中的变量被成为原变量（结构变量）。引入的松弛变量个数等于不等式约束条件的个数。（2）最优可行解为个约束条件所形成凸多边形的一个顶点，则最优可行解中至少有个变量为零。为零的变量被称为非基本变量，不为零的变量被称为基本变量。（3）目标函数用非基本变量描述时，变量的系数均大于零。计算规则（1）“指定”基本变量和非基本变量。（2）将基本变量写成非基本变量的函数。（3）将目标函数中的基本变量用上述函数代入，使目标函数中的变量全部为非基本变量。（4）判断目标函数中变量的系数，全部为正数则跳转（7）。39\n研究生课程计算机仿真讲稿（1）从目标函数中找到系数最小的非基本变量。（2）将基本变量函数中的常数项除以该函数中变量的系数。即。找到其中最小的一个。此时指定变量为基本变量，变量为非基本变量。跳转（2）。（3）结束。已知线性规划问题如下（6-14）（6-15）解指定基本变量为非基本变量为由6-15式可得（6-16）（6-17）将式6-16和6-17代入6-14可得（6-18）由于变量的系数最小并小于零，继续进行。从式6-16和6-17中得到。因此，应作为非基本变量，作为基本变量，从式6-17得到（6-19）将式6-19代入式6-16可得（6-20）将式6-19代入式6-18可得（6-21）由于变量的系数最小并小于零，继续进行。从式6-19和6-20中得到。因此，应作为非基本变量，作为基本变量，从式6-20得到（6-22）39\n研究生课程计算机仿真讲稿将式6-22代入式6-19可得（6-23）将式6-22代入式6-21可得（6-24）由于式6-24中的变量的系数均大于零，计算结束。此时非基本变量为，这两个变量取零值，将其代入式6-23、式6-22及式6-24，分别得出；；。此时；。5．单纯形表表中包含一个向量，该向量的维数等于约束条件的个数，向量中记忆基本变量的编号；表中还包含了一个约束条件矩阵的增广矩阵单纯形法计算线性规划，可以转化为矩阵运算，也称为单纯形表。下面举例说明单纯形表的构成以及计算方法。例6-7用单纯形表解下面线性规划问题。（6-25）（6-26）解首先将上述问题化为标准型（6-27）（6-28）以下构造单纯形表，见表6-439\n研究生课程计算机仿真讲稿表6-4单纯形表的构成基本变量向量约束条件增广矩阵32210001241201008540001016604000112-2-300000单纯形表的变换规则如下：（1）判断，是否全部大于零，并找到其中小于零的最小的一个，则应当作为基本变量，所在的位置。（2）计算（注：当时，认为），并找到其中大于零并且数值最小的一个，，所在位置。（3）采用行列式中的初等行变换的方法，让，。这种变换对增广矩阵中各个元素的影响如下：（6-29）（6-30）（4）令。变换结果见表6-5。表6-5单纯形表的构成（变换1）基本变量向量约束条件增广矩阵320100-0.56410010-0.525400010162010000.253-200000.759按照上述原则继续变换可得。39\n研究生课程计算机仿真讲稿表6-6单纯形表的构成（变换2）基本变量向量约束条件增广矩阵3001-200.52110010-0.525000-41282010000.25300020-0.2513按照上述原则继续变换可得。表6-7单纯形表的构成（变换3）基本变量向量约束条件增广矩阵3001-1-0.2500110000.25046000-20.51420100.5-0.125020001.50.125014可得最优解；（p6-182）39研究生课程计算机仿真讲稿clcechooff%A=input('A=')A=[22100012;1201008;40001016;04000112;-2-300000];[m,n]=size(A);m=m-1;n=n-1;%K=input('K=')K=[3456];XM=size(K);while(1)clcAKinputdd;[Amin,l]=min(A(m+1,1:n));if(Amin>=0)disp('end')disp('K=')disp(K)disp('value')disp(A(1:m,n+1)')disp('minofJ')disp(-1*A(m+1,n+1))breakendfori=1:1:mif(A(i,l)~=0),XM(i)=A(i,n+1)/A(i,l);elseXM(i)=inf;endendif(max(XM)<=0),disp('No');break;endw=0;Amin=inf;fori=1:1:mif((XM(i)>0)&(XM(i)<=Amin))w=i;Amin=XM(i);endendM=A(w,l);A(w,1:n+1)=A(w,1:n+1)./M;fori=1:1:m+1if(i==w)i=i+1;endM=A(i,l);A(i,1:n+1)=A(i,1:n+1)-A(w,1:n+1).*M;39\n研究生课程计算机仿真讲稿endK(w)=l;endechoon39研究生课程计算机仿真讲稿三、单变量的寻优方法当目标函数中的自变量为一个单一的变量时，这是寻找使得目标函数取得最优值的问题被称为单变量寻优问题。单变量寻优方法是多变量寻优方法的基础，许多多变量寻优问题往往归结为反复地求解一系列单变量函数的最优值。单变量寻优的方法可以分成两大类：分割法和插值法1．分割法分割法的基本思想为：若目标函数在已知区间存在最小点，区间长度如图6-8所示。00图6-8分割法示意图图6-9等间分割法.按照一定的规则逐渐缩小区间，当区间缩小到，区间长度小于给定的误差时，即，可以近似地认为分割法中使用的分割规则有很多种，以在缩小区间时的计算次数最少为衡量其算法最佳的原则。下面分别介绍等间法和黄金分割法，并比较它们算法。（1）等间分割法等间分割法分割区间的方法是，距离区间的两个顶点三分之二处各找到一点，距离为；距离为。通过比较和的大小来确定新的区间。当则区间为组成39\n研究生课程计算机仿真讲稿当则区间为组成计算复杂度。（2）黄金分割法如果在分割区间时，距离区间的两个顶点处各找到一点，即：=+；=。计算目标函数和的大小来确定新的区间。当时组成。此时再距离处找找，即：=+=+=+==由此可见，在第二次以后计算和时可以少计算一点。算法归纳如下：当则区间为组成，此时=。当则区间为组成，此时=。计算复杂度为(p6_186)functiony=p_186f(x)y=x^3-2*x-5;39研究生课程计算机仿真讲稿clcechooff%a=input('a=');%b=input('b=');%e1=input('e=');a=0;b=5;e1=1e-4;k=1;x1=a+0.618*(b-a);x2=b-0.618*(b-a);fx1=p6_186f(x1);fx2=p6_186f(x2);e=100;while(e>e1)if((fx1>fx2)|(fx1==fx2))b=x1;x1=x2;x2=b-0.618*(b-a);fx1=p6_186f(x1);fx2=p6_186f(x2);elsea=x2;x2=x1;x1=a+0.618*(b-a);fx1=p6_186f(x1);fx2=p6_186f(x2);ende=b-a;k=k+1;endx=(a+b)/2;disp(x)disp(k)disp(p6_186f(x))39\n研究生课程计算机仿真讲稿inputt;x2=fmin('p6_186f',0,5,[1,1e-4])echoon39研究生课程计算机仿真讲稿2．插值法（1）外推极值区间外推极值区间采用加倍步距法。首先在区间内任意找一点，设定初始步长（可以是充分小的正数），令，分别计算和并比较，有两种情况a．如果则令。如果则令。反复执行直到为止。b．如果则令。如果则令。反复执行直到为止。（2）内插求极值知道、和，以及它们对应的目标函数值为、和后，构造二次多项式，构造的方法是，设（6-37）该式应当满足（6-38）由于极值点为到数为零的点，则由6-37式可得求得（6-39）从6-38式解出和代入式6-39后求得；（6-40）39\n研究生课程计算机仿真讲稿（3）迭代计算计算出后，采用6-41式验证是否满足精度要求（6-41）满足则结束。不满足则令等于和中函数值较小的一个；另继续执行。39研究生课程计算机仿真讲稿clcechooffx0=10;h0=0.1;e1=0.0001;m=0.1;x1=x0;k=1;while(1)disp(k)f1=p6_186f(x1);x2=x1+h0;f2=p6_186f(x2);if((f2>f1)|(f2==f1))h=h0;x3=x2+h;f3=p6_186f(x3);while((f2>f3)|(f2==f3))x1=x2;f1=f2;x2=x3;f2=f3;disp('+')h=h+h0;x3=x2+h;f3=p6_186f(x3);endelseh=h0;x3=x1;f3=f1;x1=x2-h0;f2=p6_186f(x2);while((f2>f1)|(f2==f1))x3=x2;f3=f2;h=h+h0;disp('-')x1=x2-h0;f2=p6_186f(x2);endendc1=(f3-f1)/(x3-x1);c2=((f2-f1)/(x2-x1)-c1)/(x2-x3);x4=0.5*(x1+x3-c1/c2);f4=p6_186f(x4);e=(f4-f2)/f2;e=abs(e);%disp(e)if(e=fh)xx=xr;fx=fr;xr=xh;fr=fh;xh=xx;fn=fx;x(imas,:)=xh;fz(imas)=fh;enddisp('-');xs=(1-v)*xh+v*xr;fs=p6_195f(xs);if(fs