用MATLAB解最优控制问题及应用实例课件

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

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

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

12.1MATLAB工具简介12.2用MATLAB解线性二次型最优控制问题12.3用MATLAB解最优控制问题应用实例12.4小结2第十二章用MATLAB解最优控制问题及应用实例

12.1MATLAB是集数值运算、符号运算及图形处理等强大功能于一体的科学计算语言。作为强大的科学计算平台，它几乎能满足所有的计算需求。MATLAB具有编程方便、操作简单、可视化界面、优良的仿真图形环境、丰富的多学科工具箱等优点，尤其是在自动控制领域中MATLAB显示出更为强大的功能。3MATLAB是集数值运算、符号运算及图形处理等强

最优控制是在一定的约束条件下，从已给定的初始状态出发，确定最优控制作用的函数式，使目标函数为极小或极大。在设计最优控制器的过程中，运用MATLAB最优控制设计工具，会大大减小设计的复杂性。在前面的几章中，我们已经介绍了一些最优控制方法，在本章中我们将介绍一个最优控制问题的应用实例，讨论如何使用最优控制方法来设计自寻的制导导弹的最优导引律，并采用MATLAB工具实现最优导引律，通过仿真来验证最优导引律的有效性。4最优控制是在一定的约束条件下，从已给定的初始状态12.1MATLAB工具简介

1,系统模型的建立系统的状态方程为：512.1MATLAB工具简介

1,系统模型的建立5

在MATLAB中只需要将各个系数按照常规矩阵的方式输入到工作空间即可

ss(A,B,C,D)6在MATLAB中只需要将各个系数按照常规矩阵的方传递函数的零极点模型为：在MATLAB中可以采用如下语句将零极点模型输入到工作空间:zpk(Z,P,KGain)7传递函数的零极点模型为：在MATLAB中可以采用如下语句将零传递函数模型在更一般的情况下，可以表示为复数变量s的有理函数形式：8传递函数模型在更一般的情况下，可以表示为复数变量s的有理函数在MATLAB中可以采用如下语句将以上的传递函数模型输入到工作空间:G=tf(num,den);9在MATLAB中可以采用如下语句将以上的传递函数模型输入到工2,系统模型的转换把其他形式转换成状态方程模型

G1=ss(G)

把其他形式转换成零极点模型

G1=zpk(G)

把其他形式转换成一般传递函数模型

G1=tf(G)102,系统模型的转换103,系统稳定性判据求出系统所有的极点，并观察系统是否有实部大于0的极点。系统由传递函数(num,den)描述

roots(den)

系统由状态方程(A,B,C,D)描述

eig(A)113,系统稳定性判据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)124,系统的可控性与可观测性分析125,系统的时域分析对于系统的阶跃响应，控制系统工具箱中给出了一个函数step()来直接求取系统的阶跃响应，该函数的可以有如下格式来调用：

y=step(G,t)

对于系统的脉冲响应，控制系统工具箱中给出了一个函数impulse()来直接求取系统的脉冲响应，该函数的可以有如下格式来调用：

y=impulse(G,t) 135,系统的时域分析136,系统的复域与频域分析对于根轨迹的绘制，控制系统工具箱中给出了一个函数rlocus()函数来绘制系统的根轨迹，该函数的可以由如下格式来调用：

R=rlocus(G,k)146,系统的复域与频域分析14

对于Nyquist曲线的绘制，控制系统工具箱中给出了一个函数nyquist()函数，该环数可以用来直接求解Nyquist阵列，绘制出Nyquist曲线，该函数的可以由如下格式来调用：

[rx,ry]=nyquist(G,w)

对于Bode图，MATLAB控制工具箱中提供了bode()函数来求取、绘制系统的Bode图，该函数可以由下面的格式来调用

[mag,pha]=bode(G,w)15对于Nyquist曲线的绘制，控制系统工具箱中给12.2用MATLAB解线性二次型最优控制问题

一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中，为维状态向量，为维控制向量，为维输出向量。1612.2用MATLAB解线性二次型最优控制问题

一般情况寻找最优控制，使下面的性能指标最小其中，是对称半正定常数阵，是对称半正定阵,是对称正定阵。17寻找最优控制，使下面的性能指标最小其中，是对称半正

