Python SciPy signal.place_poles用法及代码示例

本文简要介绍 python 语言中 scipy.signal.place_poles 的用法。

用法:

scipy.signal.place_poles(A, B, poles, method='YT', rtol=0.001, maxiter=30)#

计算 K 使得特征值 (A - dot(B, K))=poles。

K 是增益矩阵，例如由线性系统 AX+BU 说明的设备将具有其 closed-loop 极点，即特征值 A - B*K ，尽可能接近所要求的极点。

支持 SISO、MISO 和 MIMO 系统。

参数 ：：

A, B： ndarray

线性系统的状态空间表示AX + BU。

poles： array_like

所需的实极点和/或复共轭极点。只有method="YT"(默认)支持复极点。

method: {‘YT’, ‘KNV0’}, optional：

选择哪种方法来找到增益矩阵 K。其中之一：

• ‘YT’: Yang Tits

• ‘KNV0’: Kautsky, Nichols, Van Dooren update method 0

有关算法的详细信息，请参阅引用和注释。

rtol: float, optional：

每次迭代后，特征向量的行列式A - B*K当这两个值之间的相对误差小于rol算法停止。默认值为 1e-3。

maxiter: int, optional：

计算增益矩阵的最大迭代次数。默认值为 30。

返回 ：：

full_state_feedback： 束对象

full_state_feedback 由以下部分组成：

gain_matrix 一维ndarray

闭环矩阵K如A-BK的特征值尽可能接近要求的极点。

computed_poles 一维ndarray

A-BK 对应的极点首先按升序排列为实极点，然后按字典顺序排列为复共轭。

requested_poles 一维ndarray

算法被要求放置的极点如上所示，它们可能与实现的不同。

X 二维数组

转移矩阵如X * diag(poles) = (A - B*K)*X(见注释)

rol 浮点数

实现的相对公差det(X)(见注释)。rol如果可以解决系统问题，将是NaNdiag(poles) = (A - B*K), 或 0 当优化算法不能做任何事情时，即当B.shape[1] == 1.

nb_iter int

收敛前执行的迭代次数。nb_iter如果可以解决系统问题，将是NaNdiag(poles) = (A - B*K), 或 0 当优化算法不能做任何事情时，即当B.shape[1] == 1.

注意：

Tits 和 Yang (YT)，[2] 论文是 Kautsky 等人原始论文的更新。 (KNV)论文[1]。 KNV 依赖 1 级更新来查找传输矩阵 X，使得 X * diag(poles) = (A - B*K)*X ，而 YT 使用 2 级更新。这平均会产生更稳健的解决方案(参见 [2] 第 21-22 页)，此外，YT 算法支持复杂极点，而 KNV 在其原始版本中不支持。这里只实现了 KNV 提出的更新方法 0，因此命名为 'KNV0' 。

在 Matlab 的 place 函数中使用了扩展到复极点的 KNV，YT 由 Slicot 在非自由许可下分发，名称为 robpole。目前尚不清楚 KNV0 如何扩展到复极点(Tits 和 Yang 在其论文的第 14 页上声称他们的方法不能用于将 KNV 扩展到复极点)，因此只有 YT 在此实现中支持它们。

由于极点放置问题的解决方案对于 MIMO 系统来说并不是唯一的，因此这两种方法都从一个试探性的转移矩阵开始，该矩阵以各种方式改变以增加其行列式。两种方法都已被证明可以收敛到一个稳定的解，但是根据选择初始传递矩阵的方式，它们将收敛到不同的解，因此绝对不能保证使用 'KNV0' 会产生类似于 Matlab 或任何其他方法的结果这些算法的实现。

在大多数情况下，使用默认方法'YT' 应该没问题；仅提供'KNV0' 是因为在某些特定情况下'YT' 需要它。此外，当abs(det(X)) 用作稳健性指标时，'YT' 给出的平均结果比'KNV0' 更稳健。

[2] 可通过以下 URL 获取技术报告：https://hdl.handle.net/1903/5598

参考：

[1] (1,2)

