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

1、第十二章 用MATLAB解最优控制问题及应用实例,第十二章 用MATLAB解最优控制问题及应用实例,12.1 MATLAB工具简介12.2 用MATLAB解线性二次型最优控制问题12.3 用MATLAB解最优控制问题应用实例12.4 小结,MATLAB是集数值运算、符号运算及图形处理等强大功能于一体的科学计算语言。作为强大的科学计算平台，它几乎能满足所有的计算需求。MATLAB具有编程方便、操作简单、可视化界面、优良的仿真图形环境、丰富的多学科工具箱等优点，尤其是在自动控制领域中MATLAB显示出更为强大的功能。,最优控制是在一定的约束条件下，从已给定的初始状态出发，确定最优控制作用的函数式，

2、使目标函数为极小或极大。在设计最优控制器的过程中，运用MATLAB最优控制设计工具，会大大减小设计的复杂性。在前面的几章中，我们已经介绍了一些最优控制方法，在本章中我们将介绍一个最优控制问题的应用实例，讨论如何使用最优控制方法来设计自寻的制导导弹的最优导引律，并采用MATLAB工具实现最优导引律，通过仿真来验证最优导引律的有效性。,12.1 MATLAB工具简介,1,系统模型的建立 系统的状态方程为：,在MATLAB中只需要将各个系数按照常规矩阵的方式输入到工作空间即可 ss(A,B,C,D),传递函数的零极点模型为：,在MATLAB中可以采用如下语句将零极点模型输入到工作空间:,zpk(Z,

3、P,KGain),传递函数模型在更一般的情况下，可以表示为复数变量s的有理函数形式：,在MATLAB中可以采用如下语句将以上的传递函数模型输入到工作空间:G=tf(num,den);,2,系统模型的转换 把其他形式转换成状态方程模型 G1=ss(G)把其他形式转换成零极点模型 G1=zpk(G)把其他形式转换成一般传递函数模型 G1=tf(G),3,系统稳定性判据 求出系统所有的极点，并观察系统是否有实部大于0的极点。系统由传递函数(num,den)描述 roots(den)系统由状态方程(A,B,C,D)描述 eig(A),4,系统的可控性与可观测性分析 在MATLAB的控制系统工具箱中提供

4、了ctrbf()函数。该函数可以求出系统的可控阶梯变换，该函数的调用格式为：Ac,Bc,Cc,Dc,Tc,Kc=ctrbf(A,B,C)在MATLAB的控制系统工具箱中提供了obsvf()函数。该函数可以求出系统的可观测阶梯变换，该函数的调用格式为：Ao,Bo,Co,Do,To,Ko=obsvf(A,B,C),5,系统的时域分析 对于系统的阶跃响应，控制系统工具箱中给出了一个函数step()来直接求取系统的阶跃响应，该函数的可以有如下格式来调用：y=step(G,t)对于系统的脉冲响应，控制系统工具箱中给出了一个函数impulse()来直接求取系统的脉冲响应，该函数的可以有如下格式来调用：y=

5、impulse(G,t),6,系统的复域与频域分析 对于根轨迹的绘制，控制系统工具箱中给出了一个函数rlocus()函数来绘制系统的根轨迹，该函数的可以由如下格式来调用：R=rlocus(G,k),对于Nyquist曲线的绘制，控制系统工具箱中给出了一个函数nyquist()函数，该环数可以用来直接求解Nyquist阵列，绘制出Nyquist曲线，该函数的可以由如下格式来调用：rx,ry=nyquist(G,w)对于Bode图，MATLAB控制工具箱中提供了bode()函数来求取、绘制系统的Bode图，该函数可以由下面的格式来调用 mag,pha=bode(G,w),12.2 用MATLAB解

6、线性二次型最优控制问题,一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中，为 维状态向量，为 维控制向量，为维输出向量。,寻找最优控制，使下面的性能指标最小,其中，是 对称半正定常数阵，是 对称半正定阵,是 对称正定阵。,我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：,这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。,方法一：,求解代数黎卡提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令，则可以写出下面的迭代公式,%*MATLAB程序*%I=e

7、ye(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;,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，方便我们以后调用来求解代数黎卡提方程。,方法二：在MATLAB的控制系统工具箱中提供了求解代数黎卡提方程的函数lqr()，其调用的格式为：K,P,E=lqr

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

9、=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)计算）。,采用care函数的优点在于可以设置P的终值条件，例如我们可以在下面的程序中设置P的终值条件为0.2;0.2。P,E,K,RR=care(A,B,Q,R,0.2;0.2,eye(size(A)采用lqr

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

11、9406 21.7131 21.7131 11.2495E=-7.2698-2.4798,方法三：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),运行结果：P=67.9406 21.7131 21.7131 11.2495E=-7.2698-2.4798K=13.0276 6.7496RR=2.8458e-015,以上的三种方法的运行结果相同。我们可以得到，最优控制变量与状态变量之间的关系：在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状

12、态曲线，其程序如下：,%*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;,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),0 2.5 0 0.1),axis(ax(2),0 2.5-1 0),

