声振论坛 首页 基础理论 控制理论 查看内容

Lyapunov、Sylvester和Riccati方程的Matlab求解

2017-1-24 20:40| 发布者: Vickyvictoria| 查看: 489| 评论: 4|原作者: iApollo|来自: 新浪iApollo的博客

摘要: Lyapunov、Sylvester和Riccati方程是控制系统常用到的几个方程,应用和计算比较广泛。
Lyapunov、Sylvester和Riccati方程是控制系统常用到的几个方程,应用和计算比较广泛。

一、Lyapunov方程

1、连续Lyapunov方程

连续Lyapunov方程可以表示为:
01.gif
Lyapunov方程来源与微分方程稳定性理论,其中要求C为对称正定的n×n方阵,从而可以证明解X亦为n×n对称矩阵,这类方程直接求解比较困难,不过有了Matlab中lyap()函数,就简单多了。

>> A=[1 2 3;4 5 6;7 8 0]
>> C=-[10 5 4;5 6 7;4 7 9]
>> X=lyap(A,C)

2、Lyapunov方程的解析解

利用Kroncecker乘积的表示方法,可以将Lyapunov方程写为:
02.gif
可见,方程有唯一解的条件并不局限与C对称正定,只要满足非奇异即可保证方程唯一解。同时也打破了传统观念,C必须对称正定的。

  1. function x=lyap2(A,C)
  2. %Lyapunov方程的符号解法
  3. n=size(C,1);
  4. A0=kron(A,eye(n))+kron(eye(n),A);
  5. c=C(:);
  6. x0=-inv(A0)*c;
  7. x=reshape(x0,n,n)
复制代码

恩下面看一个示例,体会下符号解法:

>>A=[1 2 3;4 5 6;7 8 0];
>>C=-[10 5 4;5 6 7;4 7 9];
>>x=lyap2(sym(A),sym(C))

x =
[ -71/18,   35/9,   7/18]
[   35/9,  -25/9,    2/9]
[   7/18,    2/9,   -1/9]

3、离散Lyapunov方程

离散Lyapunov方程的一般形式为:
03.gif

Matlab中直接提供了dlyap()函数求解该方程:X=dlyap(A,Q)

其实,如果A矩阵非奇异,则等式两边同时右乘06.gif得到:
03.gif
就可以将其变换成连续的Sylvester方程
05.gif
而Sylvester方程是广义Lyapunov方程,故离散的Lyapunov方程还可以使用下面的方法求解:
  1. B=-inv(A’)
  2. C=Q*inv(A’)
  3. X=lyap(A,B,C)
复制代码

下面总结下我们上面的讲到的知识点:
  • X=lyap(A,C)                                 连续Lyapunov方程数值解法
  • X=lyap2(A,C)                               连续Lyapunov方程符号解法
  • X=lyap(A,B,C)                                广义Lyapunov方程,即Sylvester方程
  • X=dlyap(A,Q)或者X=lyap(A,-inv(A’),Q*inv(A’))    离散Lyapunov方程
二、Sylvester方程Matlab求解
Sylvester方程的一般形式为:
07.gif
该方程又称为广义的Lyapunov方程,式中A是n×n方阵,B是m×m方阵,X和C是n×m矩阵。Matlab控制工具箱提供了直接的求解该方程的lyap()函数:
  1. A=[8 1 6;3 5 7;4 9 2]
  2. B=[2 3;4 5]
  3. C=[1 2;3 4;5 6]
  4. X=lyap(A,B,C)
复制代码

同理,我们使用Kronecker乘机的形式将原方程进行如下变化:
08.gif

故可以编写Sylvester方程的解析解函数:
  1. function X=lyap3(A,B,C)
  2. %Sylvester方程的解析解法
  3. %rewrited by dynamic
  4. %more information
  5. If nargin==2,C=B;B=A';end
  6. [nr,nc]=size(C);
  7. A0=kron(A,eye(nc))+kron(eye(nr),B');
  8. try
  9.     C1=C';
  10.     X0=-inv(A0)*C1(:);
  11.     X=reshape(X0,nc,nr);
  12. catch
  13.     error('Matlabsky提醒您:矩阵奇异!');
  14. end
复制代码

使用上面的数据,我们试验下该解析解法
>>X=lyap3(sym(A),B,C)
X =
[   2853/14186,    557/14186,  -9119/14186]
[  11441/56744,   8817/56744, -50879/56744]

三、Riccati方程的Matlab求解
Riccati方程是一类很著名的二次型矩阵形式,其一般形式为:
09.gif
由于含有矩阵X的二次项,所有Riccati方程求解要Lyapunov方程更难,Matlab控制工具箱提供了are()函数,可以直接求解该函数:
  1. A=[-2 1 -3;-1 0 -2;0 -1 -2]
  2. B=[2 2 -2;-1 5 -2;-1 1 2]
  3. C=[5 -4 4;1 0 4;1 -1 5]
  4. X=are(A,B,C)
复制代码
发表评论

最新评论

引用 zhangzy 2017-1-25 08:43
全都是线性代数的知识  
引用 eastar 2017-2-3 14:27
zhangzy 发表于 2017-1-25 08:43
全都是线性代数的知识

确实  线性代数用的比较多 毕竟计算机处理的已矩阵居多
引用 sh_lin30 2017-2-6 16:19
不错 收藏啦
引用 Raspberry 2017-2-8 08:41
请问楼主如何求解形如Dx=Ax+xA'+xB'Bx+CC'的微分Riccati方程

查看全部评论(4)

返回顶部