我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：18我们用最小值原理求解上述问题，可以把上述问题归结可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成，这样可以保证系统状态渐进的趋近于零值，这样可以得出矩阵趋近于常值矩阵，且，这样上述黎卡提矩阵微分方程可以简化成为：19可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。20这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并方法一：求解代数黎卡提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令，则可以写出下面的迭代公式21方法一：求解代数黎卡提方程的算法有很多，下面我们介绍一种简单2222如果收敛于一个常数矩阵，即，则可以得出代数黎卡提方程的解为：上面的迭代算法可以用MATLAB来实现：23如果收敛于一个常数矩阵，即%***************MATLAB程序***************%I=eye(size(A));iA=inv(I-A);E=iA*(I+A);G=2*iA^2*B;H=R+B'*iA'*Q*iA*B;W=Q*iA*B;P0=zeros(size(A));i=0;24%***************MATLAB程序******while(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，方便我们以后调用来求解代数黎卡提方程。25while(1),i=i+1;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方法二：26这里的求解是建立在MATLAB的控制系统工具箱中给出的一个基于Schur变换的黎卡提方程求解函数are()基础上的，该函数的调用格式为：

X=are(M,T,V)27这里的求解是建立在MATLAB的控制系统工具箱中给出的一个基其中，矩阵满足下列代数黎卡提方程，are是AlgebraicRiccatiEquation的缩写。对比前面给出的黎卡提方程，可以容易得出28其中，矩阵满足下列代数黎卡提方程，are是Al方法三：我们也可以采用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(sum(diag(Res’*Res)))，或者用Norm(Res’,fro’)计算）。29方法三：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采用care函数的优点在于可以设置P的终值条件例12-1

线性系统为：，其目标函数是：确定最优控制。31例12-1 线性系统为：解：方法一：A=[01;-5,-3];B=[0;1];Q=[500200;200100];R=1.6667;mylqK=inv(R)*B'*PPE32解：32运行结果：K=13.02766.7496P=67.940621.713121.713111.2495E=-0.11110.2222-1.1111-0.777833运行结果：33方法二:A=[01;-5,-3];B=[0;1];Q=[500200;200100];R=1.6667;[K,P,E]=lqr(A,B,Q,R)34方法二:34运行结果：K=13.02766.7496P=67.940621.713121.713111.2495E=-7.2698-2.479835运行结果：35方法三：A=[01;-5,-3];B=[0;1];Q=[500200;200100];R=1.6667;[P,E,K,RR]=care(A,B,Q,R,zeros(size(B)),eye(size(A)))36方法三：36运行结果：P=67.940621.713121.713111.2495E=-7.2698-2.4798K=13.02766.7496RR=2.8458e-01537运行结果：37

以上的三种方法的运行结果相同。我们可以得到，最优控制变量与状态变量之间的关系：在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：38以上的三种方法的运行结果相同。我们可以得%***************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;39%***************MATLAB程序******[ap,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),[02.500.1]),axis(ax(2),[02.5-10])40[ap,bp,cp,dp]=augstate(ap,bp,C运行结果：

图12-1最优控制曲线与最优状态曲线41运行结果：41

该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量输出，y(4)是系统控制变量输出。用plotyy函数进行双坐标显示，并设置相应的坐标范围。42该程序采用augstate函数将状态变量作为输出

以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠。第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法。但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我们只能使用第三种方法来求解，例如设置P的终值条件为[0.2;0.2]。43以上三种方法中，第一种方法易于理解黎卡提方程的解程序如下：%***************MATLAB程序***************%A=[01;-5,-3];B=[0;1];Q=[500200;200100];R=1.6667;[P,E,K,RR]=care(A,B,Q,R,[0.2;0.2],eye(size(A)))44程序如下：44运行结果：P=67.723321.568521.568511.0961E=-7.3052-2.4723K=13.06086.7775RR=1.2847e-014最优控制变量与状态变量之间的关系：45运行结果：45例12-2

无人飞行器的最优高度控制，飞行器的控制方程如下46例12-2 无人飞行器的最优高度控制，飞行器的控制方程如下

是飞行器的高度；是油门输入；设计控制律使得如下指标最小初始状态。绘制系统状态与控制输入，对如下给定的矩阵进行仿真分析.47是飞行器的高度；是油门输入；a).b).c).d).4848解：线性二次型最优控制指标如下：其中Q和R分别是对状态变量和控制量的加权矩阵，线性二次型最优控制器设计如下：1）、Q=diag（1，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为

k1=[0.70712.07722.0510]，u（t）=–k1*x（t）；所画状态响应曲线及控制输入响应曲线如下图12-2所示：49解：线性二次型最优控制指标如下：495050图12-2状态响应曲线及控制输入响应曲线51图12-2状态响应曲线及控制输入响应曲线512）、Q=diag（1，0，0），R=2000时，由MATLAB求得最优状态反馈矩阵为k2=[0.02240.25170.4166]，u（t）=–k2*x（t）；

