Python で連続時間 Lyapunov 方程式を解く
制御工学でシステムの安定性理論等で使用される 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]