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

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

一般情况的线性二次问题可表示如下：设线性时变系统的方程为其中，为维状态向量，为维控制向量，为维输出向量。

寻找最优控制，使下面的性能指标最小其中，是对称半正定常数阵，是对称半正定阵,

是对称正定阵。第1页/共42页第一页，共43页。2解线性二次型最优控制问题我们用最小值原理求解上述问题，可以把上述问题归结为求解如下黎卡提（Riccati）矩阵微分方程：可以看出，上述的黎卡提矩阵微分方程求解起来非常困难，所以我们往往求出其稳态解。例如目标函数中指定终止时间可以设置成，这样可以保证系统状态渐近的趋近于零值，这样可以得出矩阵趋近于常值矩阵，且，这样上述黎卡提矩阵微分方程可以简化成为：这个方程称为代数黎卡提方程。代数黎卡提方程的求解非常简单，并且其求解只涉及到矩阵运算，所以非常适合使用MATLAB来求解。第2页/共42页第二页，共43页。3解线性二次型最优控制问题方法一：求解代数黎卡提方程的算法有很多，下面介绍一种简单的迭代算法来解该方程。令，则可以写出下面的迭代公式如果收敛于一个常数矩阵，即，则可以得出代数黎卡提方程的解为：其中第3页/共42页第三页，共43页。4解线性二次型最优控制问题%*************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;第4页/共42页第四页，共43页。5解线性二次型最优控制问题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，方便我们以后调用来求解代数黎卡提方程。第5页/共42页第五页，共43页。6解线性二次型最优控制问题方法二：在MATLAB的控制系统工具箱中提供了求解代数黎卡提方程的函数lqr()，其调用的格式为：［K,P,E］=lqr(A,B,Q,R)式中，输入矩阵为A,B,Q,R，其中(A,B)为给定的对象状态方程模型，(Q,R)分别为加权矩阵Q和R；返回矩阵K为状态反馈矩阵，P为代数黎卡提方程的解，E为闭环系统的极点。第6页/共42页第六页，共43页。7解线性二次型最优控制问题这里的求解是建立在MATLAB的控制系统工具箱中给出的一个基于Schur分解的黎卡提方程求解函数are()基础上的，该函数的调用格式为：

X=are(M,T,V)其中，矩阵满足下列代数黎卡提方程，are是AlgebraicRiccatiEquation的缩写。对比前面给出的黎卡提方程，可以容易得出第7页/共42页第七页，共43页。8解线性二次型最优控制问题方法三：我们也可以采用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’)计算第8页/共42页第八页，共43页。9解线性二次型最优控制问题采用care函数的优点在于可以设置P的终值条件。例如，可以在下面的程序中设置P的终值条件为[0.2;0.2]。

