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

第十二章 用MATLAB解最优控制问题及应用实例,,第十二章 用MATLAB解最优控制问题及应用实例,12.1 MATLAB工具简介 12.2 用MATLAB解线性二次型最优控制问题 12.3 用MATLAB解最优控制问题应用实例 12.4 小结,MATLAB是集数值运算、符号运算及图形处理等强大功能于一体的科学计算语言作为强大的科学计算平台，它几乎能满足所有的计算需求MATLAB具有编程方便、操作简单、可视化界面、优良的仿真图形环境、丰富的多学科工具箱等优点，尤其是在自动控制领域中MATLAB显示出更为强大的功能最优控制是在一定的约束条件下，从已给定的初始状态出发，确定最优控制作用的函数式，使目标函数为极小或极大在设计最优控制器的过程中，运用MATLAB最优控制设计工具，会大大减小设计的复杂性 在前面的几章中，我们已经介绍了一些最优控制方法，在本章中我们将介绍一个最优控制问题的应用实例，讨论如何使用最优控制方法来设计自寻的制导导弹的最优导引律，并采用MATLAB工具实现最优导引律，通过仿真来验证最优导引律的有效性12.1 MATLAB工具简介,,1, 系统模型的建立 系统的状态方程为：,在MATLAB中只需要将各个系数按照常规矩阵的方式输入到工作空间即可 ss(A,B,C,D),,,传递函数的零极点模型为：,在MATLAB中可以采用如下语句将零极点模型输入到工作空间:,,zpk(Z,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的控制系统工具箱中提供了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解线性二次型最优控制问题,一般情况的线性二次问题可表示如下： 设线性时变系统的方程为 其中， 为 维状态向量， 为 维控制向量， 为维输出向量。

寻找最优控制，使下面的性能指标最小,,,,,,,其中， 是 对称半正定常数阵， 是 对称半正定阵, 是 对称正定阵我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：,,,这个方程称为代数黎卡提方程代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解方法一：,,求解代数黎卡提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令 ，则可以写出下面的迭代公式,,,,,,%***************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; end end P=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是Algebraic Riccati Equation的缩写 对比前面给出的黎卡提方程，可以容易得出,,,,方法三： 我们也可以采用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,,线性系统为： ， 其目标函数是：,,确定最优控制解： 方法一： A=[0 1;-5,-3]; B=[0;1]; Q=[500 200;200 100]; R=1.6667; mylq K=inv(R)*B'*P P E,运行结果： K = 13.0276 6.7496 P = 67.9406 21.7131 21.7131 11.2495 E = -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.7496 P = 67.9406 21.7131 21.7131 11.2495 E =-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.2495 E = -7.2698 -2.4798 K =13.0276 6.7496 RR = 2.8458e-015,以上的三种方法的运行结果相同。

我们可以得到，最优控制变量与状态变量之间的关系： 在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：,,%***************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]),运行结果： 图12-1 最优控制曲线与最优状态曲线,该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。

因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量输出，y(4)是系统控制变量输出用plotyy函数进行双坐标显示，并设置相应的坐标范围以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我们只能使用第三种方法来求解，例如设置P的终值条件为[0.2;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.0961 E =-7.3052 -2.4723 K =13.0608 6.7775 RR =1.2847e-014 最优控制变量与状态变量之间的关系：,,例12-2,无人飞行器的最优高度控制，飞行器的控制方程如下,,,,是飞行器的高度； 是油门输入；设计控制律使得如下指标最小,,,,初始状态 。

绘制系统状态与控制输入，对如下给定的 矩阵进行仿真分析.,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-3所示：,,图12-3 状态。