制御工学の最適レギュレータ問題等で登場する Riccati 方程式を Python で解く方法のメモ。 Riccati 方程式の解は SciPy の solve_continuous_are 関数で求められます。ただし、この方法では解を求めるだけなので、状態フィードバックゲインを求めたり閉ループ系の極を確認するのは別途計算が必要です。

状態フィードバックゲインを求めたり、閉ループ系の極まで求めたい場合は、Python-Control の care か lqr 関数を使用した方が簡単です。それぞれの関数は出力変数の順序が異なるのが注意点です。

以下の例では Riccati 方程式の解 $S$、状態フィードバックゲイン $K$、閉ループ系の極 $E$ をそれぞれの関数で求めています。

import numpy as np
import scipy as sp
import control as ct

A = np.array([[0, 1],
              [0, -1]])
B = np.array([[0],
              [1]])
Q = np.diag([2, 1])
R = np.array([[1]])

S1 = sp.linalg.solve_continuous_are(A, B, Q, R)
K1 = np.linalg.inv(R).dot(B.T).dot(S1)
E1 = np.linalg.eigvals(A-B.dot(K1))
print("S1 =\n", S1)
print("K1 =\n", K1)
print("E1 =\n", E1)

S2, E2, K2 = ct.care(A, B, Q, R)
print("\nS2 =\n", S2)
print("K2 =\n", K2)
print("E2 =\n", E2)

K3, S3, E3 = ct.matlab.lqr(A, B, Q, R)
print("\nS3 =\n", S3)
print("K3 =\n", K3)
print("E3 =\n", E3)

S1 =
 [[ 3.10754795  1.41421356]
 [ 1.41421356  1.19736823]]
K1 =
 [[ 1.41421356  1.19736823]]
E1 =
 [-1.09868411+0.45508986j -1.09868411-0.45508986j]

S2 =
 [[ 3.10754795  1.41421356]
 [ 1.41421356  1.19736823]]
K2 =
 [[ 1.41421356  1.19736823]]
E2 =
 [-1.09868407+0.45508987j -1.09868407-0.45508987j]

S3 =
 [[ 3.10754795  1.41421356]
 [ 1.41421356  1.19736823]]
K3 =
 [[ 1.41421356  1.19736823]]
E3 =
 [-1.09868407+0.45508987j -1.09868407-0.45508987j]