所画状态响应曲线及控制输入响应曲线如下图12-3所示：522）、Q=diag（1，0，0），R=2000时，由MATL5353图12-3状态响应曲线及控制输入响应曲线54图12-3状态响应曲线及控制输入响应曲线543）、Q=diag（10，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为

k3=[2.23614.38923.3077]，u（t）=–k3*x（t）；所画状态响应曲线及控制输入响应曲线如下图12-4所示：553）、Q=diag（10，0，0），R=2时，由MATLAB5656图12-4状态响应曲线及控制输入响应曲线57图12-4状态响应曲线及控制输入响应曲线574）、Q=diag（1，100，0），R=2时，由MATLAB求得最优状态反馈矩阵为

k4=[0.70717.61124.6076]，u（t）=–k4*x（t）；所画状态响应曲线及控制输入响应曲线如下图12-5所示：584）、Q=diag（1，100，0），R=2时，由MATLA5959图12-5状态响应曲线及控制输入响应曲线60图12-5状态响应曲线及控制输入响应曲线60由1），2），3），4）可分析如下：图12-3与图12-2相比，当Q不变，R增大时，各相应曲线达到稳态所需时间增长，即响应变慢；但波动幅值变小，反馈矩阵变小；图12-4与图12-2和图12-3相比，当Q对角线上第1个元素增大时，各相应曲线达到稳态所需时间变短，即响应快；但波动幅值变大，反馈矩阵增大；61由1），2），3），4）可分析如下：图12-3与图12-2相由图12-5可知，当Q对角线上第2个元素增大时，状态x1，x2曲线达到稳态所需时间较长，即响应较慢，平缓的趋于零；状态x3，控制输入u达到稳态所需时间短，即响应快；状态x2，x3波动幅值较小，比图12-2和图12-4小，比图12-3稍大，控制输入u波动幅值比图12-2和图12-4小，比图12-3大；反馈矩阵最大。62由图12-5可知，当Q对角线上第2个元素增大时，状态x1，x

综上所述可得结论：Q=diag（1，0，0），R=2时，系统各方面响应较好。矩阵Q变大时，反馈矩阵变大；当Q的对角线上第1个元素变大时，各曲线波动幅值变大，达到稳态所需时间变短；

63综上所述可得结论：Q=diag（1，0，0），R=2时

当Q的对角线上第2个元素变大时，各曲线波动幅值变小；达到稳态所需时间，状态x1，x2增长，状态x3，控制输入u变短；当R变大时，反馈矩阵变小；各曲线波动幅值变小；达到稳态所需时间变长。所以根据实际的系统允许，我们应该适当选择Q和R。64当Q的对角线上第2个元素变大时，各曲线波动幅值变小；%***************MATLAB程序***************%a=[010;001;00-1/2];b=[0;0;1/2];c=[100;010;001];d=[0;0;0];figure(1)q=[100;000;000];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);65%***************MATLAB程序******n=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,holdonforj=1:nu(j,:)=-k*(x(j,:))';end66n=length(x(:,3));66figure(2)plot(T,u);xlabel('time-s');ylabel('response');title('图(1.b)Q=diag(1,0,0),R=2时控制输入u的响应曲线')grid,holdon%**************************67figure(2)67figure(3)qa=[100;000;000];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));68figure(3)68Ta=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,holdonforj=1:naua(j,:)=-ka*(xa(j,:))';end69Ta=0:60/na:60-60/na;69figure(4)plot(Ta,ua);xlabel('time-s');ylabel('response');title('图(2.b)Q=diag(1,0,0),R=2000时控制输入u的响应曲线')grid,holdon%%%*******************************70figure(4)70figure(5)qb=[1000;000;000];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));71figure(5)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,holdonforj=1:nbub(j,:)=-kb*(xb(j,:))';end72Tb=0:20/nb:20-20/nb;72figure(6)plot(Tb,ub);xlabel('time-s');ylabel('response');title('图(3.b)Q=diag(10,0,0),R=2时控制输入u的响应曲线')grid,holdon%%%*************73figure(6)73figure(7)qc=[100;01000;000];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));74figure(7)74Tc=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,holdonforj=1:ncuc(j,:)=-kc*(xc(j,:))';end75Tc=0:50/nc:50-50/nc;75figure(8)plot(Tc,uc);xlabel('time-s');ylabel('response');title('图(4.b)Q=diag(1,100,0),R=2时控制输入u的响应曲线')grid,holdon76figure(8)7612.3用MATLAB解最优控制问题应用实例

