chapter12-用MATLAB解最优控制问题及应用实例

第十二章用MATLAB解最优控制问题及应用实例第十二章用MATLAB解最优控制问题及应用实例12.1MATLAB工具简介12.2用MATLAB解线性二次型最优控制问题12.3用MATLAB解最优控制问题应用实例12.4小结MATLAB是集数值运算、符号运算及图形处理等强大功能于一体的科学计算语言。作为强大的科学计算平台，它几乎能满足所有的计算需求。MATLAB具有编程方便、操作简单、可视化界面、优良的仿真图形环境、丰富的多学科工具箱等优点，尤其是在自动控制领域中MATLAB显示出更为强大的功能。最优控制是在一定的约束条件下，从已给定的初始状态出发，确定最优控制作用的函数式，使目标函数为极小或极大。在设计最优控制器的过程中，运用MATLAB最优控制设计工具，会大大减小设计的复杂性。在前面的几章中，我们已经介绍了一些最优控制方法，在本章中我们将介绍一个最优控制问题的应用实例，讨论如何使用最优控制方法来设计自寻的制导导弹的最优导引律，并采用MATLAB工具实现最优导引律，通过仿真来验证最优导引律的有效性。12.1MATLAB工具简介DuCxyBuAxx1,系统模型的建立系统的状态方程为：在MATLAB中只需要将各个系数按照常规矩阵的方式输入到工作空间即可ss(A,B,C,D)],,,;;,,,;,,,[],,,;;,,,;,,,[],,,;;,,,;,,,[],,,;;,,,;,,,[212222111211212222111211212222111211212222111211qpqqppqnqqnnnpnnppnnnnnndddddddddDcccccccccCbbbbbbbbbBaaaaaaaaaA传递函数的零极点模型为：)())(()())(()(2121nmpspspszszszsKsG在MATLAB中可以采用如下语句将零极点模型输入到工作空间:];;;;[];;;;[;2121nmpppPzzzZKKGainzpk(Z,P,KGain)传递函数模型在更一般的情况下，可以表示为复数变量s的有理函数形式：nnnnnmmmmasasasasbsbsbsbsG122111121)(在MATLAB中可以采用如下语句将以上的传递函数模型输入到工作空间:G=tf(num,den);];,,,,,1[];,,,,[121121nnmmaaaadenbbbbnum2,系统模型的转换把其他形式转换成状态方程模型G1=ss(G)把其他形式转换成零极点模型G1=zpk(G)把其他形式转换成一般传递函数模型G1=tf(G)3,系统稳定性判据求出系统所有的极点，并观察系统是否有实部大于0的极点。系统由传递函数(num,den)描述roots(den)系统由状态方程(A,B,C,D)描述eig(A)4,系统的可控性与可观测性分析在MATLAB的控制系统工具箱中提供了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=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解线性二次型最优控制问题一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中，为维状态向量，为维控制向量，为维输出向量。)()()()()(tUtBtXtAtX)()()(tXtCtY()Xtn()Utm)(tYl寻找最优控制，使下面的性能指标最小fttTTffTdttUtRtUtetQtetPeteuJ0)()()()()()(21)()(21)(l其中，是对称半正定常数阵，是对称半正定阵,是对称正定阵。llll)(tQ)(tRmmP我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：1()()()()()()()()()()()TTKtKtAtAtKtKtBtRtBtKtQt可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成，这样可以保证系统状态渐进的趋近于零值，这样可以得出矩阵趋近于常值矩阵，且，这样上述黎卡提矩阵微分方程可以简化成为：ft)(tK0)(tK1()()()()()()()()()()0TTKtAtAtKtKtBtRtBtKtQt这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。方法一：求解代数黎卡提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令，则可以写出下面的迭代公式0011TTTTiiiiiEEEGWGGHEGWQAIAIE1BAIG12BAIQAIBRHTT11)(BAIQW1如果收敛于一个常数矩阵，即，则可以得出代数黎卡提方程的解为：上面的迭代算法可以用MATLAB来实现：1iii11112AIAIPiT%***************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;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(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是AlgebraicRiccatiEquation的缩写。对比前面给出的黎卡提方程，可以容易得出VTM,,0VXTXXMMXTAMTBBRT1QV方法三：我们也可以采用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’)计算）。采用care函数的优点在于可以设置P的终值条件，例如我们可以在下面的程序中设置P的终值条件为[0.2;0.2]。[P,E,K,RR]=care(A,B,Q,R,[0.2;0.2],eye(size(A)))采用lqr()函数不能设置代数黎卡提方程的边界条件。例12-1线性系统为：，其目标函数是：uxx1035100]6667.1[10020020050021dtuuxxJTT确定最优控制。解：方法一：A=[01;-5,-3];B=[0;1];Q=[500200;200100];R=1.6667;mylqK=inv(R)*B'*PPE运行结果：K=13.02766.7496P=67.940621.713121.713111.2495E=-0.11110.2222-1.1111-0.7778方法二:A=[01;-5,-3];B=[0;1];Q=[500200;200100];R=1.6667;[K,P,E]=lqr(A,B,Q,R)运行结果：K=13.02766.7496P=67.940621.713121.713111.2495E=-7.2698-2.4798方法三：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)))运行结果：P=67.940621.713121.713111.2495E=-7.2698-2.4798K=13.02766.7496RR=2.8458e-015以上的三种方法的运行结果相同。我们可以得到，最优控制变量与状态变量之间的关系：在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：)(6.7496)(13.0276)(21*txtxtu%***************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),[02.500.1]),axis(ax(2),[02.5-10])运行结果：图12-1最优控制曲线与最优状态曲线该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量输出，y(4)是系统控制变量输出。用plotyy函数进行双坐标显示，并设置相应的坐标范围。以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠。第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法。但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我