李雅普诺夫方程的新解法

塞拉斯, 宋扬, 孙华飞

(北京理工大学 数学与统计学院, 北京 100081)

众所周知, 一个不稳定的系统在工程上是无法使用的, 稳定性是控制理论核心的研究内容之一. 李雅普诺夫方程在线性系统的稳定性研究中具有重要的作用. 通常, 李雅普诺夫方程没有解析解, 只有数值解. 如利用定义在正定矩阵流形上的仿射黎曼度量, 利用测地距离作为目标函数来求解李亚普诺夫方程[1-4]. 本文利用建立在对称正定矩阵上的对数欧氏度量诱导的测地距离作为目标函数, 给出求解李雅普诺夫方程的几何求解方法.

设SPD(n)表示对称正定矩阵的全体, 可以验证它是一般线性群GL(n,R)的子流形. 对称正定流形拥有良好的几何性质, 在图像处理、信号处理以及与统计相关的领域具有广泛的应用. 如果把对称正定矩阵仅仅看成矩阵来对待, 也就是说仅利用它的代数性质来使用未免有些浪费, 而在它们上面再赋予几何的性质使之成为流形, 就可以充分挖掘它的价值. 在SPD(n)上可以定义不同的黎曼度量, 使之成为黎曼度量, 但是如何针对实际问题的需求选择适当的黎曼度量对解决问题至关重要.

对于时不变的线性系统

(1)

u(x(t))=xT(t)Px(t),

(2)

ATP+PA+Q=0.

(3)

文中作者利用几何的方法给出方程式(3)的求解过程. 基本想法是在SPD(n)上建立一个以测地距离为目标函数的测量SPD(n)上两点间最短距离的函数, 当距离函数趋近于0时便获得了求解的结果.

1 称正定矩阵流形的黎曼几何结构

在本节中，简要介绍定义在对称正定矩阵流形SPD(n)上的各种黎曼度量，使得SPD(n)成为黎曼流形. 可以在SPD(n)定义不同的度量, 因而其所呈现的几何结构也不同. 例如定义欧氏度量

〈X,Y〉A=tr(XTY),

(4)

其中tr表示矩阵的迹. 对于SPD(n)上的任意两点A,B, 距离为

(5)

但是, 这样定义的内积所获得的测地线是一个直线γ(t)=Α+(B-A)t, 这里的参数t的范围有限制, 换句话说该直线不能连接SPD(n)上的任意两点. 这样的度量不能用于计算SPD(n)上任意两点的距离.

为了克服上面度量的弱点, 可以定义所谓的仿射的黎曼度量[5-6]. 对于SPD(n)的切空间TASPD(n)上任意的两点X,Y,A∈SPD(n), 定义

〈X,Y〉A=tr(A-1XA-1Y),

(6)

在该度量下的测地线可以表示为

γ(t)=A1/2exp(tA-1/2SA-1/2)A1/2,

(7)

式(7)表示过点A,切向量为S的曲线,参数t没有限制,γ(t)为连接SPD(n)上任意两点的测地线. 利用该测地线, 可以得到连接SPD(n)上任意两点A,B的测地距离为

(8)

其中λi为矩阵A-1B的特征值, 它们是恒正的. 利用式(8), 也可以对式(3)求解, 但是利用上述方法在计算上复杂度比较高, 原因是在仿射黎曼度量下的黎曼流形SPD(n)是一个截面曲率为非正的弯曲空间.

为了降低计算的复杂度, 在SPD(n)上定义新的乘法

A·B=exp[log(A)+log(B)],

(9)

式中A,B是SPD(n)中的两个矩阵. 以验证在该乘法定义下, SPD(n)是一个群, 群结构和拓扑结构是相容的, 从而SPD(n)是一个李群, 而且是一个可交换的李群. 进一步,SPD(n)上存在双不变度量, 使得SPD(n)与欧氏空间TASPD(n)等距. 这样做的意义在于, 尽管SPD(n)本身不是一般线性群GL(n,R)的李子群, 但是在这样的乘法定义下成为李群, 而且在等距意义下保曲率的性质使得SPD(n)成为一个平坦的欧氏空间.

在SPD(n)上定义新的黎曼度量--对数欧氏度量

〈X,Y〉A=〈dlogAX,dlogAY〉I,

(10)

式中:d表示微分;I为单位矩阵. 经计算可得知, 连接SPD(n)上任意两点A,B的测地线可以表示为

γ(t)=exp[(1-t)logA+tlogB],

(11)

此时连接A,B两点的测地距离函数为

dL(A,B)=‖logA-logB‖.

(12)

对于给定的对称正定矩阵Q, 利用式(12)有

dL[Q,-(ATP+PA)]=‖logQ-

log(-(ATP+PA))‖,

(13)

2 对称正定矩阵流形上的梯度算法[7-11]

对于一般的光滑黎曼流形(M,g),以及定义在A上的函数f:M→R,可以利用自然梯度给出求解函数f(θ)的最小值的迭代公式

(14)

利用黎曼度量式(6), 对于建立在SPD(n)上的目标函数fR:SPD(n)→R, 其中

fR(P)=dR[Q,-(ATP+PA)],

(15)

可以获得求目标函数fR的最小值点的迭代公式

(16)

利用欧式内积(4), 对于建立在SPD(n)上的目标函数fE:SPD(n)→R, 其中

fE(P)=dE(Q,-(ATP+PA)),

(17)

可以获得求目标函数fE的最小值点的迭代公式

(18)

根据SPD(n)上的对数欧式度量(10), 建立目标函数

fL(P)=dL(Q,-(ATP+PA)),

(19)

当式(19)趋近于0时, 就可以获得式(3)的解. 希望求该目标函数的最小值来获得正定矩阵P. 利用梯度算法可以得到计算式(19)的迭代公式

(20)

3 模拟仿真

在本节中，将利用例子，通过模拟仿真比较LGDA法来验证的算法的优越性.

对于给定的对称正定矩阵

以及初始矩阵P0=0.5I, 用LGDA、NGDA、GDA法分别求解方程(3), 误差矩阵的模长小于0.01时算法停止. LGDA法求解得到的对称正定矩阵为

迭代次数为3步. NGDA法求解得到的对称正定矩阵为

迭代次数为7步. GDA法不收敛,通过图1,发现

相比于GDA法, NGDA和LGDA有更好的收敛性, 并且LDGA的收敛速度比NGDA更快.

4 结 论

利用定义在对称正定矩阵流形上的欧氏对数度量, 给出了李雅普诺夫矩阵的几何上的求解算法, 模拟仿真结果验证了该方法的优越性. 进一步的研究将着眼于两个方面, 一方面是利用该方法给出代数Riccati方程的求解方法, 研究最优控制问题; 另一方面是利用建立在对称正定矩阵上的主丛结构, 给出新的测地距离函数, 来求解李雅普诺夫方程以及Riccati方程, 获得更好的求解方法; 进一步, 将尝试研究线性时变系统的李雅普诺夫方程和Riccati方程的求解问题, 并研究相关的系统的稳定性和最优控制问题.