12.3.1导弹运动状态方程的建立12.3.2最优导引律的求解与仿真验证7712.3用MATLAB解最优控制问题应用实例

12.3.在现有的自寻的导弹中，大都采用比例导引法。假设导弹和目标在同一平面内运动，按比例导引制导律，假设导弹的速度向量的旋转角速度垂直于瞬时的弹目视线，并且正比于导弹与目标之间的视线角速率，假设目标的法向加速度为零，那么可得：

(12-1)78在现有的自寻的导弹中，大都采用比例导引法。假设导弹和目标在同其中，为导弹的速度与基准方向的夹角，为导弹与目标连线与基准方向的夹角，称为视线角，是视线角速率，是比例常数，称为导航比，通常为3~6。比例导引的实质是使导弹向着减小的方向运动，抑制视线旋转，也就是使导弹的相对速度对准目标，保证导弹向着前置碰撞点飞行。79其中，为导弹的速度与基准方向的夹角，为导弹与目标

比例导引法是经典的导引方法。下面我们从最优控制理论的观点来研究自寻的导弹的最优导引规律问题。80比例导引法是经典的导引方法。下面我们从12.3.1导弹运动状态方程的建立

导弹与目标的运动关系是非线性的，如果把导弹与目标的运动方程相对于理想弹道线性化，可得导弹运动的线性状态方程.8112.3.1导弹运动状态方程的建立

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

(12-2)82假设导弹和目标在同一平面内运动，如图12-6所示。选为图12-6导弹和目标运动几何关系图83图12-6导弹和目标运动几何关系图83假定和比较小，因此，则将上式对t求导，并根据导弹和目标的关系（如图12-6所示）可得

(12-3)

(12-4)84假定和比较小，因此，将上式对t求导，并根据导以表示，表示（即），则

(12-5)

(12-6)式中表示目标的横向加速度，表示导弹横向加速度，分别以和表示，那么

(12-7)85以表示，表示（即），则式中导弹的横向加速度为一控制量。一般将控制信号加给舵机，舵面偏转后产生弹体攻角，而后产生横向加速度。如果忽略舵机和弹体的惯性，而且假设控制量的单位与加速度单位相同，则可用控制量来表示，也就是令

(12-8)

所以(12-7)式为：

(12-9)86导弹的横向加速度为一控制量。一般将控制信号加给舵机，舵这样可得导弹运动状态方程为：

(12-10)

(12-11)87这样可得导弹运动状态方程为：87可写成矩阵的形式：

(12-12)式中，

，，，。(12-13)如果不考虑目标的机动，即，则在这种情况下，式(12-12)变成：

(12-14)88可写成矩阵的形式：88下面来考虑(12-4)式，该式可写成

(12-15)

其中表示导弹相对目标的接近速度。由于的值都比较小，可近似表示导弹与目标之间的距离。设为导弹与目标的遭遇时刻(即导弹与目标相碰撞或两者之间的距离为最短的时刻)，则在某一瞬时，导弹与目标的距离可近似用下式表示：

(12-16)89下面来考虑(12-4)式，该式可写成89又考虑到对于导弹制导来说，最基本的要求是脱靶量越小越好，因此，应该选择最优控制量，使得下面的指标函数为最小。

(12-17)90又考虑到对于导弹制导来说，最基本的要求是脱靶量越小越好，因此然而，当要求一个反馈形式的控制时，按上式列出的问题很难求解。所以我们以时刻，即时的值作为脱靶量，要求值越小越好。另外，由于舵偏角受到限制，导弹结构能够承受的最大载荷也受到限制，所以控制信号也应该受到限制。91然而，当要求一个反馈形式的控制时，按上式列出的问题很难求解。因此，我们选择以下形式的二次型指标函数：

(12-18)式中，。(12-19)即

(12-20)给定初始条件，应用最优控制理论，可以求出使为最小的。92因此，我们选择以下形式的二次型指标函数：92

由于系统是线性的，指标函数是二次型的，因此，求最优控制规律就可以认为是一个求解线性二次型的过程。对于线性二次型问题，可采用变分法、极小值原理、动态规划或其他方法求得最优控制93由于系统是线性的，指标函数是二次型的，(12-21)式中满足下列黎卡提矩阵微分方程

(12-22)

的终端条件为

(12-23)因此求解线性二次型问题的关键是求解黎卡提矩阵微分方程。94

12.3.2最优导引律的求解与仿真验证

当不考虑弹体惯性时，而且假定目标不机动，即，导弹运动状态方程为

