Lyapunov方程求解(附件)
广西大学实验报告纸
学院：电气工程学院 专业：自动化 成绩： 组员：陈平忠（1302120238） 黄智榜（1302120237） 班级：
实验地点：808实验室 2015年12月 18日 实验内容：Lyapunov 方程求解
【实验目的】
1、掌握求解Lyapunov 方程的一种方法，了解并使用MATLAB 中相应函数。
【实验设备与软件】
1、硬件：PC 机一台；软件：MATLAB/Simulink 。
【实验原理】
1、线性定常系统渐进稳定的Lyapunov 方程判据
线性定常连续系统为渐进稳定的充要条件是：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程：
Q PA P A T -=+
该条件在传递函数最小实现下等价于：全部特征根都是负实数或实部为负的复数，亦即全部根都位于左半复平面。
线性定常离散系统为渐进稳定的充要条件：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程：
Q P PG G T -=-
该条件在传递函数最小实现下等价于：全部特征根的摸均小于1，即都在单位圆内。
2、在MATLAB 控制工具箱中，函数lyap 和dlyap 用来求解lyapunov 方程。
P =lyap （T A ,Q )可解连续时间系统的lyapunov 方程，其中，Q 和A 为具有相同维数的方阵（A 是系统矩阵）。
如果Q 是对称的，则解P 也是对称的。
P =dlyap （T G ，Q ）可解离散时间系统的lyapunov 方程，其中，Q 和G 为具有相同维数的方阵(G 是系统矩阵)。
如果Q 是对称的，则解P 也是对称的。
3、连续情况下的最小相位系统：系统的零点均在左半复平面，但系统首先是稳定的，其他情况为非最小相位系统。
【实验内容、方法、过程与分析】 题目1实验内容：
输入连续状态空间模型()∑=D C,B,A,:
[]0,110
0,0001,01
00
001000014283==⎥⎥⎥⎥
⎦
⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢
⎢⎣⎡----=D C B A
（1）选取正定矩阵⎥
⎥
⎥⎥⎦
⎤⎢⎢⎢
⎢⎣⎡=100001000010
0001Q ,求稳定性判别矩阵P ，判定系统是否稳定。
Q PA P A T -=+
（2）求线性系统阶跃响应曲线，并判定是否为最小相位系统， （3）求系统的实现，判定是否是最小实现并比较。
题目1实验过程及结果分析：
根据题意，在实验中，先通过运算可以得出结果，根据结果做出如下的.m 文件 程序：
①、由实验.m 文件程序运行后结果：
A=[-3 -8 -2 -4;1 0 0 0;0 1 0 0;0 0 1 0]; B=[1;0;0;0]; C=[0 0 1 1]; D=0;
Q=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]; p=lyap(A',Q) y=ss(A,B,C,D) [V,X]=eig(A) step(y)
得到正定矩阵P ：
②、由题意得出系统的响应曲线：
由图可知：该系统是渐进稳定的。
求特征根
x =
Columns 1 through 2
-1.4737 + 2.2638i 0.0000 + 0.0000i
0.0000 + 0.0000i -1.4737 - 2.2638i
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i
Columns 3 through 4
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i
-0.0263 + 0.7399i 0.0000 + 0.0000i
0.0000 + 0.0000i -0.0263 - 0.7399i
由结果可以得出，此系统特征根的实部全部都为负数，亦全部的根都在左边平面。
所以该系统为最小相位系统。
所以，根据题意，更改A矩阵，求其阶跃响应曲线，并进行比较得：
之前的A矩阵：
更改之前的特征值：
x =
Columns 1 through 2
-1.4737 + 2.2638i 0.0000 + 0.0000i
0.0000 + 0.0000i -1.4737 - 2.2638i
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i Columns 3 through 4
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i -0.0263 + 0.7399i 0.0000 + 0.0000i
0.0000 + 0.0000i -0.0263 - 0.7399i 更改前的阶跃响应：
更改之后的A矩阵：
更改之后的特征值： X =
4.7926 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -1.7297 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0315 + 0.6939i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0315 - 0.6939i 更改后的阶跃响应：
对比特征值可知，更改矩阵A 后特征根有一个为正数，即在右半平面； 对比阶跃响应图可知，更改矩阵A 后，其阶跃响应为发散的。
2、输入离散状态空间模型()∑=D C,H,G,
[]0,100,321,001323031==⎥⎥
⎥⎦
⎤
⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---=D C H G
(1)选定正定矩阵⎥⎥
⎥⎦
⎤
⎢⎢⎢⎣⎡=100010001Q ,求稳定性判别矩阵P 。
(2)请定义离散情况下的最小相位系统。
(3)求线性系统阶跃响应曲线，并按你所定义的判别矩阵是否为最小相位系统。
根据题意，在实验中，先通过运算可以得出结果，根据结果做出如下的c 文件程序：
①、由实验c 文件程序运行后结果： G=[1 3 0;-3 -2 -3;1 0 0]; H=[1;2;3]; C=[0 0 1]; D=0;
Q=[1 0 0;0 1 0;0 0 1]; P=dlyap(G',Q) y=ss(G,H,C,D) [V,X]=eig(G) step(y)
得到矩阵P ：
（2）、请定义离散情况下的相位系统
对于线性定常离散系统，全部特征根的模均小于1，即都在单位圆内，才能认为是最小相位系统。
由CV明显可看出不满足上述条件,且通过图形可知，系统不稳定。
现改变G 的值：
由图可知，阶跃响应最终稳定，满足线性定常离散系统的条件，即极点均为于单位圆内。
【实验总结】
2、通过本次实验了解并掌握了Lyapunov方程的一种用MATLAB求解的方法，并熟悉了线
性定常系统渐进稳定的Lyapunov方程判据和求解lyapunov方程的一些函数。
lyapunov方程的求解
lyapunov方程的求解
听说Lyapunov方程了吗？就是那个让一堆数学大师头疼的东西。
不过别担心，咱们就用大白话聊聊。
说简单点，Lyapunov方程就像是给系统稳定性拍了个“X光”，能看出系统内部的问题。
你想想看，要是你的自行车轮子不稳，骑
起来就得摇摇晃晃，对吧？这就是因为稳定性没搞好。
而Lyapunov
方程就是帮我们找到那个能让系统稳如泰山的“魔法公式”。
话说回来，求解Lyapunov方程可不是件轻松的事儿。
你得有点
数学功底，还得有点耐心和毅力。
有时候，解这个方程就像是解一
个复杂的拼图游戏，得把各个碎片拼在一起，才能看到完整的图画。
不过，好消息是，现在有了电脑和数学软件，求解Lyapunov方
程变得容易多了。
就像是你有了一个超级助手，帮你处理那些繁琐
的计算和推理。
这样一来，你就能更快地找到答案，也不用那么头
疼了。
所以啊，虽然Lyapunov方程听起来有点吓人，但只要咱们用对
方法，就能轻松搞定它。
就像是你面对一个看似复杂的问题，只要找到了解决方法，就能迎刃而解。
这就是数学的魅力所在！。
lyapunov方程求数值解
lyapunov方程求数值解
Lyapunov方程是控制理论中的一个重要方程，用于求解线性系
统的稳定性。
Lyapunov方程的一般形式为AX + XA^T = -Q，其中A
是系统的状态矩阵，X是要求解的对称正定矩阵，Q是一个对称正定
矩阵。
Lyapunov方程的解决对于确定系统的稳定性和性能至关重要。
要求解Lyapunov方程的数值解，通常可以采用以下方法之一：
1. Schur分解法，这是一种常用的数值方法，它将状态矩阵A
分解为一个正交矩阵和一个上三角矩阵的乘积。
然后，可以将Lyapunov方程转化为一个更容易求解的形式，进而求解X的数值解。
2. 离散时间Lyapunov方程的数值解，对于离散时间系统，可
以利用迭代法或者数值线性代数方法来求解Lyapunov方程的数值解。
3. MATLAB等数学软件，许多数学软件包括MATLAB都提供了专
门用于求解Lyapunov方程的函数或工具箱，可以直接利用这些工具
来求解数值解。
无论采用哪种方法，都需要注意数值解的稳定性和精度，尤其
是在系统维度较大时。
此外，还需要对所得到的数值解进行验证，确保其满足Lyapunov方程的定义和性质。
总之，求解Lyapunov方程的数值解需要结合数值方法和专业工具，以确保得到准确可靠的结果。
李雅普诺夫（Lyapunov）函数分析.pdf
讲义81. 李雅普诺夫（Lyapunov ）函数分析本讲中，对于一些有*E (,)0t S r w ⎡⎤=⎣⎦的*γ，我们研究1(,)t t t t t r r S r w γ+=+的收敛性。
回顾一下确定性实例中的Lyapunov 函数分析，我们选取了函数()V r 使得** ()0, ,()()0, , ()0.T V r r V r S r r r V r •≥∀•∇<≠•∇=如收敛性的论证为：我们发现()t V r 随时间减小并且有下限，因此，()t V r 收敛。
对V 和S 采用技术条件，可以证明*t r r →。
现在转到随机实例，用t F 表示到t 时刻的过程历史记录，显然，t F 可表示为{},,,,,,.t l l t r l t w l t l t γ=≤<≤F注意，步长t γ依赖于随机的历史记录，而步依赖于扰动t w 。
定义欧几里德范数122()T V V V =。
定理1 假设V ∃使得（a ）()0, ,V r r ≥∀（b ）L ∃使得22()()V r V r L r r ∇−∇≤−(李普希茨连续Lipschitz continuity) （c ）12,K K ∃使得221222E (,)(),t t t t S r w K K V r ⎡⎤≤+∇⎣⎦F（d ）c ∃使得22()E (,)().T t t t t t V r S r w c V r ∇⎡⎤≤−∇⎣⎦F 则，如t γ满足0t t γ∞==∞∑和20t t γ∞=<∞∑，有z ()t V r 收敛。
z lim ()0t t V r →∞∇=z 每一个t r 的极限点r 满足()0V r ∇=我们将证明某特例的收敛性，该特例对于一些*r 有2*122()V r r r =−。
定理2 假设2*122()V r r r=−满足（a ）12,K K ∃使得2122E (,)(),t t t t S r w K K V r ⎡⎤≤+⎣⎦F（b ）c ∃使得()E (,)().T t t t t t V r S r w cV r ∇⎡⎤≤−⎣⎦F则，如t γ满足0t t γ∞==∞∑和20t t γ∞=<∞∑，有*t r r →, w.p.1（以概率1）我们用下面的上鞅收敛定理证明定理2。
李雅普诺夫方程求解
李雅普诺夫方程求解李雅普诺夫方程是一个非线性偏微分方程，具体形式如下：ut + uux + αuxx = 0其中，u(x,t)为未知函数，α为常数。
它的物理意义是描述一维非粘性流体中的波动行为。
该方程的解析解一般较难求解，但是可以通过一些数值方法进行近似求解。
求解李雅普诺夫方程的一种经典方法是使用有限差分法。
该方法将连续的一维空间离散化成N个点，同时将时间轴也进行离散化，得到一个网格结构。
在这个网格上，我们可以用差分方程来逼近方程的求解。
具体来说，我们可以使用简单的方法，比如向前欧拉方法（即前向差分法）或者向后欧拉方法（即后向差分法），也可以使用更高阶的方法，比如Crank-Nicolson方法。
无论使用什么方法，都需要注意网格的选择。
如果网格太粗，求解结果的精度会降低；如果网格太细，计算时间会增加，同时出现数值不稳定的现象。
通常情况下，我们需要通过试探，确定合适的网格大小。
求解李雅普诺夫方程的另外一种方法是使用数值模拟法。
该方法可以对方程进行更加精细的求解，同时可以考虑更加复杂和现实的情形。
数值模拟法的基本思想是将流体划分成一个个微小的体积元，同时考虑它们之间的相互作用和力的作用。
在这个基础上，我们可以模拟出流体在某一时刻的状态，并利用时间迭代，得到流体在未来各个时刻的状态。
数值模拟法的缺点是计算速度较慢，同时也难以处理特定的边界条件。
但是，它适用于各种不同的物理问题，并且也可以处理更加复杂的流体现象。
总的来说，李雅普诺夫方程是一个非常重要的理论问题。
虽然它的解析解较为复杂，但是通过数值方法和物理模拟，我们可以有效地求解它，同时深入研究一维非粘性流体的波动行为。
Lyapunov方程求解(附件)
广西大学实验报告纸学院：电气工程学院 专业：自动化 成绩： 组员：陈平忠（1302120238） 黄智榜（1302120237） 班级：实验地点：808实验室 2015年12月 18日 实验内容：Lyapunov 方程求解【实验目的】1、掌握求解Lyapunov 方程的一种方法，了解并使用MATLAB 中相应函数。
【实验设备与软件】1、硬件：PC 机一台；软件：MATLAB/Simulink 。
【实验原理】1、线性定常系统渐进稳定的Lyapunov 方程判据线性定常连续系统为渐进稳定的充要条件是：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程：Q PA P A T -=+该条件在传递函数最小实现下等价于：全部特征根都是负实数或实部为负的复数，亦即全部根都位于左半复平面。
线性定常离散系统为渐进稳定的充要条件：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程：Q P PG G T -=-该条件在传递函数最小实现下等价于：全部特征根的摸均小于1，即都在单位圆内。
2、在MATLAB 控制工具箱中，函数lyap 和dlyap 用来求解lyapunov 方程。
P =lyap （T A ,Q )可解连续时间系统的lyapunov 方程，其中，Q 和A 为具有相同维数的方阵（A 是系统矩阵）。
如果Q 是对称的，则解P 也是对称的。
P =dlyap （T G ，Q ）可解离散时间系统的lyapunov 方程，其中，Q 和G 为具有相同维数的方阵(G 是系统矩阵)。
如果Q 是对称的，则解P 也是对称的。
3、连续情况下的最小相位系统：系统的零点均在左半复平面，但系统首先是稳定的，其他情况为非最小相位系统。
【实验内容、方法、过程与分析】 题目1实验内容：输入连续状态空间模型()∑=D C,B,A,:[]0,1100,0001,0100001000014283==⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=D C B A（1）选取正定矩阵⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=1000010000100001Q ,求稳定性判别矩阵P ，判定系统是否稳定。
模型降阶中的Lyapunov方程求解方法
Ｏ ｒ ｅ Ｉ ｅ ｃ ｎｎｅ ｔ ｄ ｒ ｎｔ ｒ ｏ ｃ Ｍ ａ ｒ ｍ ｏ ｌｎｇ ｃｏ ｄｅｉ Ａ ｌ ｏｒｔ ｇ ｉｈｍ ）。 ［］
算 法 降低 Ｌ ａ ｕ ｏ ｙ ｐ ｎ ｖ方 程 求 解 的 计 算 复 杂 性 ， 种 这
方法 既能综 合 上述 方 法 的优 点 ， 兼 顾 了 速度 且 保 又 持 了系统 特性 ．
处 理该 问题 ． 平 衡截 断 法 的关 键 步 骤是 求 解 两 个 Ｌ ａ ｕ ｏ ｙｐ ｎｖ
量 的维度 相匹 配． 系统 的控 制 矩 阵 和 观察 矩 阵分 此
别 是 下 面两个 Ｌ ａ ｕ ｏ ｙ ｐ ｎ ｖ方 程 的解 Ｐ和 Ｑ
｜ Ｐ Ｐ Ａ ＢＢ 一 ０
方程 ， 以此得 到系 统 的控 制 和 观察 矩 阵 ＿ ． 而 ， ４然 ］ 直 接 求解 Ｌ ａ ｕ ｏ ｙ ｐ ｎ ｖ方 程所 需 的计 算 量 将 是难 以承 受 的． 这种 直接 求解 法 最 初 在 控 制领 域 也 仅 仅 用 来 处
文献标识码 ： Ａ
中 图分 类 号 ： ３ ． ； ７ ． Ｏ ２ １ １ ＴＰ ２ １ １
方 程 的低秩 近似 解 ， 使用 Ｐ Ｉ Ｒ ＭＡ（ ａ ｓ ｅＲｅ ｕ ｅ— Ｐ ｓｉ ｄ ｃｄ ｖ
０ 引 言
在 对一 个 复杂 的 线 性 系 统进 行 仿 真 的时 候 ， 为 了减少 计算 时 间和 存 储 量 的需 求 ， 常会 对 原 系 统 通 进 行 降阶处 理 ， 以期 望 得 到 一 个 低 阶 近似 模 型 来 仿
厂 Ａ Ｂ ］ ． 厂Ｔ ， ＡＴ Ｔ ］ Ｂ
Ｄ Ｊ
究 者在 此算 法 中 引入 迭 代 思 想 产 生 了 ＡＤ 算 法 Ｉ 和 Ｌ ｗ— ａ ｋ Ｌ ｏ Ｒ ｎ （ Ｒ）Ｓ ｔ ｍｉ ｈ算 法ｌ ］ 尽 管 这 些 算 法 ６ ．
matlab编写的Lyapunov指数计算程序
% | y x -b |
%
% Then, the variational equation has a form:
%
% F = J*Y
% where Y is a square matrix with the same dimension as J.
% Corresponding m-file:
for l=1:n1 gsc(k)=gsc(k)+y(n1*j+l)*y(n1*k+l); end;
end;
for k=1:n1
for l=1:(j-1)
y(n1*j+k)=y(n1*j+k)-gsc(l)*y(n1*l+k);
end;
end;
znorm(j)=0.0;
for k=1:n1 znorm(j)=znorm(j)+y(n1*j+k)^2; end;
% stept - step on t-variable for Gram-Schmidt renormalization procedure.
% tend - finish value of time
% ystart - start point of trajectory of ODE system.
while calculation_progress == 1
[T,Y] = integrator(DS(1).method_int,@ode_lin,[t tt],y,options,P,n,neq,n_exp);
first_call = 0;
if calculation_progress == 99, break; end;
Matlab方程的求解
Matlab的Lyapunov、Sylvester和Riccati方程的Matlab求解一、连续Lyapunov方程连续Lyapunov方程可以表示为Lyapunov方程来源与微分方程稳定性理论，其中要求C为对称正定的n×n方阵，从而可以证明解X亦为n×n对称矩阵，这类方程直接求解比较困难，不过有了Matlab中lyap()函数，就简单多了。
>> A=[1 2 3;4 5 6;7 8 0]A =1 2 34 5 67 8 0>> C=-[10 5 4;5 6 7;4 7 9]C =-10 -5 -4-5 -6 -7-4 -7 -9>> X=lyap(A,C)X =-3.9444 3.8889 0.38893.8889 -2.7778 0.22220.3889 0.2222 -0.1111二、Lyapunov方程的解析解利用Kroncecker乘积的表示方法，可以将Lyapunov方程写为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]三、离散Lyapunov方程离散Lyapunov方程的一般形式为Matlab中直接提供了dlyap()函数求解该方程：X=dlyap(A,Q)其实，如果A矩阵非奇异，则等式两边同时右乘得到就可以将其变换成连续的Sylvester方程而Sylvester方程是广义Lyapunov方程，故离散的Lyapunov方程还可以使用下面的方法求解B=-inv(A’)C=Q*in v(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矩阵。
离散时间的lyapunov方程
离散时间的lyapunov方程
离散时间的Lyapunov方程在控制理论和动力系统中起着重要作用。
Lyapunov方程是由俄罗斯数学家亚历山大·米哈伊洛维奇·李
亚普诺夫在1892年提出的，用于研究非线性系统的稳定性。
离散时
间的Lyapunov方程描述了动力系统在离散时间下的稳定性和收敛性。
离散时间的Lyapunov方程可以表示为：
\[ V(x_{k+1}) V(x_k) = -\alpha (x_k) \]
其中，\( V(x_k) \)是Lyapunov函数，\( x_k \)是系统在第
k个离散时间点的状态，\( \alpha (x_k) \)是一个非负函数。
Lyapunov函数通常被选为系统状态的二次型函数，以便通过Lyapunov方程来分析系统的稳定性。
通过离散时间的Lyapunov方程，我们可以判断系统在不同状态
下的稳定性。
如果对于所有的状态\( x_k \)，都存在一个非负函数
\( \alpha (x_k) \)，使得Lyapunov方程成立，那么系统在离散时
间下是稳定的。
这种分析方法在控制系统的设计和动力系统的研究
中具有重要的应用价值。
离散时间的Lyapunov方程为我们提供了一种有效的工具，用于研究非线性系统在离散时间下的稳定性和收敛性。
通过对Lyapunov 函数和非负函数的选择，我们可以对系统的稳定性进行定量分析，并设计出有效的控制策略。
因此，离散时间的Lyapunov方程对于控制理论和动力系统的研究具有重要的理论和实际意义。
Lyapunov指数计算
但是在不断地迭代过程中， dF 的值要不断变化。为了表示从整体看两相邻初始 dx
状态分离的情况，必须对时间（或迭代次数）取平均。为此，设平均每次迭代所
引起的指数分离中的指数为 λ ，则初始两点 x0 与 x0 + ε 相距 ε ，经 n 次迭代后，
点变为 F n (x0 ) 与 F n (x0 + ε ) ，相距 εenλ ，即：
n
∑ 成的超球体的体积由 λ1 + λ2 + + λk 决定。而 λi 则决定整个 n 维超球体的体 i =1
积收缩的快慢。
由以上分析可以看出：
（1）任何吸引子，不论其奇怪与否，至少有一个 Lyapunov 指数小于零，否 则轨线就不可能收缩为吸引子。
（2）稳定定态和周期运动（以及准周期运动）不可能有正的 Lyapunov 指数。 稳定定态的 Lyapunov 指数都是负的；周期运动的 λ1 = 0 ，其余的 Lyapunov 指数 也都是负的。
eenλ(x0 ) = F n (x0 + ε ) - F n (x0 ) 取极限 ε → 0 ， n → ∞ ，由上式得到：
（1）