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

转载：http://blog.sina.com.cn/s/blog_82662ce70100sh9h.html



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

一、Lyapunov方程

1、连续Lyapunov方程

连续Lyapunov方程可以表示为





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

>> A=[1 2 3;45 6;7 8 0]

C=-[10 5 4;5 6 7;4 7 9]

>> X=lyap(A,C)

2、Lyapunov方程的解析解

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





可见，方程有唯一解的条件并不局限与C对称正定，只要满足非奇异即可保证方程唯一解。同时也打破了传统观念，C必须对称正定的。

function x=lyap2(A,C)

%Lyapunov方程的符号解法

n=size(C,1);

A0=kron(A,eye(n))+kron(eye(n),A);

c=C(:);

x0=-inv(A0)*c;

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方程的一般形式为：





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

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


得到：





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





而Sylvester方程是广义Lyapunov方程，故离散的Lyapunov方程还可以使用下面的方法求解：

B=-inv(A’)

C=Q*inv(A’)

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方程的一般形式为





该方程又称为广义的Lyapunov方程，式中A是n×n方阵，B是m×m方阵，X和C是n×m矩阵。Matlab控制工具箱提供了直接的求解该方程的lyap()函数：

A=[8 1 6;3 5 7;4 9 2]

B=[2 3;4 5]

C=[1 2;3 4;5 6]

X=lyap(A,B,C)

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





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

function X=lyap3(A,B,C)

%Sylvester方程的解析解法

%rewrited by dynamic

%more information

If nargin==2,C=B;B=A';end

[nr,nc]=size(C);

A0=kron(A,eye(nc))+kron(eye(nr),B');

try

C1=C';

X0=-inv(A0)*C1(:);

X=reshape(X0,nc,nr);

catch

error('Matlabsky提醒您：矩阵奇异！');

end

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

>>X=lyap3(sym(A),B,C)

X =

[ 2853/14186, 557/14186, -9119/14186]

[ 11441/56744, 8817/56744,-50879/56744]

三、Riccati方程的Matlab求解

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





由于含有矩阵X的二次项，所有Riccati方程求解要Lyapunov方程更难，Matlab控制工具箱提供了are()函数，可以直接求解该函数：

A=[-2 1 -3;-1 0 -2;0 -1 -2]

B=[2 2 -2;-1 5 -2;-1 1 2]

C=[5 -4 4;1 0 4;1 -1 5]

X=are(A,B,C)