如何用python-control（IVP问题）模拟系统传递函数的时间响应？

Posted 2023-03-14

【中文标题】如何用python-control（IVP问题）模拟系统传递函数的时间响应？

【英文标题】：How to simulate the time response of a system transfer function with python-control (IVP problem)?

【发布时间】：2021-02-12 00:00:40

【问题描述】：

我试图演示如何使用系统传递函数的定义和python-control 模块来“求解”（模拟解）微分方程初值问题 (IVP)。事实上，我在控制方面真的是个新手。

我以这个简单的微分为例：y'' - 4y' + 13y = 0，具有以下初始条件：y(0) = 1 和 y'(0) = 0。

我手动实现了这个传递函数： Y(s) = (s - 4)/(s^2 - 4*s + 13).

所以，在 Python 中，我正在编写这段小代码（请注意，y_ans 是这个差分 IVP 的答案为 seen here）：

import numpy as np
import control as ctl
import matplotlib.pyplot as plt

t = np.linspace(0., 1.5, 100)
sys = ctl.tf([1.,-4.],[1.,-4.,13.])
T, yout, _ = ctl.forced_response(sys, T=t, X0=[1, 0])

y_ans = lambda x: 1/3*np.exp(2*x)*(3*np.cos(3*x) - 2*np.sin(3*x))

plt.plot(t, y_ans(t), '-.', color='gray', alpha=0.5, linewidth=3, label='correct answer')
plt.plot(T, yout, 'r', label='simulated')
plt.legend()

这段代码让我得到这张图：

但是当我在yout 前面插入一个负号时，我得到了一个我想要的匹配：

plt.plot(T, -yout, 'r', label='simulated') ### yout with negative sign

我做错了什么？ Python 控制文档对我来说不是很清楚。另外，我不知道我对control.forced_response 的X0 参数的解释是否正确。是否有可能按照我的意图做到这一点？

欢迎任何对这个主题有所了解的人投稿。

编辑

设置X0 = [0,0] 给了我这个图表：

【问题讨论】：

似乎您在forced_response 中对初始条件进行了两次编码，一次在传递函数的分子中，然后在X0 中再次编码。尝试设置X0=[0,0]。您可能已经为(s-4)^2/(s^2-4s+13) 实施了解决方案。

谢谢，@LutzLehmann，但是当我设置X0=[0,0] 或者如果我让它设置为默认值时，结果是横坐标上的一条归零的平线（我在问题中包含了图表，编辑部分)。

【参考方案1】：

我认为这里最好的做法是将您的系统转换为状态空间并看看发生了什么（有许多可能的状态空间表示）：

sys_tf = ctl.tf([1.,-4.],[1.,-4.,13.])
sys_ss = ctl.tf2ss(sys_tf)
print(sys_ss)

输出：

A = [[ 4.00000000e+00  1.30000000e+00]
     [-1.00000000e+01  1.33226763e-15]]

B = [[-1.]
     [ 0.]]

C = [[-1.  -0.4]]

D = [[0.]]

我们要找到x(0)，这样y(0) = Cx(0) = 1 和y'(0) = CAx(0) = 0。

我们可以写出这些方程并手动求解，也可以使用线性代数：

A = np.vstack([sys_ss.C, sys_ss.C @ sys_ss.A])
b = np.array([[1], [0]])
x0 = np.linalg.solve(A, b)
print(x0)

给予：

[[-1.00000000e+00]
 [-1.36642834e-15]]

因此，这应该有效：

T, yout = ctl.forced_response(sys_ss, T=t, X0=[-1, 0])

此外，由于您只对初始条件（即u(t)=0）的瞬态响应感兴趣，您可以使用initial_response 函数：

T, yout = ctl.initial_response(sys_ss, T=t, X0=[-1, 0])

【讨论】：

如果检查最新版本 0.9.0 中的原始代码，则会发现问题中的行为不再存在。

确实如此，尽管我能够通过对代码稍作更改来复制问题中的结果：T, yout = ctl.forced_response(sys, T=t, X0=[1, 0])。因此，我认为我的回答仍然是相关的。

【参考方案2】：

感谢@LutzLehmann 的评论，我一直在思考“编码两次”的含义。所以，回到第一方，我意识到这个传递函数包含输入（时间或斜坡）和初始条件。它实际上是一个输出。我需要某种逆拉普拉斯变换，或者，正如我开始认为的那样，我只需要按原样模拟它，无需进一步的信息。

因此，我设法使用了脉冲输入（拉普拉斯变换等于 1），并且能够得到与我的 tf 及时模拟的输出。

import numpy as np
import control as ctl
import matplotlib.pyplot as plt

t = np.linspace(0., 1.5, 100)
sys = ctl.tf([1.,-4.],[1.,-4.,13.])
T, yout = ctl.impulse_response(sys, T=t) # HERE is what I wanted

y_ans = lambda x: 1/3*np.exp(2*x)*(3*np.cos(3*x) - 2*np.sin(3*x))

plt.plot(t, y_ans(t), '-.', color='gray', alpha=0.5, linewidth=3, label='correct answer')
plt.plot(T, yout, 'r', label='simulated')
plt.legend()

现在我想我可以展示如何使用 python-control 来间接模拟微分方程的答案。 :-D