《用MATLAB解最优控制问题及应用实例ppt课件.ppt》由会员分享，可在线阅读，更多相关《用MATLAB解最优控制问题及应用实例ppt课件.ppt（128页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第十二章第十二章 用用MATLABMATLAB解最解最优控制问题及应用实例优控制问题及应用实例1.第十二章第十二章 用用MATLABMATLAB解最优解最优控制问题及应用实例控制问题及应用实例12.1 MATLAB12.1 MATLAB工具简介工具简介12.2 12.2 用用MATLABMATLAB解线性二次型最优控制问题解线性二次型最优控制问题12.3 12.3 用用MATLABMATLAB解最优控制问题应用实例解最优控制问题应用实例12.4 12.4 小结小结2 MATLAB是集数值运算、符号运算及图形处理等强大功能于一体的科学计算语言。作为强大的科学计算平台，它几乎能满足所有的计算需求。
2、MATLAB具有编程方便、操作简单、可视化界面、优良的仿真图形环境、丰富的多学科工具箱等优点，尤其是在自动控制领域中MATLAB显示出更为强大的功能。3 最优控制是在一定的约束条件下，从已给定的初始状态出发，确定最优控制作用的函数式，使目标函数为极小或极大。在设计最优控制器的过程中，运用MATLAB最优控制设计工具，会大大减小设计的复杂性。 在前面的几章中，我们已经介绍了一些最优控制方法，在本章中我们将介绍一个最优控制问题的应用实例，讨论如何使用最优控制方法来设计自寻的制导导弹的最优导引律，并采用MATLAB工具实现最优导引律，通过仿真来验证最优导引律的有效性。412.1 MATLAB12.1
3、 MATLAB工具简介工具简介DuCxyBuAxx 1, 系统模型的建立 系统的状态方程为：5 在MATLAB中只需要将各个系数按照常规矩阵的方式输入到工作空间即可 ss(A,B,C,D),;,;,;,;,;,;,;,;,212222111211212222111211212222111211212222111211qpqqppqnqqnnnpnnppnnnnnndddddddddDcccccccccCbbbbbbbbbBaaaaaaaaaA6传递函数的零极点模型为：)()()()()(2121nmpspspszszszsKsG在MATLAB中可以采用如下语句将零极点模型输入到工作空间:;21
4、21nmpppPzzzZKKGainzpk(Z,P,KGain)7传递函数模型在更一般的情况下，可以表示为复数变量s的有理函数形式：nnnnnmmmmasasasasbsbsbsbsG122111121)(8在MATLAB中可以采用如下语句将以上的传递函数模型输入到工作空间:G=tf(num,den);, 1 ;,121121nnmmaaaadenbbbbnum92, 系统模型的转换 把其他形式转换成状态方程模型 G1=ss(G) 把其他形式转换成零极点模型 G1=zpk(G) 把其他形式转换成一般传递函数模型 G1=tf(G)103, 系统稳定性判据 求出系统所有的极点，并观察系统是否有实部
5、大于0的极点。 系统由传递函数 (num,den) 描述 roots(den) 系统由状态方程 (A,B,C,D) 描述 eig(A)114, 系统的可控性与可观测性分析 在MATLAB的控制系统工具箱中提供了ctrbf()函数。该函数可以求出系统的可控阶梯变换，该函数的调用格式为： Ac,Bc,Cc,Dc,Tc,Kc=ctrbf(A,B,C) 在MATLAB的控制系统工具箱中提供了obsvf()函数。该函数可以求出系统的可观测阶梯变换，该函数的调用格式为： Ao,Bo,Co,Do,To,Ko=obsvf(A,B,C)125, 系统的时域分析 对于系统的阶跃响应，控制系统工具箱中给出了一个函数
6、step()来直接求取系统的阶跃响应，该函数的可以有如下格式来调用： y=step(G,t) 对于系统的脉冲响应，控制系统工具箱中给出了一个函数impulse()来直接求取系统的脉冲响应，该函数的可以有如下格式来调用： y=impulse (G,t)136, 系统的复域与频域分析 对于根轨迹的绘制，控制系统工具箱中给出了一个函数rlocus()函数来绘制系统的根轨迹，该函数的可以由如下格式来调用： R=rlocus(G,k)14 对于Nyquist曲线的绘制，控制系统工具箱中给出了一个函数nyquist()函数，该环数可以用来直接求解Nyquist阵列，绘制出Nyquist曲线，该函数的可以由
7、如下格式来调用： rx,ry=nyquist(G,w) 对于Bode图，MATLAB控制工具箱中提供了bode()函数来求取、绘制系统的Bode图，该函数可以由下面的格式来调用 mag,pha=bode(G,w)1512.2 12.2 用用MATLABMATLAB解线性二次型最优控制问题解线性二次型最优控制问题一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中， 为 维状态向量， 为 维控制向量， 为维输出向量。 )()()()()(tUtBtXtAtX)()()(tXtCtY( )X tn( )U tm)(tYl16寻找最优控制，使下面的性能指标最小fttTTffTdttUtRtU
8、tetQtetPeteuJ0)()()()()()(21)()(21)(l其中， 是 对称半正定常数阵， 是 对称半正定阵, 是 对称正定阵。l ll l)(tQ)(tRmmP17 我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：1( )( ) ( )( ) ( )( ) ( )( )( ) ( )( )TTK tK t AtA t K tK t Bt R t B t K tQt18可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成 ，这样可以保证系统状态渐进的趋近于零值，这样可以得
9、出矩阵趋近于常值矩阵 ，且 ，这样上述黎卡提矩阵微分方程可以简化成为：ft)(tK0)(tK1( ) ( )( )( )( ) ( )( )( )( )( )0TTK t A tAt K tK t B t Rt Bt K tQ t19这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。20方法一：方法一：求解代数黎卡提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令 ，则可以写出下面的迭代公式00 11TTTTiiiiiEEEGWGGHEGWQ21AIAIE 1BAIG12BAIQAIBRHTT11)(BAIQ
10、W122如果 收敛于一个常数矩阵，即 ，则可以得出代数黎卡提方程的解为：上面的迭代算法可以用MATLAB来实现：1iii 11112AIAIPiT23%*MATLAB程序*%I=eye(size(A);iA=inv(I-A);E=iA*(I+A);G=2*iA2*B;H=R+B*iA*Q*iA*B;W=Q*iA*B;P0=zeros(size(A);i=0;24while(1),i=i+1; P=E*P0*E-(E*P0*G+W)*inv(G*P0*G+H)*(E*P0*G+W)+Q; if(norm(P-P0)eps),break; else,P0=P; endendP=2*iA*P*iA;
11、我们把这个文件命名为mylq.m，方便我们以后调用来求解代数黎卡提方程。25方法二：方法二： 在MATLAB的控制系统工具箱中提供了求解代数黎卡提方程的函数lqr()，其调用的格式为： K,P,E=lqr(A,B,Q,R) 式中输入矩阵为A,B,Q,R，其中(A,B)为给定的对象状态方程模型，(Q,R)分别为加权矩阵Q和R；返回矩阵K为状态反馈矩阵，P为代数黎卡提方程的解，E为闭环系统的零极点。26这里的求解是建立在MATLAB的控制系统工具箱中给出的一个基于Schur变换的黎卡提方程求解函数are()基础上的，该函数的调用格式为： X=are(M,T,V)27其中， 矩阵满足下列代数黎卡提方
12、程，are是Algebraic Riccati Equation的缩写。对比前面给出的黎卡提方程，可以容易得出VTM,0VXTXXMMXTAMTBBRT1QV28方法三：方法三： 我们也可以采用care()函数对连续时间代数黎卡提 方程求解，其调用方法如下：P,E,K,RR=care(A,B,Q,R,zeros(size(B),eye(size(A) 式中输入矩阵为A,B,Q,R，其中(A,B)为给定的对象状态方程模型，(Q,R)分别为加权矩阵Q和R；返回矩阵P为代数黎卡提方程的解，E为闭环系统的零极点，K为状态反馈矩阵，RR是相应的留数矩阵Res的Frobenius范数（其值为：sqrt(s
13、um(diag(Res*Res)，或者用Norm(Res, fro)计算）。29 采用care函数的优点在于可以设置P的终值条件，例如我们可以在下面的程序中设置P的终值条件为0.2;0.2。 P,E,K,RR=care(A,B,Q,R,0.2;0.2,eye(size(A) 采用lqr()函数不能设置代数黎卡提方程的边界条件。30例例12-112-1 线性系统为： ，其目标函数是：uxx10351006667. 1 10020020050021dtuuxxJTT确定最优控制。31解：解：方法一：A=0 1;-5,-3;B=0;1;Q=500 200;200 100;R=1.6667;mylqK
14、=inv(R)*B*PPE32运行结果：K = 13.0276 6.7496P = 67.9406 21.7131 21.7131 11.2495E = -0.1111 0.2222 -1.1111 -0.777833方法二:A=0 1;-5,-3;B=0;1;Q=500 200;200 100;R=1.6667;K,P,E=lqr(A,B,Q,R)34运行结果：K = 13.0276 6.7496P = 67.9406 21.7131 21.7131 11.2495E =-7.2698 -2.479835方法三：A=0 1;-5,-3;B=0;1;Q=500 200;200 100;R=1.
15、6667;P,E,K,RR=care(A,B,Q,R,zeros(size(B),eye(size(A) 36运行结果：P = 67.9406 21.7131 21.7131 11.2495E = -7.2698 -2.4798K =13.0276 6.7496RR = 2.8458e-01537 以上的三种方法的运行结果相同。我们可以得到，最优控制变量与状态变量之间的关系： 在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：)(6.7496)(13.0276)(21*txtxtu38%*MATLAB程序*%figure(pos,50,50,200,150
16、,color,w);axes(pos,0.15,0.14,0.72,0.72)ap=A-B*K;bp=B;C=1,0;D=0;39ap,bp,cp,dp=augstate(ap,bp,C,D);cp=cp;-K;dp=dp;0;G=ss(ap,bp,cp,dp);y,t,x=step(G);plotyy(t,y(:,2:3),t,y(:,4)ax,h1,h2=plotyy(t,y(:,2:3),t,y(:,4);axis(ax(1),0 2.5 0 0.1),axis(ax(2),0 2.5 -1 0)40运行结果：运行结果： 图图12-1 12-1 最优控制曲线与最优状态曲线最优控制曲线与最
17、优状态曲线41 该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量输出，y(4)是系统控制变量输出。用plotyy函数进行双坐标显示，并设置相应的坐标范围。42 以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠。第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法。但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我们只能使用第三种方法来求解，例如设置P的终值条件为0.2;0.2。43程序如下：%*MA
18、TLAB程序*%A=0 1;-5,-3;B=0;1;Q=500 200;200 100;R=1.6667;P,E,K,RR=care(A,B,Q,R,0.2;0.2,eye(size(A)44运行结果：P =67.7233 21.5685 21.5685 11.0961E =-7.3052 -2.4723K =13.0608 6.7775RR =1.2847e-014最优控制变量与状态变量之间的关系：)(6.7775)(13.0608)(21*txtxtu45例例12-212-2 无人飞行器的最优高度控制，飞行器的控制方程如下 )(2/100)()()(2/100100010)()()(tut
19、hththththth 46 是飞行器的高度； 是油门输入；设计控制律使得如下指标最小)(th)(tudttuththththththtutxJ02)(2)()()(000000001)(),(),(21)(),( 初始状态 。绘制系统状态与控制输入，对如下给定的 矩阵进行仿真分析.Tththth0 , 0 ,10)(),()( ，,QR47a).b).c).d).100000 ,R2000Q2000R,000000001Q2,0000000010RQ10001000 ,2000QR48解：解：线性二次型最优控制指标如下： 其中Q和R分别是对状态变量和控制量的加权矩阵， 线性二次型最优控制器设
20、计如下：1）、Q=diag（1，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为 k1=0.7071 2.0772 2.0510， u（t）=k1*x（t）； 所画状态响应曲线及控制输入响应曲线如下图12-2所示：4950图图12-2 12-2 状态响应曲线及控制输入响应曲线状态响应曲线及控制输入响应曲线512）、Q=diag（1，0，0），R=2000时，由MATLAB求得最优状态反馈矩阵为k2=0.0224 0.2517 0.4166， u（t）= k2*x（t）； 所画状态响应曲线及控制输入响应曲线如下图12-3所示：5253图图12-3 12-3 状态响应曲线及控制输入响应曲线
21、状态响应曲线及控制输入响应曲线543）、Q=diag（10，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为 k3=2.2361 4.3892 3.3077，u（t）= k3*x（t）； 所画状态响应曲线及控制输入响应曲线如下图12-4所示：5556图图12-4 12-4 状态响应曲线及控制输入响应曲线状态响应曲线及控制输入响应曲线574）、Q=diag（1，100，0），R=2时，由MATLAB求 得最优状态反馈矩阵为 k4=0.7071 7.6112 4.6076， u（t）= k4*x（t）；所画状态响应曲线及控制输入响应曲线如下图12-5所示：5859图图12-5 12-5 状
22、态响应曲线及控制输入响应曲线状态响应曲线及控制输入响应曲线60由1），2），3），4）可分析如下：图12-3与图12-2相比，当Q不变，R增大时，各相应曲线达到稳态所需时间增长，即响应变慢；但波动幅值变小，反馈矩阵变小；图12-4与图12-2和图12-3相比，当Q对角线上第1个元素增大时，各相应曲线达到稳态所需时间变短，即响应快；但波动幅值变大，反馈矩阵增大；61由图12-5可知，当Q对角线上第2个元素增大时，状态x1，x2曲线达到稳态所需时间较长，即响应较慢，平缓的趋于零；状态x3，控制输入u达到稳态所需时间短，即响应快；状态x2，x3波动幅值较小，比图12-2和图12-4小，比图12-3稍
23、大，控制输入u波动幅值比图12-2和图12-4小，比图12-3大；反馈矩阵最大。62 综上所述可得结论：Q=diag（1，0，0），R=2时，系统各方面响应较好。 矩阵Q变大时，反馈矩阵变大；当Q的对角线上第1个元素变大时，各曲线波动幅值变大，达到稳态所需时间变短； 63 当Q的对角线上第2个元素变大时，各曲线波动幅值变小；达到稳态所需时间，状态x1，x2增长，状态x3，控制输入u变短； 当R变大时，反馈矩阵变小；各曲线波动幅值变小；达到稳态所需时间变长。 所以根据实际的系统允许，我们应该适当选择Q和R。64%*MATLAB程序*%a=0 1 0;0 0 1;0 0 -1/2;b=0;0;1/
24、2;c=1 0 0;0 1 0;0 0 1;d=0;0;0;figure(1)q=1 0 0;0 0 0;0 0 0;r=2;k,p,e=lqr(a,b,q,r)x0=10;0;0;a1=a-b*k;y,x=initial(a1,b,c,d,x0,20);65n=length(x(:,3);T=0:20/n:20-20/n;plot(T,x(:,1),black,T,x(:,2),red,T,x(:,3),green);xlabel(time-s);ylabel(response);title(图(1.a) Q=diag(1,0,0),R=2时状态响应曲线)grid,hold onfor j=
25、1:n u(j,:)=-k*(x(j,:);end66figure(2)plot(T,u);xlabel(time-s);ylabel(response);title(图(1.b) Q=diag(1,0,0),R=2时控制输入u的响应曲线)grid,hold on%*67figure(3)qa=1 0 0;0 0 0;0 0 0;ra=2000;ka,pa,ea=lqr(a,b,qa,ra)x0=10;0;0;aa1=a-b*ka;ya,xa=initial(aa1,b,c,d,x0,60);na=length(xa(:,3);68Ta=0:60/na:60-60/na;plot(Ta,xa(
26、:,1),black,Ta,xa(:,2),red,Ta,xa(:,3),green);xlabel(time-s);ylabel(response);title(图(2.a) Q=diag(1,0,0),R=2000时状态响应曲线)grid,hold onfor j=1:na ua(j,:)=-ka*(xa(j,:);end69figure(4)plot(Ta,ua);xlabel(time-s);ylabel(response);title(图(2.b) Q=diag(1,0,0),R=2000时控制输入u的响应曲线)grid,hold on%*70figure(5)qb=10 0 0;0
27、 0 0;0 0 0;rb=2;kb,pb,eb=lqr(a,b,qb,rb)x0=10;0;0;ab1=a-b*kb;yb,xb=initial(ab1,b,c,d,x0,20);nb=length(xb(:,3);71Tb=0:20/nb:20-20/nb;plot(Tb,xb(:,1),black,Tb,xb(:,2),red,Tb,xb(:,3),green);xlabel(time-s);ylabel(response);title(图(3.a) Q=diag(10,0,0),R=2时状态响应曲线)grid,hold onfor j=1:nb ub(j,:)=-kb*(xb(j,:)
28、;end72figure(6)plot(Tb,ub);xlabel(time-s);ylabel(response);title(图(3.b) Q=diag(10,0,0),R=2时控制输入u的响应曲线)grid,hold on%*73figure(7)qc=1 0 0;0 100 0;0 0 0;rc=2;kc,pc,ec=lqr(a,b,qc,rc)x0=10;0;0;ac1=a-b*kc;yc,xc=initial(ac1,b,c,d,x0,50);nc=length(xc(:,3);74Tc=0:50/nc:50-50/nc;plot(Tc,xc(:,1),black,Tc,xc(:,
29、2),red,Tc,xc(:,3),green);xlabel(time-s);ylabel(response);title(图(4.a) Q=diag(1,100,0),R=2时状态响应曲线)grid,hold onfor j=1:nc uc(j,:)=-kc*(xc(j,:);end75figure(8)plot(Tc,uc);xlabel(time-s);ylabel(response);title(图(4.b) Q=diag(1,100,0),R=2时控制输入u的响应曲线)grid,hold on7612.3 12.3 用用MATLABMATLAB解最优控制问题应用实例解最优控制问题应
30、用实例12.3.1 12.3.1 导弹运动状态方程的建立导弹运动状态方程的建立12.3.2 12.3.2 最优导引律的求解与仿真验证最优导引律的求解与仿真验证77在现有的自寻的导弹中，大都采用比例导引法。假设导弹和目标在同一平面内运动，按比例导引制导律，假设导弹的速度向量的旋转角速度 垂直于瞬时的弹目视线，并且正比于导弹与目标之间的视线角速率 ，假设目标的法向加速度为零，那么可得： (12-1)q qN78其中， 为导弹的速度与基准方向的夹角， 为导弹与目标连线与基准方向的夹角，称为视线角， 是视线角速率， 是比例常数，称为导航比，通常为36。比例导引的实质是使导弹向着 减小的方向运动，抑制视
31、线旋转，也就是使导弹的相对速度对准目标，保证导弹向着前置碰撞点飞行。qq Nq 79 比例导引法是经典的导引方法。下面我们从最优控制理论的观点来研究自寻的导弹的最优导引规律问题。8012.3.1 12.3.1 导弹运动状态方程的建立导弹运动状态方程的建立 导弹与目标的运动关系是非线性的，如果把导弹与目标的运动方程相对于理想弹道线性化，可得导弹运动的线性状态方程. 81假设导弹和目标在同一平面内运动，如图12-6所示。选 为固定坐标。导弹速度向量 与 轴成 角，目标速度向量 为与 轴成 角。导弹与目标的连线 与轴 成 角。假定导弹以尾追的方式攻击目标。坐标轴 和 的方向可以任意选择，使 和 都比
32、较小。再假定导弹和目标均匀速飞行，也就是说 和 均为恒值。使用相对坐标状态变量，设 为导弹与目标在 轴方向上的距离偏差， 为导弹与目标在 轴方向上的距离偏差，即 (12-2)oxyMVoyTVoyTMToyqoxoyT,qMVTVxoxyoyMTMTyyyxxx82xxyyMxTxMyTyMVTV0qTRTMH图图12-6 12-6 导弹和目标运动几何关系图导弹和目标运动几何关系图83假定 和 比较小，因此，则1cos, 1cos,sin,sinTTTT将上式对t求导，并根据导弹和目标的关系（如图12-6所示）可得 (12-3) coscossinsinMTTMTMTTMTVVyyyVVxxx
33、MTMTMTTMTVVyyyVVxxx (12-4)84以 表示 ， 表示 （即 ），则 (12-5) (12-6)1xx2xx 1x 21xx MTTVVxx2式中 表示目标的横向加速度， 表示导弹横向加速度，分别以 和 表示，那么 (12-7)TTVMVTaMaMTaax285导弹的横向加速度 为一控制量。一般将控制信号加给舵机，舵面偏转后产生弹体攻角 ，而后产生横向加速度 。如果忽略舵机和弹体的惯性，而且假设控制量的单位与加速度单位相同，则可用控制量 来表示 ，也就是令 (12-8) 所以(12-7)式为： (12-9)MaMauMaMauuaxT286这样可得导弹运动状态方程为： (1
34、2-10) (12-11)21xx Taux287可写成矩阵的形式： (12-12)式中， ， ， ， 。 (12-13)如果不考虑目标的机动，即 ，则在这种情况下，式(12-12)变成： (12-14)TDaBuAXX21xxX0010A10B10D0TaBuAXX88下面来考虑(12-4)式，该式可写成 (12-15) 其中 表示导弹相对目标的接近速度。由于 的值都比较小， 可近似表示导弹与目标之间的距离。设 为导弹与目标的遭遇时刻(即导弹与目标相碰撞或两者之间的距离为最短的时刻)，则在某一瞬时 ，导弹与目标的距离 可近似用下式表示： (12-16)(TMVVyCTMVVVT,和qyftt
35、y)()()(ttVttVVtyfCfTM89又考虑到对于导弹制导来说，最基本的要求是脱靶量越小越好，因此，应该选择最优控制量 ，使得下面的指标函数为最小。 (12-17)u22)()()()(fMfTfMfTtytytxtxJ90然而，当要求一个反馈形式的控制时，按上式列出的问题很难求解。所以我们以 时刻，即 时 的值 作为脱靶量，要求 值越小越好。另外，由于舵偏角受到限制，导弹结构能够承受的最大载荷也受到限制，所以控制信号 也应该受到限制。 ftt 0)()(ffCfttVty)(1ftx)(1ftxu91因此，我们选择以下形式的二次型指标函数： (12-18)式中 ， 。 (12-19)
36、即 (12-20)给定初始条件 ，应用最优控制理论，可以求出使 为最小的 。dtRuuQXXtCXtXJfttTTffT021212100ccC0000QdtRutCXtXJfttffT022121 0tXJu92 由于系统是线性的，指标函数是二次型的，因此，求最优控制规律就可以认为是一个求解线性二次型的过程。 对于线性二次型问题，可采用变分法、极小值原理、动态规划或其他方法求得最优控制93 (12-21)式中 满足下列黎卡提矩阵微分方程 (12-22) 的终端条件为 (12-23)因此求解线性二次型问题的关键是求解黎卡提矩阵微分方程。 KXBRtuT1*QKBKBRKAKAKTT1KKCtK
37、f9412.3.2 12.3.2 最优导引律的求解与仿真验证最优导引律的求解与仿真验证当不考虑弹体惯性时，而且假定目标不机动，即，导弹运动状态方程为 (12-24)BuAXX95指标函数为 (12-25)式中 ， , ， 。 dtRuuQXXtCXtXJfttTTffT021210010A10B2100ccC0000Q96给出 时刻， 的初值 ，采用极小值原理可求得最优控制 为 (12-26)0tt 21xx 和 0201txtx和 tu* 2321 21 211212*34211 22231312fffffffcc ttcc ttc ttxcc ttxRRu tc ttc ttcc ttRR
38、RR 97在指标函数中，如不考虑导弹的相对运动速度 项，则可令 。 变成 (12-27)以 除上式的分子和分母，得 (12-28)2x02c tu* RttcRxttcxttctufff313122111*1c 31221*333ttcRxttxtttufff98为了使脱靶量为最小，应选取 ，则 (12-29)根据图12-6可得 1c ttxttxtuff221*3ttVxyxtgqfc1199当 比较小时， ，则 (12-30) (12-31)qqtgq ttVxqfc1ttxttxVttVxttxqffcfcf2212111100将上式代入(12-29)式，可得 (12-32)在上式中，
39、的单位是加速度的单位 。把 与导弹速度向量 的旋转角速度 联系起来，则有 (12-33) qVtuc3*u2秒米DVqVVVuMcM3u101从(12-32)和(12-33)式可以看出，当不考虑弹体惯性时，最优导引规律就是比例导引，其导航比为 。这证明了比例导引是一种很好的导引方法。最优导引规律的形成可用图12-7来表示。McVV3102 下面将对最优导引律进行MATLAB仿真，并给出源代码和仿真结果。 cV3s1s1ttVfc121ttVfcq*u1x2x图图12-7 12-7 最优导引方框图最优导引方框图103xxyyMxTxMyTyMVTV0HqTTMRTaMaMTLHE图图12-8 1
40、2-8 最优导引攻击几何平面最优导引攻击几何平面104最优导引攻击几何关系如图12-8所示，在这里讨论的目标和导弹均认为是二维拦截几何平面上的质点，分别以速度 和 运动。导弹的初始位置为相对坐标系的参考点，导弹初始速度矢量指向目标的初始位置， 为导弹的指令（垂直于视线）。TVMVMa105其中： (12-34) (12-35) (12-36) 为目标速度在 轴上的分解， 是目标的角度。导弹和目标之间的接近速度为： (12-37)TTTVaTTTyVVcosTTTxVVsinTcTMVR ,TxTyVV, x y106目标的速度分量可由其位置变化得到： (12-38)同样地，我们可以得到导弹的位
41、置和速度的微分方程： ， (12-39) ， (12-40)TxTxTyTyVRVR,MxMxaVMyMyaVMxMxVRMyMyVR107上面几式中的下标x，y分别表示在x和y轴上的分量。 是导弹在地球坐标系的加速度分量。为了得到导弹的加速度分量，我们必须得到弹目的相对位移： (12-41) (12-42)MyMxaa,MxTxTMxRRRMyTyTMyRRR108TMyTMxRRq1tanMxTxTMxVVVMyTyTMyVVV从图12-8中，根据三角关系我们可以得到视线角： (12-43)如果定义地球坐标系的速度分量为： (12-44) (12-45)109我们可以根据视线角的公式求导后
42、得到视线角速率： (12-46) (12-47)所以我们不难得出弹目的接近速度为： (12-48)TMTMyTMxTMxTMyRVRVRq22122)(TMyTMxTMRRRTMTMyTMyTMxTMxTMcRVRVRRV)(110根据最优导引制导律： (12-49)可得到导弹的加速的分量为： (12-50) (12-51) (12-52)qVVMc3)(tan1MyMxVVqVaMMxcosqVaMMysin111 以上列出了两维的最优导引制导的必要方程，但是使用最优导引制导的导弹并不是直接向着目标发射的，而是向着一个能够导引导弹命中目标的方向发射，考虑了视线角之后可以得到导弹的指向角L。从
43、图12-8中我们可以看出，如果导弹进入了碰撞三角区（如果目标和导弹同时保持匀速直线运动，导弹必定会命中目标），这时利用正弦公式可以得到指向角的表达式：MTTVqVL)sin(sin1(12-53)112 但是实际上导弹不可能能确切地在碰撞三角区发射，所以不能精确地得到拦截点。因为我们不知道目标将会如何机动，所以拦截点位置只能大概地估计。事实上，这也是需要导航系统的原因！初始时刻导弹偏离碰撞三角的角度称之为指向角误差（Head-Error）。考虑了导弹初始时刻的指向角和指向角误差之后，导弹的初始速度分量可以表示为：)cos()0(HeadErrorLqVVMMy)sin()0(HeadError
44、LqVVMMx(12-54)(12-55)113使用MATLAB编程，具体代码如下：%*MATLAB程序*%最优制导律仿真，初始化系统的参数clear all; %清除所有内存变量global SignVc;pi=3.14159265;Vm=1000;Vt=500;%导弹和目标的速度HeadError=0; %指向角误差114ThetaT=pi; %目标的速度方向Rmx=0;Rmy=0; %导弹的位置Rtx=5000;Rty=10000;%目标的位置At=0; %目标法向加速度Vtx=Vt*sin(ThetaT);%目标的速度分量Vty=Vt*cos(ThetaT);Rtmx=Rtx-Rmx;
45、 %弹目相对距离Rtmy=Rty-Rmy;AmMax=15*9.8; %导弹的最大机动能力为15GRtm=sqrt(Rtmx2+Rtmy2);115SightAngle=atan(Rtmx/Rtmy); %视线角LeadAngle=asin(Vt*sin(SightAngle-ThetaT)/Vm); %指向角Vmx=Vm*sin(SightAngle-LeadAngle+HeadError); %导弹的速度分量Vmy=Vm*cos(SightAngle-LeadAngle+HeadError);Vtmx=Vtx-Vmx;Vtmy=Vty-Vmy;%弹目的相对运动速度116Vc=-(Rtmx*
46、Vtmx+Rtmy*Vtmy)/Rtm;SignVc=sign(Vc); %Vc的符号Time=0;TimeStep=0.1; %时间和时间步长file=fopen(output.txt,w);%将数据写入文件%循环while(1)%Vc改变符号仿真结束if(sign(Vc) = SignVc) break;else117if(RtmAmMax)Am=AmMax;end%加速度分量Amx=Am*cos(SightAngle);Amy=-Am*sin(SightAngle);Time=Time+TimeStep;%目标位置119Rtx=Rtx+TimeStep*Vtx;Rty=Rty+TimeS
47、tep*Vty;dThetaT=At/Vt;ThetaT=ThetaT+TimeStep*dThetaT;Vtx=Vt*sin(ThetaT);%目标的速度分量Vty=Vt*cos(ThetaT);%导弹位置Rmx=Rmx+TimeStep*Vmx;Rmy=Rmy+TimeStep*Vmy;%导弹速度120Vmx=Vmx+TimeStep*Amx;Vmy=Vmy+TimeStep*Amy;Vm=sqrt(Vmx2+Vmy2);%弹目相对位移Rtmx=Rtx-Rmx;Rtmy=Rty-Rmy;%上一步的脱靶量Rtm0=Rtm;Rtm=sqrt(Rtmx2+Rtmy2);SightAngle=at
48、an(Rtmx/Rtmy); %视线角%弹目相对速度121Vtmx=Vtx-Vmx;Vtmy=Vty-Vmy;Vc=-(Rtmx*Vtmx+Rtmy*Vtmy)/Rtm;%数据写如文件fprintf(file,%f %f %f %f %f %f %f n,Time,Rmx,Rmy,Rtx,Rty,sqrt(Amx2+Amy2), Rtm);endendstatus=fclose(file);122仿真结果：图图12-9 12-9 指向误差，目标不机动的攻击情况指向误差，目标不机动的攻击情况0123图图12-1012-10 指向误差，目标不机动时导弹的加速度指向误差，目标不机动时导弹的加速度01
49、24图图12-1112-11 指向误差，目标不机动时攻击情况指向误差，目标不机动时攻击情况20125图图12-1212-12 指向误差，目标不机动时导弹加速度指向误差，目标不机动时导弹加速度2012612.4 12.4 小结小结 本章首先简单介绍了一下MATLAB控制系统分析工具的使用，接着介绍了如何用MATLAB工具设计线性二次型最优控制器。线性二次型最优控制可以归结为求解黎卡提（Riccati）方程。最后，在此基础上，我们给出了两个运用MATLAB求解线性二次型最优控制问题的例子。127 此外，本章还讨论了如何使用最优控制方法来设计自寻的制导导弹的最优导引律。通过最优二次型控制问题的求解，我们发现：当不考虑弹体惯性时，应用最优控制理论所得到的最优导引规律与目前广泛采用的比例导引法相一致。因此，从现代控制理论的观点来看，比例导引法是一种很好的导引方法。最后，我们还给出了最优导引律的仿真验证程序和仿真曲线。128