用MATLAB解线性二次型最优控制问题

1用MATLAB解线性二次型最优控制问题2解线性二次型最优控制问题一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中，为维状态向量，为维控制向量，为维输出向量。)()()()()(tUtBtXtAtX)()()(tXtCtY()Xtn()Utm)(tYl寻找最优控制，使下面的性能指标最小011()()()()()()()()()22ftTTTfftJUetPetetQtetUtRtUtdt其中，是对称半正定常数阵，是对称半正定阵,是对称正定阵。llll)(tQ)(tRmmP3解线性二次型最优控制问题我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：1()()()()()()()()()()()TTKtKtAtAtKtKtBtRtBtKtQt可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成，这样可以保证系统状态渐近的趋近于零值，这样可以得出矩阵趋近于常值矩阵，且，这样上述黎卡提矩阵微分方程可以简化成为：ft)(tK0)(tK10TTKAAKKBRBKQ这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。4解线性二次型最优控制问题方法一：求解代数黎卡提方程的算法有很多，下面介绍一种简单的迭代算法来解该方程。令，则可以写出下面的迭代公式0011TTTTiiiiiEEEGWGGHEGWQAIAIE1BAIG12BAIQAIBRHTT11)(BAIQW1如果收敛于一个常数矩阵，即，则可以得出代数黎卡提方程的解为：1iii11112AIAIPiT其中5解线性二次型最优控制问题%*************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;6解线性二次型最优控制问题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，方便我们以后调用来求解代数黎卡提方程。7解线性二次型最优控制问题方法二：在MATLAB的控制系统工具箱中提供了求解代数黎卡提方程的函数lqr()，其调用的格式为：［K,P,E］=lqr(A,B,Q,R)式中，输入矩阵为A,B,Q,R，其中(A,B)为给定的对象状态方程模型，(Q,R)分别为加权矩阵Q和R；返回矩阵K为状态反馈矩阵，P为代数黎卡提方程的解，E为闭环系统的极点。8解线性二次型最优控制问题这里的求解是建立在MATLAB的控制系统工具箱中给出的一个基于Schur分解的黎卡提方程求解函数are()基础上的，该函数的调用格式为：X=are(M,T,V)其中，矩阵满足下列代数黎卡提方程，are是AlgebraicRiccatiEquation的缩写。对比前面给出的黎卡提方程，可以容易得出VTM,,0VXTXXMMXTMATBBRT1VQ9解线性二次型最优控制问题方法三：我们也可以采用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’)计算10解线性二次型最优控制问题采用care函数的优点在于可以设置P的终值条件。例如，可以在下面的程序中设置P的终值条件为[0.2;0.2]。[P,E,K,RR]=care(A,B,Q,R,[0.2;0.2],eye(size(A)))采用lqr()函数不能设置代数黎卡提方程的边界条件。11解线性二次型最优控制问题例已知线性系统为目标函数uxx1035100]6667.1[10020020050021dtuuxxJTT确定最优控制。12解线性二次型最优控制问题解：方法一：A=[01;-5,-3];B=[0;1];Q=[500200;200100];R=1.6667;mylqK=inv(R)*B'*PPE=eig(A-B*K)运行结果：K=13.02766.7496P=67.940621.713121.713111.2495E=-2.4798-7.2698uxx1035100]6667.1[10020020050021dtuuxxJTT13解线性二次型最优控制问题方法二: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=-2.4798-7.269814解线性二次型最优控制问题方法三：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=3.4487e-01615解线性二次型最优控制问题以上的三种方法的运行结果相同。于是可以得到，最优控制变量与状态变量之间的关系：在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：*12()13.0276()6.7496()utxtxt16解线性二次型最优控制问题%***************MATLAB程序***************%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);figure('pos',[50,50,200,150],'color','w');axes('pos',[0.15,0.14,0.72,0.72])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])17解线性二次型最优控制问题运行结果：例图1最优控制曲线与最优状态曲线x1x2u*18解线性二次型最优控制问题该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量输出，y(4)是系统控制变量输出。用plotyy函数进行双坐标显示，并设置相应的坐标范围。以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠。第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法。但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我们只能使用第三种方法来求解，例如设置P的终值条件为[0.2;0.2]。19解线性二次型最优控制问题程序如下：%***************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)))20解线性二次型最优控制问题运行结果：P=67.723321.568521.568511.0961E=-7.3052-2.4723K=13.06086.7775RR=1.2847e-014最优控制变量与状态变量之间的关系：)(6.7775)(13.0608)(21*txtxtu21解线性二次型最优控制问题例无人飞行器的最优高度控制，飞行器的控制方程如下)(2/100)()()(2/100100010)()()(tuthththththth是飞行器的高度；是油门输入；设计控制律使得如下指标最小)(th)(tu20()100()1[(),()][(),(),()]000()()2()2000()()hthtJxtuthththththtutdththt22解线性二次型最优控制问题初始状态。绘制系统状态与控制输入，对如下给定的矩阵进行仿真分析.Tththth]0,0,10[)](),()([，,QR(a）100000,R2000Q2000R,000000001Q2,0000000010RQ10001000,2000QR(b）(c）(d）23解线性二次型最优控制问题解：线性二次型最优控制器设计如下：1）、Q=diag（1，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为k1=[0.70712.07722.0510]，u(t)=k1*x(t)；所画状态响应曲线及控制输入曲线如图所示。例图2状态响应曲线及控制输入曲线01020304050-50510time-sresponse图(1.a)Q=diag(1,0,0),R=2时状态响应曲线01020304050-8-6-4-202time-sresponse图(1.b)Q=diag(1,0,0),R=2时控制输入u的曲线24解线性二次型最优控制问题2）、Q=diag（1，0，0），R=2000时，由MATLAB求得最优状态反馈矩阵为k2=[0.02240.25170.4166]，u(t)=–k2*x(t)；所画状态响应曲线及控制输入曲线如图所示。例图3状态响应曲线及控制输入曲线01020304050-20246810time-sresponse图(2.a)Q=diag(1,0,0),R=2000时状态响应曲线01020304050-0.3-0.2-0.100.10.2time-sresponse图(2.b)Q=diag(1,0,0),R=2000时控制输入u的曲线25解线性二次型最优控制问题3）、Q=diag（10，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为k3=[2.23614.38923.3077]，u(t)=–k3*x(t)；所画状态响应曲线及控制输入曲线如图所示。例图4状态响应曲线及控制输入曲线01020304050-50510time-sresponse图(3.a)Q=diag(10,0,0),R=2时状态响应曲线01020304050-30-20-10010time-sresponse图(3