Python で連続時間代数 Riccati 方程式を解く
制御工学の最適レギュレータ問題等で登場する 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]