制御工学でシステムの安定性理論等で使用される Lyapunov 方程式を Python で解く方法のメモ。システムが漸近安定であることの必要十分条件は Lyapunov 方程式 \(\mathbf{A}\mathbf{X}+\mathbf{X}\mathbf{A}^{T}+\mathbf{Q}=0\) が任意の正定行列 \(\mathbf{Q}>0\) に対し、唯一の正定解を持つことです。正定行列は対称行列(またはエルミート行列)であるので、全ての固有値が正の実数値という性質を持ちます。

SciPy と Python-Control に Lyapunov 方程式を解く関数があります。 SciPy で解く場合は Lyapunov 方程式が \(\mathbf{A}\mathbf{X}+\mathbf{X}\mathbf{A}^{T}=\mathbf{Q}\) の形式であることに注意が必要です。

なお、より一般的な Lyapunov 方程式である Sylvester 方程式 \(\mathbf{A}\mathbf{X}+\mathbf{X}\mathbf{B}+\mathbf{Q}=0\) も SciPy と Python-Control で解くことができます。

以下の例のように Lyapunov 方程式の解 \(S\) を求め、解が正定行列であることを確認することで、システムが漸近安定であることが確認できます。

# coding: utf-8
import numpy as np
import scipy as sp
import control as ct

A = np.array([[1, 2],
              [-3, -4]])
Q = np.array([[3, 1],
              [1, 1]])

S1 = sp.linalg.solve_lyapunov(A, -Q)
S2 = sp.linalg.solve_sylvester(A, A.conj().transpose(), -Q)
S3 = ct.lyap(A, Q)
S4 = ct.lyap(A, A.conj().transpose(), Q)
print("S1 =\n", S1)
print("\nS2 =\n", S2)
print("\nS3 =\n", S3)
print("\nS4 =\n", S4)

E = sp.linalg.eigvals(S1)
print("\nE =\n", E)

S1 =
 [[ 6.16666667 -3.83333333]
 [-3.83333333  3.        ]]

S2 =
 [[ 6.16666667 -3.83333333]
 [-3.83333333  3.        ]]

S3 =
 [[ 6.16666667 -3.83333333]
 [-3.83333333  3.        ]]

S4 =
 [[ 6.16666667 -3.83333333]
 [-3.83333333  3.        ]]

E =
 [ 8.73078905+0.j  0.43587762+0.j]