[P,E,K,RR]=care(A,B,Q,R,[0.2;0.2],eye(size(A)))采用lqr()函数不能设置代数黎卡提方程的边界条件。第9页/共42页第九页，共43页。10解线性二次型最优控制问题例已知线性系统为目标函数确定最优控制。第10页/共42页第十页，共43页。11解线性二次型最优控制问题解：方法一：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.2698第11页/共42页第十一页，共43页。12解线性二次型最优控制问题方法二: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.2698第12页/共42页第十二页，共43页。13解线性二次型最优控制问题方法三：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-016第13页/共42页第十三页，共43页。14解线性二次型最优控制问题以上的三种方法的运行结果相同。于是可以得到，最优控制变量与状态变量之间的关系：在以上程序的基础上，可以得到在最优控制的作用下的最优控制曲线与最优状态曲线，其程序如下：第14页/共42页第十四页，共43页。15解线性二次型最优控制问题%***************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])第15页/共42页第十五页，共43页。16解线性二次型最优控制问题运行结果：例图1最优控制曲线与最优状态曲线x1x2u*第16页/共42页第十六页，共43页。17解线性二次型最优控制问题该程序采用augstate函数将状态变量作为输出变量，用于显示；输出项作为最优控制的输出。因此，阶跃响应输出y中，y(1)是系统输出，y(2)和y(3)是状态变量输出，y(4)是系统控制变量输出。用plotyy函数进行双坐标显示，并设置相应的坐标范围。以上三种方法中，第一种方法易于理解黎卡提方程的解法，其解法简单但是并不可靠。第二种方法比起另两种方法使用方便，不易出错，所以我们推荐使用这种方法。但是采用lqr()函数不能设置代数黎卡提方程的边界条件，所以，如果题目设置了P的终值条件，我们只能使用第三种方法来求解，例如设置P的终值条件为[0.2;0.2]。第17页/共42页第十七页，共43页。18解线性二次型最优控制问题程序如下：%***************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)))第18页/共42页第十八页，共43页。19解线性二次型最优控制问题运行结果：P=67.723321.568521.568511.0961E=-7.3052-2.4723K=13.06086.7775RR=1.2847e-014最优控制变量与状态变量之间的关系：第19页/共42页第十九页，共43页。20解线性二次型最优控制问题例无人飞行器的最优高度控制，飞行器的控制方程如下是飞行器的高度；是油门输入；设计控制律使得如下指标最小第20页/共42页第二十页，共43页。21解线性二次型最优控制问题初始状态。绘制系统状态与控制输入，对如下给定的矩阵进行仿真分析.(a）(b）(c）(d）第21页/共42页第二十一页，共43页。22解线性二次型最优控制问题解：线性二次型最优控制器设计如下：1）、Q=diag（1，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为

k1=[0.70712.07722.0510]，u(t)=k1*x(t)；所画状态响应曲线及控制输入曲线如图所示。例图2状态响应曲线及控制输入曲线第22页/共42页第二十二页，共43页。23解线性二次型最优控制问题2）、Q=diag（1，0，0），R=2000时，由MATLAB求得最优状态反馈矩阵为k2=[0.02240.25170.4166]，u(t)=–k2*x(t)；所画状态响应曲线及控制输入曲线如图所示。例图3状态响应曲线及控制输入曲线第23页/共42页第二十三页，共43页。24解线性二次型最优控制问题3）、Q=diag（10，0，0），R=2时，由MATLAB求得最优状态反馈矩阵为

k3=[2.23614.38923.3077]，u(t)=–k3*x(t)；所画状态响应曲线及控制输入曲线如图所示。例图4状态响应曲线及控制输入曲线第24页/共42页第二十四页，共43页。25解线性二次型最优控制问题4）、Q=diag（1，100，0），R=2时，由MATLAB求得最优状态反馈矩阵为

k4=[0.70717.61124.6076]，u(t)=–k4*x(t)；所画状态响应曲线及控制输入曲线如图所示。例图5状态响应曲线及控制输入曲线第25页/共42页第二十五页，共43页。26解线性二次型最优控制问题由1），2），3），4）可分析如下：例图3与例图2相比，当Q不变，R增大时，各相应曲线达到稳态所需时间增长，即响应变慢；但波动幅值变小，反馈矩阵变小；例图4与例图2和例图3相比，当Q对角线上第1个元素增大时，各相应曲线达到稳态所需时间变短，即响应快；但波动幅值变大，反馈矩阵增大；由例图5可知，当Q对角线上第2个元素增大时，状态x1，x2曲线达到稳态所需时间较长，即响应较慢，平缓的趋于零；状态x3，控制输入u达到稳态所需时间短，即响应快；状态x2，x3波动幅值较小，比例图2和例图4小，比例图3稍大，控制输入u波动幅值比例图2和例图4小，比例图3大；反馈矩阵最大。第26页/共42页第二十六页，共43页。27解线性二次型最优控制问题综上所述可得结论：Q=diag（1，0，0），R=2时，系统各方面响应较好。矩阵Q变大时，反馈矩阵变大；当Q的对角线上第1个元素变大时，各曲线波动幅值变大，达到稳态所需时间变短；当Q的对角线上第2个元素变大时，各曲线波动幅值变小，达到稳态所需时间、状态x1、x2增长，状态x3，控制输入u变短。当R变大时，反馈矩阵变小；各曲线波动幅值变小；达到稳态所需时间变长。所以根据实际的系统允许，我们应该适当选择Q和R。程序第27页/共42页第二十七页，共43页。28解线性二次型最优控制问题矩阵Q变大时，反馈矩阵变大；当Q的对角线上第1个元素变大时，各曲线波动幅值变大，达到稳态所需时间变短。结论第28页/共42页第二十八页，共43页。29解线性二次型最优控制问题矩阵Q变大时，反馈矩阵变大；当Q的对角线上第2个元素变大时，各曲线波动幅值变小，达到稳态所需时间、状态x1、x2增长，状态x3，控制输入u变短。结论第29页/共42页第二十九页，共43页。30解线性二次型最优控制问题当R变大时，反馈矩阵变小；各曲线波动幅值变小；达到稳态所需时间变长。结论第30页/共42页第三十页，共43页。31解线性二次型最优控制问题%********MATLAB程序***myUAV.m******%a=[010;001;00-1/2];b=[0;0;1/2];c=[100;010;001];d=[0;0;0];figure(1)q=[100;000;000];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);第31页/共42页第三十一页，共43页。32解线性二次型最优控制问题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时状态响应曲线')gridforj=1:nu(j,:)=-k*(x(j,:))';endfigure(2)plot(T,u);xlabel('time-s');ylabel('response');title('图(1.b)Q=diag(1,0,0),R=2时控制输入u的曲线')grid第32页/共42页第三十二页，共43页。33解线性二次型最优控制问题%**************************figure(3)qa=[100;000;000];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));第33页/共42页第三十三页，共43页。34解线性二次型最优控制问题Ta=0:60/na:60-60/na;plot(Ta,xa(:,1),'black',Ta,xa(:,2),'red',Ta,xa(:,3),'green');xlabel('time-s');ylabel('response');title('图(2.a)Q=diag(1,0,0),R=2000时状态响应曲线')gridforj=1:naua(j,:)=-ka*(xa(j,:))';end第34页/共42页第三十四页，共43页。35解线性二次型最优控制问题figure(4)plot(Ta,ua);xlabel('time-s');ylabel('response');title('图(2.b)Q=diag(1,0,0),R=2000时控制输入u的曲线')grid%%%*******************************第35页/共42页第三十五页，共43页。36解线性二次型最优控制问题figure(5)qb=[1000;000;000];rb=2;[kb,pb,eb]=lqr(a,b,qb,rb)x0=[10;0;0];ab1=a-b*kb;[yb,xb]=initial(ab1,b,c,d,x0,20);nb=length(xb(:,3));第36页/共42页第三十六页，共43页。37解线性二次型最优控制问题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时状态响应曲线')gridforj=1:nbub(j,:)=-kb*(xb(j,:))';end第37页/共42页第三十七页，共43页。38解线性二次型最优控制问题figure(6)plot(Tb,ub);xlabel('time-s');ylabel('response');title('图(3.b)