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