13、运行结果：图12-1 最优控制曲线与最优状态曲线,该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量输出，y(4)是系统控制变量输出。用plotyy函数进行双坐标显示，并设置相应的坐标范围。,以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠。第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法。但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我们只能使用第三种方法来求解，例如设置P的终值条件为0.2;

14、0.2。,程序如下：%*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;0.2,eye(size(A),运行结果：P=67.7233 21.5685 21.5685 11.0961E=-7.3052-2.4723K=13.0608 6.7775RR=1.2847e-014最优控制变量与状态变量之间的关系：,例12-2,无人飞行器的最优高度控制，飞行器的控制方程如下,是飞行器的高度；是油门输入；设计控制律使得如下指标最小,初始状态。绘制系统状态与控制输入，对如下给定的 矩阵进行仿真分

15、析.,a).b).c).d).,解：线性二次型最优控制指标如下：其中Q和R分别是对状态变量和控制量的加权矩阵，线性二次型最优控制器设计如下：1）、Q=diag（1，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为 k1=0.7071 2.0772 2.0510，u（t）=k1*x（t）；所画状态响应曲线及控制输入响应曲线如下图12-2所示：,图12-2 状态响应曲线及控制输入响应曲线,2）、Q=diag（1，0，0），R=2000时，由MATLAB求得最优状态反馈矩阵为k2=0.0224 0.2517 0.4166，u（t）=k2*x（t）；所画状态响应曲线及控制输入响应曲线如下图12

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

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

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

