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

1 用MATLAB解线性二次型最优控制问题 2 解线性二次型最优控制问题 一般情况的线性二次问题可表示如下 设线性时变系统的方程为其中 为维状态向量 为维控制向量 为维输出向量 寻找最优控制 使下面的性能指标最小 其中 是对称半正定常数阵 是对称半正定阵 是对称正定阵 3 解线性二次型最优控制问题 我们用最小值原理求解上述问题 可以把上述问题归结为求解如下黎卡提 Riccati 矩阵微分方程 可以看出 上述的黎卡提矩阵微分方程求解起来非常困难 所以我们往往求出其稳态解 例如目标函数中指定终止时间可以设置成 这样可以保证系统状态渐近的趋近于零值 这样可以得出矩阵趋近于常值矩阵 且 这样上述黎卡提矩阵微分方程可以简化成为 这个方程称为代数黎卡提方程 代数黎卡提方程的求解非常简单 并且其求解只涉及到矩阵运算 所以非常适合使用MATLAB来求解 4 解线性二次型最优控制问题 方法一 求解代数黎卡提方程的算法有很多 下面介绍一种简单的迭代算法来解该方程 令 则可以写出下面的迭代公式 如果收敛于一个常数矩阵 即 则可以得出代数黎卡提方程的解为 其中 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的缩写 对比前面给出的黎卡提方程 可以容易得出 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 解线性二次型最优控制问题 例已知线性系统为 目标函数 确定最优控制 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 2698 13 解线性二次型最优控制问题 方法二 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 14 解线性二次型最优控制问题 方法三 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 15 解线性二次型最优控制问题 以上的三种方法的运行结果相同 于是可以得到 最优控制变量与状态变量之间的关系 在以上程序的基础上 可以得到在最优控制的作用下的最优控制曲线与最优状态曲线 其程序如下 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最优控制曲线与最优状态曲线 x1 x2 u 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最优控制变量与状态变量之间的关系 21 解线性二次型最优控制问题 例无人飞行器的最优高度控制 飞行器的控制方程如下 是飞行器的高度 是油门输入 设计控制律使得如下指标最小 22 解线性二次型最优控制问题 初始状态 绘制系统状态与控制输入 对如下给定的矩阵进行仿真分析 a b c d 23 解线性二次型最优控制问题 解 线性二次型最优控制器设计如下 1 Q diag 1 0 0 R 2时 由MATLAB求得最优状态反馈矩阵为k1 0 70712 07722 0510 u t k1 x t 所画状态响应曲线及控制输入曲线如图所示 例图2状态响应曲线及控制输入曲线 24 解线性二次型最优控制问题 2 Q diag 1 0 0 R 2000时 由MATLAB求得最优状态反馈矩阵为k2 0 02240 25170 4166 u t k2 x t 所画状态响应曲线及控制输入曲线如图所示 例图3状态响应曲线及控制输入曲线 25 解线性二次型最优控制问题 3 Q diag 10 0 0 R 2时 由MATLAB求得最优状态反馈矩阵为k3 2 23614 38923 3077 u t k3 x t 所画状态响应曲线及控制输入曲线如图所示 例图4状态响应曲线及控制输入曲线 26 解线性二次型最优控制问题 4 Q diag 1 100 0 R 2时 由MATLAB求得最优状态反馈矩阵为k4 0 70717 61124 6076 u t k4 x t 所画状态响应曲线及控制输入曲线如图所示 例图5状态响应曲线及控制输入曲线 27 解线性二次型最优控制问题 由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大 反馈矩阵最大 28 解线性二次型最优控制问题 综上所述可得结论 Q diag 1 0 0 R 2时 系统各方面响应较好 矩阵Q变大时 反馈矩阵变大 当Q的对角线上第1个元素变大时 各曲线波动幅值变大 达到稳态所需时间变短 当Q的对角线上第2个元素变大时 各曲线波动幅值变小 达到稳态所需时间 状态x1 x2增长 状态x3 控制输入u变短 当R变大时 反馈矩阵变小 各曲线波动幅值变小 达到稳态所需时间变长 所以根据实际的系统允许 我们应该适当选择Q和R 程序 29 解线性二次型最优控制问题 矩阵Q变大时 反馈矩阵变大 当Q的对角线上第1个元素变大时 各曲线波动幅值变大 达到稳态所需时间变短 结论 30 解线性二次型最优控制问题 矩阵Q变大时 反馈矩阵变大 当Q的对角线上第2个元素变大时 各曲线波动幅值变小 达到稳态所需时间 状态x1 x2增长 状态x3 控制输入u变短 结论 31 解线性二次型最优控制问题 当R变大时 反馈矩阵变小 各曲线波动幅值变小 达到稳态所需时间变长 结论 32 解线性二次型最优控制问题 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 33 解线性二次型最优控制问题 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 34 解线性二次型最优控制问题 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 35 解线性二次型最优控制问题 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 36 解线性二次型最优控制问题 figure 4 plot Ta ua xlabel time s ylabel response title 图 2 b Q diag 1 0 0 R 2000时控制输入u的曲线 grid 37 解线性二次型最优控制问题 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 38 解线性二次型最优控制问题 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