《第十二章用MATLAB解最优控制问题及应用实例演示文稿》由会员分享，可在线阅读，更多相关《第十二章用MATLAB解最优控制问题及应用实例演示文稿（127页珍藏版）》请在金锄头文库上搜索。
1、第十二章用MATLAB解最优控制问题及应用实例演示文稿第一页，共一百二十七页。ppt第十二章用MATLAB解最优控制问题及应用实例第二页，共一百二十七页。第十二章第十二章 用用MATLABMATLAB解最优控解最优控制问题及应用实例制问题及应用实例12.1 MATLAB12.1 MATLAB工具简介工具简介12.2 12.2 用用MATLABMATLAB解线性二次型最优控制问题解线性二次型最优控制问题12.3 12.3 用用MATLABMATLAB解最优控制问题应用实例解最优控制问题应用实例12.4 12.4 小结小结第三页，共一百二十七页。 MATLAB是集数值运算、符号运算及图形处理等强大
2、功能于一体的科学计算语言。作为强大的科学计算平台，它几乎能满足所有的计算需求。MATLAB具有编程方便、操作简单、可视化界面、优良的仿真图形环境、丰富的多学科工具箱等优点，尤其是在自动控制领域中MATLAB显示出更为强大的功能。第四页，共一百二十七页。 最优控制是在一定的约束条件下，从已给定的初始状态出发，确定最优控制作用的函数式，使目标函数为极小或极大。在设计最优控制器的过程中，运用MATLAB最优控制设计工具，会大大减小设计的复杂性。 在前面的几章中，我们已经介绍了一些最优控制方法，在本章中我们将介绍一个最优控制问题的应用实例，讨论如何使用最优控制方法来设计自寻的制导导弹的最优导引律，并采
3、用MATLAB工具实现最优导引律，通过仿真来验证最优导引律的有效性。第五页，共一百二十七页。12.1 MATLAB12.1 MATLAB工具简介工具简介1, 系统模型的建立 系统的状态方程为：第六页，共一百二十七页。 在MATLAB中只需要将各个系数按照常规矩阵的方式输入到工作空间即可 ss(A,B,C,D)第七页，共一百二十七页。传递函数的零极点模型为：在MATLAB中可以采用如下语句将零极点模型输入到工作空间:zpk(Z,P,KGain)第八页，共一百二十七页。传递函数模型在更一般的情况下，可以表示为复数变量s的有理函数形式：第九页，共一百二十七页。在MATLAB中可以采用如下语句将以上的
4、传递函数模型输入到工作空间: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()函数。该函数可以求
5、出系统的可控阶梯变换，该函数的调用格式为： 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()来直接求取系统的脉冲响应，该函数的可以有如下格式来调用：
6、y=impulse (G,t)第十四页，共一百二十七页。6, 系统的复域与频域分析 对于根轨迹的绘制，控制系统工具箱中给出了一个函数rlocus()函数来绘制系统的根轨迹，该函数的可以由如下格式来调用： R=rlocus(G,k)第十五页，共一百二十七页。 对于Nyquist曲线的绘制，控制系统工具箱中给出了一个函数nyquist()函数，该环数可以用来直接求解Nyquist阵列，绘制出Nyquist曲线，该函数的可以由如下格式来调用： rx,ry=nyquist(G,w) 对于Bode图，MATLAB控制工具箱中提供了bode()函数来求取、绘制系统的Bode图，该函数可以由下面的格式来调用
7、 mag,pha=bode(G,w)第十六页，共一百二十七页。12.2 12.2 用用MATLABMATLAB解线性二次型最优控制问题解线性二次型最优控制问题一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中， 为 维状态向量， 为 维控制向量， 为维输出向量。 第十七页，共一百二十七页。寻找最优控制，使下面的性能指标最小其中， 是 对称半正定常数阵， 是 对称半正定阵, 是 对称正定阵。第十八页，共一百二十七页。 我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：第十九页，共一百二十七页。可以看出，上述的黎卡提矩阵微分方程求解起来非常困
8、难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成 ，这样可以保证系统状态渐进的趋近于零值，这样可以得出矩阵趋近于常值矩阵 ，且 ，这样上述黎卡提矩阵微分方程可以简化成为：第二十页，共一百二十七页。这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。第二十一页，共一百二十七页。方法一：方法一：求解代数黎卡提方程的算法有很多，下面我们介绍一种简单的迭代算法来解该方程，令 ，则可以写出下面的迭代公式第二十二页，共一百二十七页。第二十三页，共一百二十七页。如果 收敛于一个常数矩阵，即 ，则可以得出代数黎卡提方程的解
9、为：上面的迭代算法可以用MATLAB来实现：第二十四页，共一百二十七页。%*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;第二十五页，共一百二十七页。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，方便我们以后调用来求
10、解代数黎卡提方程。第二十六页，共一百二十七页。方法二：方法二： 在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)第二十八页，共一百二十七页。其中， 矩阵满
11、足下列代数黎卡提方程，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(
12、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-112-1 线性系统为： ，其目标函数是：确定最优控制。第三十二页，共一百二十七页。解：解：方法一：A=0 1;-5,-3;B=0;1;Q=500 200;200 100;R=1.6667;mylqK=inv(R)*B
13、*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.9406 21.7131 21.7131 11.2495E =-7.2698 -2.4798第三十六页，共一百二十七页。方法三：A
14、=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第三十八页，共一百二十七页。 以上的三种方法的运行结果相同。我们可以得到，最优控制变量与状态变量之间的关系： 在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：第三十
15、九页，共一百二十七页。%*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)
16、,0 2.5 -1 0)第四十一页，共一百二十七页。运行结果：运行结果： 图图12-1 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)
