solve_discrete_are#
- scipy.linalg.solve_discrete_are(a, b, q, r, e=None, s=None, balanced=True)[源代码]#
求解离散时间代数 Riccati 方程 (DARE)。
DARE 的定义如下:
\[A^HXA - X - (A^HXB) (R + B^HXB)^{-1} (B^HXA) + Q = 0\]解存在的限制条件是:
单位圆盘外的所有 \(A\) 的特征值都应是可控的。
相关的辛矩阵束(参见“备注”)的特征值应远离单位圆。
此外,如果
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布尔值
一个布尔值,指示是否对数据执行平衡步骤。默认设置为 True。
- 返回:
- x(M, M) ndarray
离散代数 Riccati 方程的解。
- 引发:
- LinAlgError
对于无法分离矩阵束的稳定子空间的情况。详细信息请参见“备注”部分和参考文献。
另请参见
solve_continuous_are
求解连续代数 Riccati 方程
备注
该方程通过构造扩展辛矩阵束来求解,如 [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 乘 m 的矩阵,包含张成稳定子空间的特征向量,具有 2m 行,并被划分为两个 m 行矩阵。更多详细信息请参见 [1] 和 [2]。
为了提高 QZ 分解的精度,该矩阵束会经过一个平衡步骤,其中 \(H\) 和 \(J\) 的行/列的绝对值之和(去除对角线元素后)按照 [3] 中给出的方法进行平衡。如果数据存在少量数值噪声,平衡可能会放大其影响,因此需要进行一些清理。
0.11.0 版本新增。
参考文献
[1] (1,2)P. van Dooren , “A Generalized Eigenvalue Approach For Solving Riccati Equations.”, SIAM Journal on Scientific and Statistical Computing, Vol.2(2), DOI:10.1137/0902010
[2]A.J. Laub, “A Schur Method for Solving Algebraic Riccati Equations.”, Massachusetts Institute of Technology. Laboratory for Information and Decision Systems. LIDS-R ; 859. 在线获取 : http://hdl.handle.net/1721.1/1301
[3]P. Benner, “Symplectic Balancing of Hamiltonian Matrices”, 2001, SIAM J. Sci. Comput., 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