J.考茨基，N.K. Nichols 和 P. van Dooren，“线性状态反馈中的鲁棒极点分配”，国际控制杂志，卷。 1985 年，第 41 页，第 1129-1155 页。

[2] (1,2,3)

A.L. Tits 和 Y. Yang，“通过状态反馈进行稳健极点分配的全局收敛算法”，IEEE Transactions on Automatic Control，Vol。 41，第 1432-1452 页，1996。

例子：

一个使用 KNV 和 YT 算法演示实际极点放置的简单示例。这是参考 KNV 出版物 ([1]) 第 4 节中的示例编号 1：

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt

>>> A = np.array([[ 1.380,  -0.2077,  6.715, -5.676  ],
...               [-0.5814, -4.290,   0,      0.6750 ],
...               [ 1.067,   4.273,  -6.654,  5.893  ],
...               [ 0.0480,  4.273,   1.343, -2.104  ]])
>>> B = np.array([[ 0,      5.679 ],
...               [ 1.136,  1.136 ],
...               [ 0,      0,    ],
...               [-3.146,  0     ]])
>>> P = np.array([-0.2, -0.5, -5.0566, -8.6659])

现在使用 KNV 方法 0、默认 YT 方法和 YT 方法计算 K，同时强制算法迭代 100 次，并在每次调用后打印一些结果。

>>> fsf1 = signal.place_poles(A, B, P, method='KNV0')
>>> fsf1.gain_matrix
array([[ 0.20071427, -0.96665799,  0.24066128, -0.10279785],
       [ 0.50587268,  0.57779091,  0.51795763, -0.41991442]])

>>> fsf2 = signal.place_poles(A, B, P)  # uses YT method
>>> fsf2.computed_poles
array([-8.6659, -5.0566, -0.5   , -0.2   ])

>>> fsf3 = signal.place_poles(A, B, P, rtol=-1, maxiter=100)
>>> fsf3.X
array([[ 0.52072442+0.j, -0.08409372+0.j, -0.56847937+0.j,  0.74823657+0.j],
       [-0.04977751+0.j, -0.80872954+0.j,  0.13566234+0.j, -0.29322906+0.j],
       [-0.82266932+0.j, -0.19168026+0.j, -0.56348322+0.j, -0.43815060+0.j],
       [ 0.22267347+0.j,  0.54967577+0.j, -0.58387806+0.j, -0.40271926+0.j]])

X 行列式的绝对值是检查结果稳健性的一个很好的指标，'KNV0' 和'YT' 都旨在最大化它。下面对比上述结果的稳健性：

>>> abs(np.linalg.det(fsf1.X)) < abs(np.linalg.det(fsf2.X))
True
>>> abs(np.linalg.det(fsf2.X)) < abs(np.linalg.det(fsf3.X))
True

现在是复杂极点的简单示例：

>>> A = np.array([[ 0,  7/3.,  0,   0   ],
...               [ 0,   0,    0,  7/9. ],
...               [ 0,   0,    0,   0   ],
...               [ 0,   0,    0,   0   ]])
>>> B = np.array([[ 0,  0 ],
...               [ 0,  0 ],
...               [ 1,  0 ],
...               [ 0,  1 ]])
>>> P = np.array([-3, -1, -2-1j, -2+1j]) / 3.
>>> fsf = signal.place_poles(A, B, P, method='YT')

我们可以在复平面上绘制所需和计算的极点：

>>> t = np.linspace(0, 2*np.pi, 401)
>>> plt.plot(np.cos(t), np.sin(t), 'k--')  # unit circle
>>> plt.plot(fsf.requested_poles.real, fsf.requested_poles.imag,
...          'wo', label='Desired')
>>> plt.plot(fsf.computed_poles.real, fsf.computed_poles.imag, 'bx',
...          label='Placed')
>>> plt.grid()
>>> plt.axis('image')
>>> plt.axis([-1.1, 1.1, -1.1, 1.1])
>>> plt.legend(bbox_to_anchor=(1.05, 1), loc=2, numpoints=1)