19、;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);,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,hold onfor j=1:n u(j,:)=-k*(x(j,:);end,figure(2)plot(T,u);xlabel(ti

20、me-s);ylabel(response);title(图(1.b)Q=diag(1,0,0),R=2时控制输入u的响应曲线)grid,hold on%*,figure(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);,Ta=0:60/na:60-60/na;plot(Ta,xa(:,1),black,Ta,xa(:,2),red,Ta,xa(:,3),green);xlabel(time-s)

21、;ylabel(response);title(图(2.a)Q=diag(1,0,0),R=2000时状态响应曲线)grid,hold onfor j=1:na ua(j,:)=-ka*(xa(j,:);end,figure(4)plot(Ta,ua);xlabel(time-s);ylabel(response);title(图(2.b)Q=diag(1,0,0),R=2000时控制输入u的响应曲线)grid,hold on%*,figure(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

22、,xb=initial(ab1,b,c,d,x0,20);nb=length(xb(:,3);,Tb=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,:);end,figure(6)plot(Tb,ub);xlabel(time-s);ylabel(response);title(

23、图(3.b)Q=diag(10,0,0),R=2时控制输入u的响应曲线)grid,hold on%*,figure(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);,Tc=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

24、)Q=diag(1,100,0),R=2时状态响应曲线)grid,hold onfor j=1:nc uc(j,:)=-kc*(xc(j,:);end,figure(8)plot(Tc,uc);xlabel(time-s);ylabel(response);title(图(4.b)Q=diag(1,100,0),R=2时控制输入u的响应曲线)grid,hold on,12.3 用MATLAB解最优控制问题应用实例,12.3.1 导弹运动状态方程的建立12.3.2 最优导引律的求解与仿真验证,在现有的自寻的导弹中，大都采用比例导引法。假设导弹和目标在同一平面内运动，按比例导引制导律，假设导弹的速

25、度向量的旋转角速度 垂直于瞬时的弹目视线，并且正比于导弹与目标之间的视线角速率，假设目标的法向加速度为零，那么可得：(12-1),比例导引法是经典的导引方法。下面我们从最优控制理论的观点来研究自寻的导弹的最优导引规律问题。,12.3.1 导弹运动状态方程的建立,导弹与目标的运动关系是非线性的，如果把导弹与目标的运动方程相对于理想弹道线性化，可得导弹运动的线性状态方程.,图12-6 导弹和目标运动几何关系图,假定 和 比较小，因此，则,(12-4),以 表示，表示（即），则(12-5)(12-6),式中 表示目标的横向加速度，表示导弹横向加速度，分别以 和 表示，那么(12-7),这样可得导弹运

26、动状态方程为：(12-10)(12-11),由于系统是线性的，指标函数是二次型的，因此，求最优控制规律就可以认为是一个求解线性二次型的过程。对于线性二次型问题，可采用变分法、极小值原理、动态规划或其他方法求得最优控制,12.3.2 最优导引律的求解与仿真验证,当不考虑弹体惯性时，而且假定目标不机动，即，导弹运动状态方程为(12-24),下面将对最优导引律进行MATLAB仿真，并给出源代码和仿真结果。,图12-7 最优导引方框图,图12-8 最优导引攻击几何平面,以上列出了两维的最优导引制导的必要方程，但是使用最优导引制导的导弹并不是直接向着目标发射的，而是向着一个能够导引导弹命中目标的方向发射

27、，考虑了视线角之后可以得到导弹的指向角L。从图12-8中我们可以看出，如果导弹进入了碰撞三角区（如果目标和导弹同时保持匀速直线运动，导弹必定会命中目标），这时利用正弦公式可以得到指向角的表达式：,(12-53),但是实际上导弹不可能能确切地在碰撞三角区发射，所以不能精确地得到拦截点。因为我们不知道目标将会如何机动，所以拦截点位置只能大概地估计。事实上，这也是需要导航系统的原因！初始时刻导弹偏离碰撞三角的角度称之为指向角误差（Head-Error）。考虑了导弹初始时刻的指向角和指向角误差之后，导弹的初始速度分量可以表示为：,(12-54)(12-55),使用MATLAB编程，具体代码如下：%*M

28、ATLAB程序*%最优制导律仿真，初始化系统的参数clear all;%清除所有内存变量global SignVc;pi=3.14159265;Vm=1000;Vt=500;%导弹和目标的速度HeadError=0;%指向角误差,ThetaT=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=s

29、qrt(Rtmx2+Rtmy2);,SightAngle=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;%弹目的相对运动速度,Vc=-(Rtmx*Vtmx+Rtmy*Vtmy)/Rtm;SignVc=sign(Vc);%Vc的符号Time=0;TimeStep=0.1

30、;%时间和时间步长file=fopen(output.txt,w);%将数据写入文件%循环while(1)%Vc改变符号仿真结束if(sign(Vc)=SignVc)break;else,if(Rtm100)TimeStep=0.005;endSignVc=sign(Vc);%Vc的符号%视线角速率dSightAngle=(Rtmy*Vtmx-Rtmx*Vtmy)/(Rtm2);dTheta=3*Vc*dSightAngle/Vm;Theta=atan(Vmx/Vmy);%导弹加速度,Am=Vm*dTheta;%限制机动能力if(AmAmMax)Am=AmMax;end%加速度分量Amx=Am

31、*cos(SightAngle);Amy=-Am*sin(SightAngle);Time=Time+TimeStep;%目标位置,Rtx=Rtx+TimeStep*Vtx;Rty=Rty+TimeStep*Vty;dThetaT=At/Vt;ThetaT=ThetaT+TimeStep*dThetaT;Vtx=Vt*sin(ThetaT);%目标的速度分量Vty=Vt*cos(ThetaT);%导弹位置Rmx=Rmx+TimeStep*Vmx;Rmy=Rmy+TimeStep*Vmy;%导弹速度,Vmx=Vmx+TimeStep*Amx;Vmy=Vmy+TimeStep*Amy;Vm=sqr

32、t(Vmx2+Vmy2);%弹目相对位移Rtmx=Rtx-Rmx;Rtmy=Rty-Rmy;%上一步的脱靶量Rtm0=Rtm;Rtm=sqrt(Rtmx2+Rtmy2);SightAngle=atan(Rtmx/Rtmy);%视线角%弹目相对速度,Vtmx=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);,仿真结果：,图12-9 指向

33、误差，目标不机动的攻击情况,图12-10,指向误差，目标不机动时导弹的加速度,图12-11,指向误差，目标不机动时攻击情况,图12-12,指向误差，目标不机动时导弹加速度,12.4 小结,本章首先简单介绍了一下MATLAB控制系统分析工具的使用，接着介绍了如何用MATLAB工具设计线性二次型最优控制器。线性二次型最优控制可以归结为求解黎卡提（Riccati）方程。最后，在此基础上，我们给出了两个运用MATLAB求解线性二次型最优控制问题的例子。,此外，本章还讨论了如何使用最优控制方法来设计自寻的制导导弹的最优导引律。通过最优二次型控制问题的求解，我们发现：当不考虑弹体惯性时，应用最优控制理论所得到的最优导引规律与目前广泛采用的比例导引法相一致。因此，从现代控制理论的观点来看，比例导引法是一种很好的导引方法。最后，我们还给出了最优导引律的仿真验证程序和仿真曲线。,