chapter12_用MATLAB解最优控制问题及应用实例
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
ss(A,B,C,D)
传递函数的零极点模型为：
( s z1 )(s z 2 ) ( s z m ) G( s) K ( s p1 )(s p2 ) ( s pn )
在MATLAB中可以采用如下语句将零极点模型 输入到工作空间:
KGain K ; Z [ z1 ; z 2 ;; z m ]; P [ p1 ; p2 ;; pn ];
这里的求解是建立在MATLAB的控制系统工具 箱中给出的一个基于Schur变换的黎卡提方 程求解函数are()基础上的，该函数的调用 格式为： X=are(M,T,V)
其中， M , T ,V 矩阵满足下列代数黎卡提方程，are 是Algebraic Riccati Equation的缩写。
MX XM T XTX V 0
一个函数step()来直接求取系统的阶跃响应，该函数
的可以有如下格式来调用： y=step(G,t) 对于系统的脉冲响应，控制系统工具箱中给出了 一个函数impulse()来直接求取系统的脉冲响应，该 函数的可以有如下格式来调用： y=impulse (G,t)
6, 系统的复域与频域分析
对于根轨迹的绘制，控制系统工具箱中给出了一
采用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
线性系统为：
1 0 0 x x u 5 3 1
Y (t ) C (t ) X (t )
Y (t ) U (t ) 为 m 维控制向量， X (t ) 为 n 维状态向量， 其中， l 为维输出向量。
l
寻找最优控制，使下面的性能指标最小
1 T 1 tf T J (u ) e (t f ) Pe (t f ) e (t )Q(t )e(t ) U T (t ) R(t )U (t ) dt 2 2 t0
可以看出，上述的黎卡提矩阵微分方程求解起来 非常困难，所以我们往往求出其稳态解。例如目 标函数中指定终止时间可以设置成 t f ， 这样可以保证系统状态渐进的趋近于零值，这样 (t ) 0 可以得出矩阵趋近于常值矩阵 K (t ) ，且 K ， 这样上述黎卡提矩阵微分方程可以简化成为：
i 1 E i E E i G W G i G H
T T T 1 T E i G W Q
E I A
1
1
I
1
A
G 2I A B
H R B (I A ) QI A B
T T 1
W QI A B
第十二章 用MATLAB解最 优控制问题及应用实例
第十二章 用MATLAB解最优 控制问题及应用实例
12.1 12.2 MATLAB工具简介 用MATLAB解线性二次型最优控制问题
12.3
12.4
用MATLAB解最优控制问题应用实例
小结
MATLAB是集数值运算、符号运算及图形处理 等强大功能于一体的科学计算语言。作为强大的 科学计算平台，它几乎能满足所有的计算需求。 MATLAB具有编程方便、操作简单、可视化界面、 优良的仿真图形环境、丰富的多学科工具箱等优 点，尤其是在自动控制领域中MATLAB显示出更为 强大的功能。
12.1
MATLAB工具简介
1, 系统模型的建立 系统的状态方程为：
Ax Bu x y Cx Du
在MATLAB中只需要将各个系数按照常规矩阵的 方式输入到工作空间即可
A [a11 , a12 ,, a1n ; a21 , a22 ,, a2n ;; an1 , an 2 ,, ann ] B [b11 , b12 ,, b1 p ; b21 , b22 ,, b2 p ;; bn1 , bn 2 ,, bnp ] C [c11 , c12 ,, c1n ; c21 , c22 ,, c2n ;; cq1 , cq 2 ,, cqn ] D [d11 , d12 ,, d1 p ; d 21 , d 22 ,, d 2 p ;; d q1 , d q 2 ,, d qp ]
运行结果： K = 13.0276 6.7496
P = 67.9406
21.7131 E =-7.2698 -2.4798
21.7131
11.2495
方法三：
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)))
的调用格式为： [Ac,Bc,Cc,Dc,Tc,Kc]=ctrbf(A,B,C) 在MATLAB的控制系统工具箱中提供了obsvf()函 数。该函数可以求出系统的可观测阶梯变换，该函 数的调用格式为： [Ao,Bo,Co,Do,To,Ko]=obsvf(A,B,C)
5, 系统的时域分析
对于系统的阶跃响应，控制系统工具箱中给出了
最优控制是在一定的约束条件下，从已给定的 初始状态出发，确定最优控制作用的函数式，使目 标函数为极小或极大。在设计最优控制器的过程中， 运用MATLAB最优控制设计工具，会大大减小设计的 复杂性。 在前面的几章中，我们已经介绍了一些最优控 制方法，在本章中我们将介绍一个最优控制问题的 应用实例，讨论如何使用最优控制方法来设计自寻 的制导导弹的最优导引律，并采用MATLAB工具实现 最优导引律，通过仿真来验证最优导引律的有效性。
1
如果 i 1 收敛于一个常数矩阵，即 i1 i ， 则可以得出代数黎卡提方程的解为：
P2IA
T 1
i 1 I A
源自文库1
上面的迭代算法可以用MATLAB来实现：
%***************MATLAB程序***************% I=eye(size(A));
对比前面给出的黎卡提方程，可以容易得出
M A
T BR B
V Q
1
T
方法三：
我们也可以采用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’)计算）。
num [b1 , b2 ,, bm , bm1 ]; den [1, a1 , a2 ,, an1 , an ];
G=tf(num,den);
2, 系统模型的转换 把其他形式转换成状态方程模型
G1=ss(G)
把其他形式转换成零极点模型 G1=zpk(G) 把其他形式转换成一般传递函数模型 G1=tf(G)
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，方便我们以后调用来
对于Bode图，MATLAB控制工具箱中提供了 bode()函数来求取、绘制系统的Bode图，该函数 可以由下面的格式来调用
[mag,pha]=bode(G,w)
12.2
用MATLAB解线性二次型最优控制问题
一般情况的线性二次问题可表示如下： 设线性时变系统的方程为 (t ) A(t ) X (t ) B(t )U (t ) X
Q(t ) 是 l l P 是 l l 对称半正定常数阵， 其中， 对称半正定阵, R(t ) 是 m m 对称正定阵。
我们用最小值原理求解上述问题，可以把上 述问题归结为求解如下黎卡提（Riccati）矩阵微 分方程：
(t ) K (t ) A(t ) AT (t )K (t ) K (t )B(t )R1 (t )BT (t )K (t ) Q(t ) K
求解代数黎卡提方程。
方法二： 在MATLAB的控制系统工具箱中提供了求解代数 黎卡提方程的函数lqr()，其调用的格式为： ［K,P,E］=lqr(A,B,Q,R) 式中输入矩阵为A,B,Q,R，其中(A,B)为给定的 对象状态方程模型，(Q,R)分别为加权矩阵Q和R； 返回矩阵K为状态反馈矩阵，P为代数黎卡提方程的 解，E为闭环系统的零极点。
3, 系统稳定性判据
求出系统所有的极点，并观察系统是否有实部大 于0的极点。 系统由传递函数 (num,den) 描述 roots(den) 系统由状态方程 (A,B,C,D) 描述 eig(A)
4, 系统的可控性与可观测性分析
在MATLAB的控制系统工具箱中提供了ctrbf()函
数。该函数可以求出系统的可控阶梯变换，该函数
，
其目标函数是：
1 T 500 200 T J x x u [1.6667]u dt 2 0 200 100
确定最优控制。
解：
方法一：
A=[0 1;-5,-3];
B=[0;1]; Q=[500 200;200 100]; R=1.6667; mylq K=inv(R)*B'*P P
运行结果： P = 67.9406 21.7131
21.7131
E = -7.2698 -2.4798 K =13.0276
11.2495
6.7496
RR = 2.8458e-015
zpk(Z,P,KGain)
传递函数模型在更一般的情况下，可以表示为复 数变量s的有理函数形式：
b1 s m b2 s m1 bm s bm1 G( s) n n 1 n2 s a1 s a2 s an1 s an
在MATLAB中可以采用如下语句将以上的传 递函数模型输入到工作空间:
K (t ) A(t ) AT (t ) K (t ) K (t ) B(t ) R1(t ) BT (t ) K (t) Q(t ) 0
这个方程称为代数黎卡提方程。代数黎卡提方程 的求解非常简单，并且其求解只涉及到矩阵运算， 所以非常适合使用MATLAB来求解。
方法一：
求解代数黎卡提方程的算法有很多，下面我们介 绍一种简单的迭代算法来解该方程，令 0 0 ， 则可以写出下面的迭代公式
E
运行结果：
K = 13.0276 6.7496
P = 67.9406
21.7131 E = -0.1111 -1.1111
21.7131
11.2495 0.2222 -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)
个函数rlocus()函数来绘制系统的根轨迹，该函数的
可以由如下格式来调用： R=rlocus(G,k)
对于Nyquist曲线的绘制，控制系统工具箱中 给出了一个函数nyquist()函数，该环数可以用来 直接求解Nyquist阵列，绘制出Nyquist曲线，该 函数的可以由如下格式来调用： [rx,ry]=nyquist(G,w)
传递函数的零极点模型为：
( s z1 )(s z 2 ) ( s z m ) G( s) K ( s p1 )(s p2 ) ( s pn )
在MATLAB中可以采用如下语句将零极点模型 输入到工作空间:
KGain K ; Z [ z1 ; z 2 ;; z m ]; P [ p1 ; p2 ;; pn ];
这里的求解是建立在MATLAB的控制系统工具 箱中给出的一个基于Schur变换的黎卡提方 程求解函数are()基础上的，该函数的调用 格式为： X=are(M,T,V)
其中， M , T ,V 矩阵满足下列代数黎卡提方程，are 是Algebraic Riccati Equation的缩写。
MX XM T XTX V 0
一个函数step()来直接求取系统的阶跃响应，该函数
的可以有如下格式来调用： y=step(G,t) 对于系统的脉冲响应，控制系统工具箱中给出了 一个函数impulse()来直接求取系统的脉冲响应，该 函数的可以有如下格式来调用： y=impulse (G,t)
6, 系统的复域与频域分析
对于根轨迹的绘制，控制系统工具箱中给出了一
采用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
线性系统为：
1 0 0 x x u 5 3 1
Y (t ) C (t ) X (t )
Y (t ) U (t ) 为 m 维控制向量， X (t ) 为 n 维状态向量， 其中， l 为维输出向量。
l
寻找最优控制，使下面的性能指标最小
1 T 1 tf T J (u ) e (t f ) Pe (t f ) e (t )Q(t )e(t ) U T (t ) R(t )U (t ) dt 2 2 t0
可以看出，上述的黎卡提矩阵微分方程求解起来 非常困难，所以我们往往求出其稳态解。例如目 标函数中指定终止时间可以设置成 t f ， 这样可以保证系统状态渐进的趋近于零值，这样 (t ) 0 可以得出矩阵趋近于常值矩阵 K (t ) ，且 K ， 这样上述黎卡提矩阵微分方程可以简化成为：
i 1 E i E E i G W G i G H
T T T 1 T E i G W Q
E I A
1
1
I
1
A
G 2I A B
H R B (I A ) QI A B
T T 1
W QI A B
第十二章 用MATLAB解最 优控制问题及应用实例
第十二章 用MATLAB解最优 控制问题及应用实例
12.1 12.2 MATLAB工具简介 用MATLAB解线性二次型最优控制问题
12.3
12.4
用MATLAB解最优控制问题应用实例
小结
MATLAB是集数值运算、符号运算及图形处理 等强大功能于一体的科学计算语言。作为强大的 科学计算平台，它几乎能满足所有的计算需求。 MATLAB具有编程方便、操作简单、可视化界面、 优良的仿真图形环境、丰富的多学科工具箱等优 点，尤其是在自动控制领域中MATLAB显示出更为 强大的功能。
12.1
MATLAB工具简介
1, 系统模型的建立 系统的状态方程为：
Ax Bu x y Cx Du
在MATLAB中只需要将各个系数按照常规矩阵的 方式输入到工作空间即可
A [a11 , a12 ,, a1n ; a21 , a22 ,, a2n ;; an1 , an 2 ,, ann ] B [b11 , b12 ,, b1 p ; b21 , b22 ,, b2 p ;; bn1 , bn 2 ,, bnp ] C [c11 , c12 ,, c1n ; c21 , c22 ,, c2n ;; cq1 , cq 2 ,, cqn ] D [d11 , d12 ,, d1 p ; d 21 , d 22 ,, d 2 p ;; d q1 , d q 2 ,, d qp ]
运行结果： K = 13.0276 6.7496
P = 67.9406
21.7131 E =-7.2698 -2.4798
21.7131
11.2495
方法三：
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)))
的调用格式为： [Ac,Bc,Cc,Dc,Tc,Kc]=ctrbf(A,B,C) 在MATLAB的控制系统工具箱中提供了obsvf()函 数。该函数可以求出系统的可观测阶梯变换，该函 数的调用格式为： [Ao,Bo,Co,Do,To,Ko]=obsvf(A,B,C)
5, 系统的时域分析
对于系统的阶跃响应，控制系统工具箱中给出了
最优控制是在一定的约束条件下，从已给定的 初始状态出发，确定最优控制作用的函数式，使目 标函数为极小或极大。在设计最优控制器的过程中， 运用MATLAB最优控制设计工具，会大大减小设计的 复杂性。 在前面的几章中，我们已经介绍了一些最优控 制方法，在本章中我们将介绍一个最优控制问题的 应用实例，讨论如何使用最优控制方法来设计自寻 的制导导弹的最优导引律，并采用MATLAB工具实现 最优导引律，通过仿真来验证最优导引律的有效性。
1
如果 i 1 收敛于一个常数矩阵，即 i1 i ， 则可以得出代数黎卡提方程的解为：
P2IA
T 1
i 1 I A
源自文库1
上面的迭代算法可以用MATLAB来实现：
%***************MATLAB程序***************% I=eye(size(A));
对比前面给出的黎卡提方程，可以容易得出
M A
T BR B
V Q
1
T
方法三：
我们也可以采用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’)计算）。
num [b1 , b2 ,, bm , bm1 ]; den [1, a1 , a2 ,, an1 , an ];
G=tf(num,den);
2, 系统模型的转换 把其他形式转换成状态方程模型
G1=ss(G)
把其他形式转换成零极点模型 G1=zpk(G) 把其他形式转换成一般传递函数模型 G1=tf(G)
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，方便我们以后调用来
对于Bode图，MATLAB控制工具箱中提供了 bode()函数来求取、绘制系统的Bode图，该函数 可以由下面的格式来调用
[mag,pha]=bode(G,w)
12.2
用MATLAB解线性二次型最优控制问题
一般情况的线性二次问题可表示如下： 设线性时变系统的方程为 (t ) A(t ) X (t ) B(t )U (t ) X
Q(t ) 是 l l P 是 l l 对称半正定常数阵， 其中， 对称半正定阵, R(t ) 是 m m 对称正定阵。
我们用最小值原理求解上述问题，可以把上 述问题归结为求解如下黎卡提（Riccati）矩阵微 分方程：
(t ) K (t ) A(t ) AT (t )K (t ) K (t )B(t )R1 (t )BT (t )K (t ) Q(t ) K
求解代数黎卡提方程。
方法二： 在MATLAB的控制系统工具箱中提供了求解代数 黎卡提方程的函数lqr()，其调用的格式为： ［K,P,E］=lqr(A,B,Q,R) 式中输入矩阵为A,B,Q,R，其中(A,B)为给定的 对象状态方程模型，(Q,R)分别为加权矩阵Q和R； 返回矩阵K为状态反馈矩阵，P为代数黎卡提方程的 解，E为闭环系统的零极点。
3, 系统稳定性判据
求出系统所有的极点，并观察系统是否有实部大 于0的极点。 系统由传递函数 (num,den) 描述 roots(den) 系统由状态方程 (A,B,C,D) 描述 eig(A)
4, 系统的可控性与可观测性分析
在MATLAB的控制系统工具箱中提供了ctrbf()函
数。该函数可以求出系统的可控阶梯变换，该函数
，
其目标函数是：
1 T 500 200 T J x x u [1.6667]u dt 2 0 200 100
确定最优控制。
解：
方法一：
A=[0 1;-5,-3];
B=[0;1]; Q=[500 200;200 100]; R=1.6667; mylq K=inv(R)*B'*P P
运行结果： P = 67.9406 21.7131
21.7131
E = -7.2698 -2.4798 K =13.0276
11.2495
6.7496
RR = 2.8458e-015
zpk(Z,P,KGain)
传递函数模型在更一般的情况下，可以表示为复 数变量s的有理函数形式：
b1 s m b2 s m1 bm s bm1 G( s) n n 1 n2 s a1 s a2 s an1 s an
在MATLAB中可以采用如下语句将以上的传 递函数模型输入到工作空间:
K (t ) A(t ) AT (t ) K (t ) K (t ) B(t ) R1(t ) BT (t ) K (t) Q(t ) 0
这个方程称为代数黎卡提方程。代数黎卡提方程 的求解非常简单，并且其求解只涉及到矩阵运算， 所以非常适合使用MATLAB来求解。
方法一：
求解代数黎卡提方程的算法有很多，下面我们介 绍一种简单的迭代算法来解该方程，令 0 0 ， 则可以写出下面的迭代公式
E
运行结果：
K = 13.0276 6.7496
P = 67.9406
21.7131 E = -0.1111 -1.1111
21.7131
11.2495 0.2222 -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)
个函数rlocus()函数来绘制系统的根轨迹，该函数的
可以由如下格式来调用： R=rlocus(G,k)
对于Nyquist曲线的绘制，控制系统工具箱中 给出了一个函数nyquist()函数，该环数可以用来 直接求解Nyquist阵列，绘制出Nyquist曲线，该 函数的可以由如下格式来调用： [rx,ry]=nyquist(G,w)