第十二章用MATLAB解最优控制问题及应用实例

1、精选ppt1第十二章第十二章 用用MATLABMATLAB解最解最优控制问题及应用实例优控制问题及应用实例精选ppt2第十二章第十二章 用用MATLABMATLAB解最优解最优控制问题及应用实例控制问题及应用实例12.1 MATLAB12.1 MATLAB工具简介工具简介12.2 12.2 用用MATLABMATLAB解线性二次型最优控制问题解线性二次型最优控制问题12.3 12.3 用用MATLABMATLAB解最优控制问题应用实例解最优控制问题应用实例12.4 12.4 小结小结精选ppt3 MATLAB是集数值运算、符号运算及图形处理等强大功能于一体的科学计算语言。作为强大的科学计算平台

2、，它几乎能满足所有的计算需求。MATLAB具有编程方便、操作简单、可视化界面、优良的仿真图形环境、丰富的多学科工具箱等优点，尤其是在自动控制领域中MATLAB显示出更为强大的功能。精选ppt4 最优控制是在一定的约束条件下，从已给定的初始状态出发，确定最优控制作用的函数式，使目标函数为极小或极大。在设计最优控制器的过程中，运用MATLAB最优控制设计工具，会大大减小设计的复杂性。 在前面的几章中，我们已经介绍了一些最优控制方法，在本章中我们将介绍一个最优控制问题的应用实例，讨论如何使用最优控制方法来设计自寻的制导导弹的最优导引律，并采用MATLAB工具实现最优导引律，通过仿真来验证最优导引律的

3、有效性。精选ppt512.1 MATLAB12.1 MATLAB工具简介工具简介DuCxyBuAxx 1, 系统模型的建立 系统的状态方程为：精选ppt6 在MATLAB中只需要将各个系数按照常规矩阵的方式输入到工作空间即可 ss(A,B,C,D),;,;,;,;,;,;,;,;,212222111211212222111211212222111211212222111211qpqqppqnqqnnnpnnppnnnnnndddddddddDcccccccccCbbbbbbbbbBaaaaaaaaaA精选ppt7传递函数的零极点模型为：)()()()()(2121nmpspspszszszsK

4、sG在MATLAB中可以采用如下语句将零极点模型输入到工作空间:;2121nmpppPzzzZKKGainzpk(Z,P,KGain)精选ppt8传递函数模型在更一般的情况下，可以表示为复数变量s的有理函数形式：nnnnnmmmmasasasasbsbsbsbsG122111121)(精选ppt9在MATLAB中可以采用如下语句将以上的传递函数模型输入到工作空间:G=tf(num,den);, 1 ;,121121nnmmaaaadenbbbbnum精选ppt102, 系统模型的转换 把其他形式转换成状态方程模型 G1=ss(G) 把其他形式转换成零极点模型 G1=zpk(G) 把其他形式转换

5、成一般传递函数模型 G1=tf(G)精选ppt113, 系统稳定性判据 求出系统所有的极点，并观察系统是否有实部大于0的极点。 系统由传递函数 (num,den) 描述 roots(den) 系统由状态方程 (A,B,C,D) 描述 eig(A)精选ppt124, 系统的可控性与可观测性分析 在MATLAB的控制系统工具箱中提供了ctrbf()函数。该函数可以求出系统的可控阶梯变换，该函数的调用格式为： Ac,Bc,Cc,Dc,Tc,Kc=ctrbf(A,B,C) 在MATLAB的控制系统工具箱中提供了obsvf()函数。该函数可以求出系统的可观测阶梯变换，该函数的调用格式为： Ao,Bo,C

6、o,Do,To,Ko=obsvf(A,B,C)精选ppt135, 系统的时域分析 对于系统的阶跃响应，控制系统工具箱中给出了一个函数step()来直接求取系统的阶跃响应，该函数的可以有如下格式来调用： y=step(G,t) 对于系统的脉冲响应，控制系统工具箱中给出了一个函数impulse()来直接求取系统的脉冲响应，该函数的可以有如下格式来调用： y=impulse (G,t)精选ppt146, 系统的复域与频域分析 对于根轨迹的绘制，控制系统工具箱中给出了一个函数rlocus()函数来绘制系统的根轨迹，该函数的可以由如下格式来调用： R=rlocus(G,k)精选ppt15 对于Nyqui

7、st曲线的绘制，控制系统工具箱中给出了一个函数nyquist()函数，该环数可以用来直接求解Nyquist阵列，绘制出Nyquist曲线，该函数的可以由如下格式来调用： rx,ry=nyquist(G,w) 对于Bode图，MATLAB控制工具箱中提供了bode()函数来求取、绘制系统的Bode图，该函数可以由下面的格式来调用 mag,pha=bode(G,w)精选ppt1612.2 12.2 用用MATLABMATLAB解线性二次型最优控制问题解线性二次型最优控制问题一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中， 为 维状态向量， 为 维控制向量， 为维输出向量。 )()()

