【原创】Lyapunov、Sylvester和Riccati方程的Matlab求解

dynamic

Lyapunov、Sylvester和Riccati方程是控系统常用到的几个方程，应用和计算比较广泛


在这里我们只要讨论下Lypunov方程的连续方程、离散方程的数值和解析解法，其中数值解法Matlab提供的直接的lyap()和dlyap()函数

一、连续Lyapunov方程

连续Lyapunov方程可以表示为

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

1. >> A=[1 2 3;4 5 6;7 8 0]
   
2. 
   
3. A =
   
4. 
   
5.      1     2     3
   
6.      4     5     6
   
7.      7     8     0
   
8. 
   
9. >> C=-[10 5 4;5 6 7;4 7 9]
   
10. 
    
11. C =
    
12. 
    
13.    -10    -5    -4
    
14.     -5    -6    -7
    
15.     -4    -7    -9
    
16. 
    
17. >> X=lyap(A,C)
    
18. 
    
19. X =
    
20. 
    
21.    -3.9444    3.8889    0.3889
    
22.     3.8889   -2.7778    0.2222
    
23. 0.3889    0.2222   -0.1111

复制代码

二、Lyapunov方程的解析解

利用Kroncecker乘积的表示方法，可以将Lyapunov方程写为

可见，方程有唯一解的条件并不局限与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)

复制代码

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

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

复制代码

三、离散Lyapunov方程

离散Lyapunov方程的一般形式为

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

其实，如果A矩阵非奇异，则等式两边同时右乘

得到

就可以将其变换成连续的Sylvester方程

而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方程

dynamic

Sylvester方程的一般形式为

该方程又称为广义的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)
   
5. 
   
6. 
   
7. A =
   
8. 
   
9.      8     1     6
   
10.      3     5     7
    
11.      4     9     2
    
12. 
    
13. 
    
14. B =
    
15. 
    
16.      2     3
    
17.      4     5
    
18. 
    
19. 
    
20. C =
    
21. 
    
22.      1     2
    
23.      3     4
    
24.      5     6
    
25. 
    
26. 
    
27. X =
    
28. 
    
29.     0.2011    0.2016
    
30.     0.0393    0.1554
    
31.    -0.6428   -0.8966

复制代码

同理，我们使用Kronecker乘机的形式将原方程进行如下变化

故可以编写Sylvester方程的解析解函数

1. function X=lyap3(A,B,C)
   
2. %Sylvester方程的解析解法
   
3. %rewrited by dynamic
   
4. %more information http://www.matlabsky.cn
   
5. 
   
6. If nargin==2,C=B;B=A';end
   
7. [nr,nc]=size(C);
   
8. A0=kron(A,eye(nc))+kron(eye(nr),B');
   
9. try
   
10.     C1=C';
    
11.     X0=-inv(A0)*C1(:);
    
12.     X=reshape(X0,nc,nr);
    
13. catch
    
14.     error('Matlabsky提醒您：矩阵奇异！');
    
15. end

复制代码

恩，使用上面的数据，我们试验下该解析解法的

1. >>X=lyap3(sym(A),B,C)
   
2. 
   
3. X =
   
4. 
   
5. [   2853/14186,    557/14186,  -9119/14186]
   
6. [  11441/56744,   8817/56744, -50879/56744]

复制代码

dynamic

Riccati方程是一类很著名的二次型矩阵形式，其一般形式为

由于含有矩阵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)
   
5. 
   
6. 
   
7. A =
   
8. 
   
9.     -2     1    -3
   
10.     -1     0    -2
    
11.      0    -1    -2
    
12. 
    
13. 
    
14. B =
    
15. 
    
16.      2     2    -2
    
17.     -1     5    -2
    
18.     -1     1     2
    
19. 
    
20. 
    
21. C =
    
22. 
    
23.      5    -4     4
    
24.      1     0     4
    
25.      1    -1     5
    
26. 
    
27. 
    
28. X =
    
29. 
    
30.     0.9874   -0.7983    0.4189
    
31.     0.5774   -0.1308    0.5775
    
32.    -0.2840   -0.0730    0.6924

复制代码

kou1714js

感謝樓主
這個在適應性控制常常用到

yyuu

我是需要分析计算LOGISTIC数据  ,都是用来说明对初值的敏感.以下是LOGISTIC求解的程序 ,希望得到LYAPUNOV的程序.
clc
clear
close all

lambda = 3:5e-4:4;
x = 0.4*ones(1,length(lambda));

N1 = 400;                   % 前面的迭代点数
N2 = 100;                   % 后面的迭代点数

f = zeros(N1+N2,length(lambda));
for i = 1:N1+N2
    x = lambda .* x .* (1 - x);
    f(i,:) = x;
end
f = f(N1+1:end,:);
plot(lambda,f,'k.','MarkerSize',1)
xlabel('\lambda')
ylabel('x');

青山不改

回复 5# yyuu


    您有没有作过Lyapunov指数图？我也关心这个问题啊

青山不改

回复 3# dynamic


    感谢您的介绍！

xiaoyufangfei

离散时间代数riccati方程的求解
                                   -1
E'XE = A'XA - (A'XB + S)(B'XB + R)    (A'XB + S)' + Q
用dare求解riccati方程的X的代码是多少？
谢谢！

xiaoyufangfei

A=[-0.2 0 0 0,0 -0.25 0 0,0 0 -0.25 0,0 0 0 -0.25]
B=[12 -17,24 -9,-12 18,-24 10]
C=[1 0 1 0,0 1 0 1]

dynamic1

请问楼主如何求解形如Dx=Ax+xA'+xB'Bx+CC'的微分Riccati方程

注册这个号就是为了请教楼主，务必回答啊，拜谢了！