(12-24)9512.3.2最优导引律的求解与仿真验证

当不考虑弹体惯性时指标函数为

(12-25)式中

，,，。96指标函数为96给出时刻，的初值，采用极小值原理可求得最优控制为

(12-26)97给出时刻，的初值在指标函数中，如不考虑导弹的相对运动速度项，则可令。变成

(12-27)以除上式的分子和分母，得

(12-28)98在指标函数中，如不考虑导弹的相对运动速度项，则可令为了使脱靶量为最小，应选取，则

(12-29)根据图12-6可得

99为了使脱靶量为最小，应选取，则99当比较小时，，则

(12-30)

(12-31)100当比较小时，，则100将上式代入(12-29)式，可得

(12-32)在上式中，的单位是加速度的单位。把与导弹速度向量的旋转角速度联系起来，则有

(12-33)101将上式代入(12-29)式，可得101从(12-32)和(12-33)式可以看出，当不考虑弹体惯性时，最优导引规律就是比例导引，其导航比为。这证明了比例导引是一种很好的导引方法。最优导引规律的形成可用图12-7来表示。102从(12-32)和(12-33)式可以看出，当不考虑弹体惯性

下面将对最优导引律进行MATLAB仿真，并给出源代码和仿真结果。图12-7最优导引方框图103图12-7最优导引方框图103图12-8最优导引攻击几何平面104图12-8最优导引攻击几何平面104最优导引攻击几何关系如图12-8所示，在这里讨论的目标和导弹均认为是二维拦截几何平面上的质点，分别以速度和运动。导弹的初始位置为相对坐标系的参考点，导弹初始速度矢量指向目标的初始位置，为导弹的指令（垂直于视线）。105最优导引攻击几何关系如图12-8所示，在这里讨论的目标和导弹其中： (12-34)

(12-35)

(12-36)

为目标速度在轴上的分解，是目标的角度。导弹和目标之间的接近速度为：

(12-37)106其中： (12-目标的速度分量可由其位置变化得到：

(12-38)同样地，我们可以得到导弹的位置和速度的微分方程：， (12-39)

， (12-40)107目标的速度分量可由其位置变化得到：107上面几式中的下标x，y分别表示在x和y轴上的分量。是导弹在地球坐标系的加速度分量。为了得到导弹的加速度分量，我们必须得到弹目的相对位移：

(12-41) (12-42)108上面几式中的下标x，y分别表示在x和y轴上的分量。从图12-8中，根据三角关系我们可以得到视线角：

(12-43)如果定义地球坐标系的速度分量为：

(12-44)

(12-45)109从图12-8中，根据三角关系我们可以得到视线角：109我们可以根据视线角的公式求导后得到视线角速率：

(12-46) (12-47)所以我们不难得出弹目的接近速度为：

(12-48)110我们可以根据视线角的公式求导后得到视线角速率：110根据最优导引制导律：

(12-49)可得到导弹的加速的分量为：

(12-50) (12-51) (12-52)111根据最优导引制导律：111

以上列出了两维的最优导引制导的必要方程，但是使用最优导引制导的导弹并不是直接向着目标发射的，而是向着一个能够导引导弹命中目标的方向发射，考虑了视线角之后可以得到导弹的指向角L。从图12-8中我们可以看出，如果导弹进入了碰撞三角区（如果目标和导弹同时保持匀速直线运动，导弹必定会命中目标），这时利用正弦公式可以得到指向角的表达式：

(12-53)112以上列出了两维的最优导引制导的必要方程，但是使用

但是实际上导弹不可能能确切地在碰撞三角区发射，所以不能精确地得到拦截点。因为我们不知道目标将会如何机动，所以拦截点位置只能大概地估计。事实上，这也是需要导航系统的原因！初始时刻导弹偏离碰撞三角的角度称之为指向角误差（Head-Error）。考虑了导弹初始时刻的指向角和指向角误差之后，导弹的初始速度分量可以表示为：

(12-54)

(12-55)113但是实际上导弹不可能能确切地在碰撞三角区发射，所使用MATLAB编程，具体代码如下：%***************MATLAB程序***************%%最优制导律仿真，初始化系统的参数clearall; %清除所有内存变量globalSignVc;pi=3.14159265;Vm=1000;Vt=500;%导弹和目标的速度HeadError=0;%指向角误差114使用MATLAB编程，具体代码如下：114ThetaT=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(Rtmx^2+Rtmy^2);115ThetaT=pi;%目标的速度方向115SightAngle=atan(Rtmx/Rtmy);%视线角LeadAngle