8、()()(tUtBtXtAtX)()()(tXtCtY( )X tn( )U tm)(tYl精选ppt17寻找最优控制，使下面的性能指标最小fttTTffTdttUtRtUtetQtetPeteuJ0)()()()()()(21)()(21)(l其中， 是 对称半正定常数阵， 是 对称半正定阵, 是 对称正定阵。l ll l)(tQ)(tRmmP精选ppt18 我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：1( )( ) ( )( ) ( )( ) ( )( )( ) ( )( )TTK tK t AtA t K tK t Bt R t B t

9、 K tQt精选ppt19可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成 ，这样可以保证系统状态渐进的趋近于零值，这样可以得出矩阵趋近于常值矩阵 ，且 ，这样上述黎卡提矩阵微分方程可以简化成为：ft)(tK0)(tK1( ) ( )( )( )( ) ( )( )( )( )( )0TTK t A tAt K tK t B t Rt Bt K tQ t精选ppt20这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。精选ppt21方法一：方法一：求解代数黎卡

10、提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令 ，则可以写出下面的迭代公式00 11TTTTiiiiiEEEGWGGHEGWQ精选ppt22AIAIE 1BAIG12BAIQAIBRHTT11)(BAIQW1精选ppt23如果 收敛于一个常数矩阵，即 ，则可以得出代数黎卡提方程的解为：上面的迭代算法可以用MATLAB来实现：1iii 11112AIAIPiT精选ppt24%*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

11、;精选ppt25while(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;我们把这个文件命名为mylq.m，方便我们以后调用来求解代数黎卡提方程。精选ppt26方法二：方法二： 在MATLAB的控制系统工具箱中提供了求解代数黎卡提方程的函数lqr()，其调用的格式为： K,P,E=lqr(A,B,Q,R) 式中输入矩阵为A,B,Q,R，其中(A,B)为给定的对象状态方程模型，(Q,R)分别为加权矩阵Q和R；返回矩阵K为

12、状态反馈矩阵，P为代数黎卡提方程的解，E为闭环系统的零极点。精选ppt27这里的求解是建立在MATLAB的控制系统工具箱中给出的一个基于Schur变换的黎卡提方程求解函数are()基础上的，该函数的调用格式为： X=are(M,T,V)精选ppt28其中， 矩阵满足下列代数黎卡提方程，are是Algebraic Riccati Equation的缩写。对比前面给出的黎卡提方程，可以容易得出VTM,0VXTXXMMXTAMTBBRT1QV精选ppt29方法三：方法三： 我们也可以采用care()函数对连续时间代数黎卡提 方程求解，其调用方法如下：P,E,K,RR=care(A,B,Q,R,zer

13、os(size(B),eye(size(A) 式中输入矩阵为A,B,Q,R，其中(A,B)为给定的对象状态方程模型，(Q,R)分别为加权矩阵Q和R；返回矩阵P为代数黎卡提方程的解，E为闭环系统的零极点，K为状态反馈矩阵，RR是相应的留数矩阵Res的Frobenius范数（其值为：sqrt(sum(diag(Res*Res)，或者用Norm(Res, fro)计算）。精选ppt30 采用care函数的优点在于可以设置P的终值条件，例如我们可以在下面的程序中设置P的终值条件为0.2;0.2。 P,E,K,RR=care(A,B,Q,R,0.2;0.2,eye(size(A) 采用lqr()函数不能

14、设置代数黎卡提方程的边界条件。精选ppt31例例12-112-1 线性系统为： ，其目标函数是：uxx10351006667. 1 10020020050021dtuuxxJTT确定最优控制。精选ppt32解：解：方法一：A=0 1;-5,-3;B=0;1;Q=500 200;200 100;R=1.6667;mylqK=inv(R)*B*PPE精选ppt33运行结果：K = 13.0276 6.7496P = 67.9406 21.7131 21.7131 11.2495E = -0.1111 0.2222 -1.1111 -0.7778精选ppt34方法二:A=0 1;-5,-3;B=0;

15、1;Q=500 200;200 100;R=1.6667;K,P,E=lqr(A,B,Q,R)精选ppt35运行结果：K = 13.0276 6.7496P = 67.9406 21.7131 21.7131 11.2495E =-7.2698 -2.4798精选ppt36方法三：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,zeros(size(B),eye(size(A) 精选ppt37运行结果：P = 67.9406 21.7131 21.7131 11.2495E = -7.2698 -2.4798

16、K =13.0276 6.7496RR = 2.8458e-015精选ppt38 以上的三种方法的运行结果相同。我们可以得到，最优控制变量与状态变量之间的关系： 在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：)(6.7496)(13.0276)(21*txtxtu精选ppt39%*MATLAB程序*%figure(pos,50,50,200,150,color,w);axes(pos,0.15,0.14,0.72,0.72)ap=A-B*K;bp=B;C=1,0;D=0;精选ppt40ap,bp,cp,dp=augstate(ap,bp,C,D);cp

17、=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)精选ppt41运行结果：运行结果： 图图12-1 12-1 最优控制曲线与最优状态曲线最优控制曲线与最优状态曲线精选ppt42 该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量

18、输出，y(4)是系统控制变量输出。用plotyy函数进行双坐标显示，并设置相应的坐标范围。精选ppt43 以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠。第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法。但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我们只能使用第三种方法来求解，例如设置P的终值条件为0.2;0.2。精选ppt44程序如下：%*MATLAB程序*%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;

19、0.2,eye(size(A)精选ppt45运行结果：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*txtxtu精选ppt46例例12-212-2 无人飞行器的最优高度控制，飞行器的控制方程如下 )(2/100)()()(2/100100010)()()(tuthththththth 精选ppt47 是飞行器的高度； 是油门输入；设计控制律使得如下指标最小)(th)(tudttuththththt

20、hthtutxJ02)(2)()()(000000001)(),(),(21)(),( 初始状态 。绘制系统状态与控制输入，对如下给定的 矩阵进行仿真分析.Tththth0 , 0 ,10)(),()( ，,QR精选ppt48a).b).c).d).100000 ,R2000Q2000R,000000001Q2,0000000010RQ10001000 ,2000QR精选ppt49解：解：线性二次型最优控制指标如下： 其中Q和R分别是对状态变量和控制量的加权矩阵， 线性二次型最优控制器设计如下：1）、Q=diag（1，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为 k1=0.7071

21、 2.0772 2.0510， u（t）=k1*x（t）； 所画状态响应曲线及控制输入响应曲线如下图12-2所示：精选ppt50精选ppt51图图12-2 12-2 状态响应曲线及控制输入响应曲线状态响应曲线及控制输入响应曲线精选ppt522）、Q=diag（1，0，0），R=2000时，由MATLAB求得最优状态反馈矩阵为k2=0.0224 0.2517 0.4166， u（t）= k2*x（t）； 所画状态响应曲线及控制输入响应曲线如下图12-3所示：精选ppt53精选ppt54图图12-3 12-3 状态响应曲线及控制输入响应曲线状态响应曲线及控制输入响应曲线精选ppt553）、Q=di

22、ag（10，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为 k3=2.2361 4.3892 3.3077，u（t）= k3*x（t）； 所画状态响应曲线及控制输入响应曲线如下图12-4所示：精选ppt56精选ppt57图图12-4 12-4 状态响应曲线及控制输入响应曲线状态响应曲线及控制输入响应曲线精选ppt584）、Q=diag（1，100，0），R=2时，由MATLAB求 得最优状态反馈矩阵为 k4=0.7071 7.6112 4.6076， u（t）= k4*x（t）；所画状态响应曲线及控制输入响应曲线如下图12-5所示：精选ppt59精选ppt60图图12-5 12-5

23、状态响应曲线及控制输入响应曲线状态响应曲线及控制输入响应曲线精选ppt61由1），2），3），4）可分析如下：图12-3与图12-2相比，当Q不变，R增大时，各相应曲线达到稳态所需时间增长，即响应变慢；但波动幅值变小，反馈矩阵变小；图12-4与图12-2和图12-3相比，当Q对角线上第1个元素增大时，各相应曲线达到稳态所需时间变短，即响应快；但波动幅值变大，反馈矩阵增大；精选ppt62由图12-5可知，当Q对角线上第2个元素增大时，状态x1，x2曲线达到稳态所需时间较长，即响应较慢，平缓的趋于零；状态x3，控制输入u达到稳态所需时间短，即响应快；状态x2，x3波动幅值较小，比图12-2和图12

24、-4小，比图12-3稍大，控制输入u波动幅值比图12-2和图12-4小，比图12-3大；反馈矩阵最大。精选ppt63 综上所述可得结论：Q=diag（1，0，0），R=2时，系统各方面响应较好。 矩阵Q变大时，反馈矩阵变大；当Q的对角线上第1个元素变大时，各曲线波动幅值变大，达到稳态所需时间变短； 精选ppt64 当Q的对角线上第2个元素变大时，各曲线波动幅值变小；达到稳态所需时间，状态x1，x2增长，状态x3，控制输入u变短； 当R变大时，反馈矩阵变小；各曲线波动幅值变小；达到稳态所需时间变长。 所以根据实际的系统允许，我们应该适当选择Q和R。精选ppt65%*MATLAB程序*%a=0 1

25、 0;0 0 1;0 0 -1/2;b=0;0;1/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);精选ppt66n=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

26、),R=2时状态响应曲线)grid,hold onfor j=1:n u(j,:)=-k*(x(j,:);end精选ppt67figure(2)plot(T,u);xlabel(time-s);ylabel(response);title(图(1.b) Q=diag(1,0,0),R=2时控制输入u的响应曲线)grid,hold on%*精选ppt68figure(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(

27、xa(:,3);精选ppt69Ta=0:60/na:60-60/na;plot(Ta,xa(:,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,:);end精选ppt70figure(4)plot(Ta,ua);xlabel(time-s);ylabel(response);title(图(2.b) Q=diag(1,0,0),R=2

28、000时控制输入u的响应曲线)grid,hold on%*精选ppt71figure(5)qb=10 0 0;0 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);精选ppt72Tb=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(1

29、0,0,0),R=2时状态响应曲线)grid,hold onfor j=1:nb ub(j,:)=-kb*(xb(j,:);end精选ppt73figure(6)plot(Tb,ub);xlabel(time-s);ylabel(response);title(图(3.b) Q=diag(10,0,0),R=2时控制输入u的响应曲线)grid,hold on%*精选ppt74figure(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)

30、;nc=length(xc(:,3);精选ppt75Tc=0:50/nc:50-50/nc;plot(Tc,xc(:,1),black,Tc,xc(:,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,:);end精选ppt76figure(8)plot(Tc,uc);xlabel(time-s);ylabel(response);title(图(4.b) Q=diag(

31、1,100,0),R=2时控制输入u的响应曲线)grid,hold on精选ppt7712.3 12.3 用用MATLABMATLAB解最优控制问题应用实例解最优控制问题应用实例12.3.1 12.3.1 导弹运动状态方程的建立导弹运动状态方程的建立12.3.2 12.3.2 最优导引律的求解与仿真验证最优导引律的求解与仿真验证精选ppt78在现有的自寻的导弹中，大都采用比例导引法。假设导弹和目标在同一平面内运动，按比例导引制导律，假设导弹的速度向量的旋转角速度 垂直于瞬时的弹目视线，并且正比于导弹与目标之间的视线角速率 ，假设目标的法向加速度为零，那么可得： (12-1)q qN精选ppt7

32、9其中， 为导弹的速度与基准方向的夹角， 为导弹与目标连线与基准方向的夹角，称为视线角， 是视线角速率， 是比例常数，称为导航比，通常为36。比例导引的实质是使导弹向着 减小的方向运动，抑制视线旋转，也就是使导弹的相对速度对准目标，保证导弹向着前置碰撞点飞行。qq Nq 精选ppt80 比例导引法是经典的导引方法。下面我们从最优控制理论的观点来研究自寻的导弹的最优导引规律问题。精选ppt8112.3.1 12.3.1 导弹运动状态方程的建立导弹运动状态方程的建立 导弹与目标的运动关系是非线性的，如果把导弹与目标的运动方程相对于理想弹道线性化，可得导弹运动的线性状态方程. 精选ppt82假设导弹

33、和目标在同一平面内运动，如图12-6所示。选 为固定坐标。导弹速度向量 与 轴成 角，目标速度向量 为与 轴成 角。导弹与目标的连线 与轴 成 角。假定导弹以尾追的方式攻击目标。坐标轴 和 的方向可以任意选择，使 和 都比较小。再假定导弹和目标均匀速飞行，也就是说 和 均为恒值。使用相对坐标状态变量，设 为导弹与目标在 轴方向上的距离偏差， 为导弹与目标在 轴方向上的距离偏差，即 (12-2)oxyMVoyTVoyTMToyqoxoyT,qMVTVxoxyoyMTMTyyyxxx精选ppt83xxyyMxTxMyTyMVTV0qTRTMH图图12-6 12-6 导弹和目标运动几何关系图导弹和目

34、标运动几何关系图精选ppt84假定 和 比较小，因此，则1cos, 1cos,sin,sinTTTT将上式对t求导，并根据导弹和目标的关系（如图12-6所示）可得 (12-3) coscossinsinMTTMTMTTMTVVyyyVVxxxMTMTMTTMTVVyyyVVxxx (12-4)精选ppt85以 表示 ， 表示 （即 ），则 (12-5) (12-6)1xx2xx 1x 21xx MTTVVxx2式中 表示目标的横向加速度， 表示导弹横向加速度，分别以 和 表示，那么 (12-7)TTVMVTaMaMTaax2精选ppt86导弹的横向加速度 为一控制量。一般将控制信号加给舵机，舵

35、面偏转后产生弹体攻角 ，而后产生横向加速度 。如果忽略舵机和弹体的惯性，而且假设控制量的单位与加速度单位相同，则可用控制量 来表示 ，也就是令 (12-8) 所以(12-7)式为： (12-9)MaMauMaMauuaxT2精选ppt87这样可得导弹运动状态方程为： (12-10) (12-11)21xx Taux2精选ppt88可写成矩阵的形式： (12-12)式中， ， ， ， 。 (12-13)如果不考虑目标的机动，即 ，则在这种情况下，式(12-12)变成： (12-14)TDaBuAXX21xxX0010A10B10D0TaBuAXX精选ppt89下面来考虑(12-4)式，该式可写成

36、 (12-15) 其中 表示导弹相对目标的接近速度。由于 的值都比较小， 可近似表示导弹与目标之间的距离。设 为导弹与目标的遭遇时刻(即导弹与目标相碰撞或两者之间的距离为最短的时刻)，则在某一瞬时 ，导弹与目标的距离 可近似用下式表示： (12-16)(TMVVyCTMVVVT,和qyftty)()()(ttVttVVtyfCfTM精选ppt90又考虑到对于导弹制导来说，最基本的要求是脱靶量越小越好，因此，应该选择最优控制量 ，使得下面的指标函数为最小。 (12-17)u22)()()()(fMfTfMfTtytytxtxJ精选ppt91然而，当要求一个反馈形式的控制时，按上式列出的问题很难求

37、解。所以我们以 时刻，即 时 的值 作为脱靶量，要求 值越小越好。另外，由于舵偏角受到限制，导弹结构能够承受的最大载荷也受到限制，所以控制信号 也应该受到限制。 ftt 0)()(ffCfttVty)(1ftx)(1ftxu精选ppt92因此，我们选择以下形式的二次型指标函数： (12-18)式中 ， 。 (12-19)即 (12-20)给定初始条件 ，应用最优控制理论，可以求出使 为最小的 。dtRuuQXXtCXtXJfttTTffT021212100ccC0000QdtRutCXtXJfttffT022121 0tXJu精选ppt93 由于系统是线性的，指标函数是二次型的，因此，求最优控

38、制规律就可以认为是一个求解线性二次型的过程。 对于线性二次型问题，可采用变分法、极小值原理、动态规划或其他方法求得最优控制精选ppt94 (12-21)式中 满足下列黎卡提矩阵微分方程 (12-22) 的终端条件为 (12-23)因此求解线性二次型问题的关键是求解黎卡提矩阵微分方程。 KXBRtuT1*QKBKBRKAKAKTT1KKCtKf精选ppt9512.3.2 12.3.2 最优导引律的求解与仿真验证最优导引律的求解与仿真验证当不考虑弹体惯性时，而且假定目标不机动，即，导弹运动状态方程为 (12-24)BuAXX精选ppt96指标函数为 (12-25)式中 ， , ， 。 dtRuuQ

39、XXtCXtXJfttTTffT021210010A10B2100ccC0000Q精选ppt97给出 时刻， 的初值 ，采用极小值原理可求得最优控制 为 (12-26)0tt 21xx 和 0201txtx和 tu* 2321 21 211212*34211 22231312fffffffcc ttcc ttc ttxcc ttxRRu tc ttc ttcc ttRRRR 精选ppt98在指标函数中，如不考虑导弹的相对运动速度 项，则可令 。 变成 (12-27)以 除上式的分子和分母，得 (12-28)2x02c tu* RttcRxttcxttctufff313122111*1c 312

40、21*333ttcRxttxtttufff精选ppt99为了使脱靶量为最小，应选取 ，则 (12-29)根据图12-6可得 1c ttxttxtuff221*3ttVxyxtgqfc11精选ppt100当 比较小时， ，则 (12-30) (12-31)qqtgq ttVxqfc1ttxttxVttVxttxqffcfcf2212111精选ppt101将上式代入(12-29)式，可得 (12-32)在上式中， 的单位是加速度的单位 。把 与导弹速度向量 的旋转角速度 联系起来，则有 (12-33) qVtuc3*u2秒米DVqVVVuMcM3u精选ppt102从(12-32)和(12-33)式

41、可以看出，当不考虑弹体惯性时，最优导引规律就是比例导引，其导航比为 。这证明了比例导引是一种很好的导引方法。最优导引规律的形成可用图12-7来表示。McVV3精选ppt103 下面将对最优导引律进行MATLAB仿真，并给出源代码和仿真结果。 cV3s1s1ttVfc121ttVfcq*u1x2x图图12-7 12-7 最优导引方框图最优导引方框图精选ppt104xxyyMxTxMyTyMVTV0HqTTMRTaMaMTLHE图图12-8 12-8 最优导引攻击几何平面最优导引攻击几何平面精选ppt105最优导引攻击几何关系如图12-8所示，在这里讨论的目标和导弹均认为是二维拦截几何平面上的质点

42、，分别以速度 和 运动。导弹的初始位置为相对坐标系的参考点，导弹初始速度矢量指向目标的初始位置， 为导弹的指令（垂直于视线）。TVMVMa精选ppt106其中： (12-34) (12-35) (12-36) 为目标速度在 轴上的分解， 是目标的角度。导弹和目标之间的接近速度为： (12-37)TTTVaTTTyVVcosTTTxVVsinTcTMVR ,TxTyVV, x y精选ppt107目标的速度分量可由其位置变化得到： (12-38)同样地，我们可以得到导弹的位置和速度的微分方程： ， (12-39) ， (12-40)TxTxTyTyVRVR,MxMxaVMyMyaVMxMxVRMy

43、MyVR精选ppt108上面几式中的下标x，y分别表示在x和y轴上的分量。 是导弹在地球坐标系的加速度分量。为了得到导弹的加速度分量，我们必须得到弹目的相对位移： (12-41) (12-42)MyMxaa,MxTxTMxRRRMyTyTMyRRR精选ppt109TMyTMxRRq1tanMxTxTMxVVVMyTyTMyVVV从图12-8中，根据三角关系我们可以得到视线角： (12-43)如果定义地球坐标系的速度分量为： (12-44) (12-45)精选ppt110我们可以根据视线角的公式求导后得到视线角速率： (12-46) (12-47)所以我们不难得出弹目的接近速度为： (12-48

44、)TMTMyTMxTMxTMyRVRVRq22122)(TMyTMxTMRRRTMTMyTMyTMxTMxTMcRVRVRRV)(精选ppt111根据最优导引制导律： (12-49)可得到导弹的加速的分量为： (12-50) (12-51) (12-52)qVVMc3)(tan1MyMxVVqVaMMxcosqVaMMysin精选ppt112 以上列出了两维的最优导引制导的必要方程，但是使用最优导引制导的导弹并不是直接向着目标发射的，而是向着一个能够导引导弹命中目标的方向发射，考虑了视线角之后可以得到导弹的指向角L。从图12-8中我们可以看出，如果导弹进入了碰撞三角区（如果目标和导弹同时保持匀

45、速直线运动，导弹必定会命中目标），这时利用正弦公式可以得到指向角的表达式：MTTVqVL)sin(sin1(12-53)精选ppt113 但是实际上导弹不可能能确切地在碰撞三角区发射，所以不能精确地得到拦截点。因为我们不知道目标将会如何机动，所以拦截点位置只能大概地估计。事实上，这也是需要导航系统的原因！初始时刻导弹偏离碰撞三角的角度称之为指向角误差（Head-Error）。考虑了导弹初始时刻的指向角和指向角误差之后，导弹的初始速度分量可以表示为：)cos()0(HeadErrorLqVVMMy)sin()0(HeadErrorLqVVMMx(12-54)(12-55)精选ppt114使用MA

46、TLAB编程，具体代码如下：%*MATLAB程序*%最优制导律仿真，初始化系统的参数clear all; %清除所有内存变量global SignVc;pi=3.14159265;Vm=1000;Vt=500;%导弹和目标的速度HeadError=0; %指向角误差精选ppt115ThetaT=pi; %目标的速度方向Rmx=0;Rmy=0; %导弹的位置Rtx=5000;Rty=10000;%目标的位置At=0; %目标法向加速度Vtx=Vt*sin(ThetaT);%目标的速度分量Vty=Vt*cos(ThetaT);Rtmx=Rtx-Rmx; %弹目相对距离Rtmy=Rty-Rmy;AmMax=15*9.8; %导弹的最大机动能力为15GRtm=sqrt(Rtmx2+Rtmy2);精选ppt116SightAngle=atan(Rtmx/Rtmy); %视线角LeadAngle=asin