Python SciPy linalg.solve_discrete_are用法及代碼示例

本文簡要介紹 python 語言中 scipy.linalg.solve_discrete_are 的用法。

用法:

scipy.linalg.solve_discrete_are(a, b, q, r, e=None, s=None, balanced=True)#

求解discrete-time 代數 Riccati 方程 (DARE)。

DARE 定義為

\[A^HXA - X - (A^HXB) (R + B^HXB)^{-1} (B^HXA) + Q = 0\]

存在解決方案的限製是：

• All eigenvalues of \(A\) outside the unit disc, should be controllable.

• The associated symplectic pencil (See Notes), should have eigenvalues sufficiently away from the unit circle.

此外，如果 e 和 s 不都是精確的 None ，那麽 DARE 的廣義版本

\[A^HXA - E^HXE - (A^HXB+S) (R+B^HXB)^{-1} (B^HXA+S^H) + Q = 0\]

解決了。省略時，e 被假定為單位矩陣，s 被假定為零矩陣。

參數 ：：

a： (M, M) 數組

方陣

b： (M, N) 數組

輸入

q： (M, M) 數組

輸入

r： (N, N) 數組

方陣

e： (M, M) 數組，可選

非奇異方陣

s： (M, N) 類數組，可選

輸入

balanced： bool

指示是否對數據執行平衡步驟的布爾值。默認設置為 True。

返回 ：：

x： (M, M) ndarray

離散代數 Riccati 方程的解。

拋出 ：：

LinAlgError

對於無法隔離鉛筆的穩定子空間的情況。有關詳細信息，請參閱注釋部分和引用。

注意：

該方程通過形成擴展辛矩陣鉛筆來求解，如 [1] 中所述，\(H - \lambda J\) 由塊矩陣給出

[  A   0   B ]             [ E   0   B ]
[ -Q  E^H -S ] - \lambda * [ 0  A^H  0 ]
[ S^H  0   R ]             [ 0 -B^H  0 ]

並使用QZ分解方法。

在該算法中，失敗條件與產品 \(U_2 U_1^{-1}\) 的對稱性和 \(U_1\) 的條件編號相關聯。在這裏，\(U\) 是 2m-by-m 矩陣，它包含跨越 2-m 行的穩定子空間的特征向量，並劃分為兩個 m-row 矩陣。有關詳細信息，請參閱 [1] 和 [2]。

為了提高 QZ 分解的準確性，鉛筆經過一個平衡步驟，其中 \(H\) 和 \(J\) 行/列的絕對值之和(在刪除對角線條目後)按照 [3] 中給出的配方進行平衡.如果數據的數值噪聲很小，平衡可能會放大它們的影響，並且需要進行一些清理。

參考：

[1] (1,2)

P. van Dooren，“求解 Riccati 方程的廣義特征值方法。”，SIAM 科學與統計計算雜誌，第 2(2) 卷，DOI:10.1137/0902010

[2]

A.J. Laub，“求解代數 Riccati 方程的 Schur 方法。”，麻省理工學院。信息和決策係統實驗室。 LIDS-R ; 859. 在線可用：http://hdl.handle.net/1721.1/1301

[3]

P. Benner，“哈密頓矩陣的辛平衡”，2001，SIAM J. Sci。計算機, 2001, Vol.22(5), DOI:10.1137/S1064827500367993

例子：

給定 a、b、q 和 r 求解 x：

>>> import numpy as np
>>> from scipy import linalg as la
>>> a = np.array([[0, 1], [0, -1]])
>>> b = np.array([[1, 0], [2, 1]])
>>> q = np.array([[-4, -4], [-4, 7]])
>>> r = np.array([[9, 3], [3, 1]])
>>> x = la.solve_discrete_are(a, b, q, r)
>>> x
array([[-4., -4.],
       [-4.,  7.]])
>>> R = la.solve(r + b.T.dot(x).dot(b), b.T.dot(x).dot(a))
>>> np.allclose(a.T.dot(x).dot(a) - x - a.T.dot(x).dot(b).dot(R), -q)
True