3.6.3 连续时间线性系统
在上一节中，我们使用odeint( )对质量-弹簧-阻尼系统的微分方程组进行了数值积分，并且进行了PID控制模拟。该系统的微分方程为：
。通过拉普拉斯变换可以将微分方程化为容易求解的代数方程：ms2 X(s)+bsX(s)+kX(s)=F(s)。其中F(s)是f(t)的拉普拉斯变换，X(s)是x(t)的拉普拉数变换，而n次微分变成了sn。F(s)是输入信号，而X(s)是输出信号，将等式改写为输入除以输出的形式，就得到了系统的传递函数P(s)：
连续时间系统的传递函数是两个s的多项式的商。通过连续时间系统的传递函数，很容易计算某输入信号对应的输出信号。在下面的例子中，使用signal模块计算质量-弹簧-阻尼系统对阶跃信号以及正弦波信号的响应输出。❶创建lti对象，可以使用控制理论中的多种形式表示连续时间线性系统，这里使用的是传递函数分子和分母多项式的系数。多项式的系数与numpy.poly1d的约定相同，即下标为0的元素是最高次项的系数。❷调用lti.step( )方法计算系统的阶跃响应。T参数为计算响应的时间数组。❸调用signal.lsim( )计算系统对正弦波信号的响应，它的第一个参数为lti对象，也可以直接传递(numerator, denominator)。U参数为保存输入信号的数组。step( )和lsim( )计算结果中的第二项为系统的输出信号，这里忽略其余的输出。
图3-33显示阶跃响应最终稳定在x=0.05处，这时的kx=1。
图3-33 系统的阶跃响应和正弦波响应
m, b, k = 1.0, 10, 20
numerator = [1]
denominator = [m, b, k]
plant = signal.lti(numerator, denominator) ❶
t = np.arange(0, 2, 0.01)
_, x_step = plant.step(T=t) ❷
_, x_sin, _ = signal.lsim(plant, U=np.sin(np.pi*
t), T=t) ❸
传递函数的代数运算可以表示由多个连续时间系统组成的系统，例如两个系统的级联的传递函数是各个系统的传递函数的乘积。而传递函数由分子和分母两个多项式构成，因此传递函数的四则运算可以使用NumPy的poly1d相关的函数实现。下面的SYS类通过定义__mul__、__add__、__sadd__、__div__等魔法方法，让它支持四则运算。
❶feedback( )方法计算与之对应的反馈系统的传递函数。在图3-34中，P是被控制的系统，C是控制器，C的输入信号是目标信号与实际输入的差xr -x。我们从xr到x的传递函数就是这个反馈系统的传递函数。根据图示可以列出如下拉普拉斯变换之后的代数方程：
图3-34 反馈控制系统框图
整理可得：
如果将C(s)·P(s)看作系统Y(s)，那么可以得出反馈系统的传递函数为：
。
❷为了让SYS对象能作为step( )、lsim( )等函数的第一个表示系统的参数，需要定义__iter__( )魔法方法返回传递函数的分子与分母的多项式系数。
from numbers import Real
def as_sys(s):
if isinstance(s, Real):
return SYS([s], [1])
return s
class SYS(object):
def __init__(self, num, den):
self.num = num
self.den = den
def feedback(self): ❶
return self / (self + 1)
def __mul__(self, s):
s = as_sys(s)
num = np.polymul(self.num, s.num)
den = np.polymul(self.den, s.den)
return SYS(num, den)
def __add__(self, s):
s = as_sys(s)
den = np.polymul(self.den, s.den)
num = np.polyadd(np.polymul(self.num, s.den),
np.polymul(s.num, self.den))
return SYS(num, den)
def __sadd__(self, s):
return self + s
def __div__(self, s):
s = as_sys(s)
return self *
SYS(s.den, s.num)
def __iter__(self): ❷
return iter((self.num, self.den))
下面我们用SYS类计算使用PI控制器控制质量-弹簧-阻尼系统时的阶跃响应。PI控制器的传递函数为：
注意上节中介绍的PI控制器是离散时间的，使用累加器近似计算积分器的输出，而本节采用连续时间系统的系统响应模拟控制系统。
❶质量-弹簧-阻尼系统的传递函数为plant，❷PI控制器的传递函数为pi_ctrl，为了step( )不抛出LinAlgError异常，这里将PI控制器的传递函数的分母常数项设置为一个非常小的值。❸计算反馈系统的传递函数feedback。由图3-35可以看出Ki为0时，系统的输出位移与目标位移之间存在一定的差距，Kp越大差距越小，但是会出现过冲现象。适当调节Kp与Ki可以减弱过冲现象，但是仍然会有超过目标位移的时刻。
图3-35 使用PI控制器的控制系统的阶跃响应
M, b, k = 1.0, 10, 20
plant = SYS([1], [M, b, k]) ❶
pi_settings = [(10, 1e-10), (200, 1e-10),
(200, 100), (50, 100)]
fig, ax = pl.subplots(figsize=(8, 3))
for pi_setting in pi_settings:
pi_ctrl = SYS(pi_setting, [1, 1e-6]) ❷
feedback = (pi_ctrl *
plant).feedback( ) ❸
_, x = signal.step(feedback, T=t)
label = "$K_p={:d}, K_i={:3.0f}$".format(*
pi_setting)
ax.plot(t, x, label=label)
ax.legend(loc="best", ncol=2)
ax.set_xlabel(u"时间(秒)")
ax.set_ylabel(u"位移(米)")
为了计算施加于质量的控制力，可以将误差信号传递给lsim( )计算控制器的输出：
_, f, _ = signal.lsim(pi_ctrl, U=1-x, T=t)
为了彻底消除过冲现象，需要使用PID控制，PID控制器的传递函数为：
下面计算PID控制器构成的反馈系统的阶跃响应。由于PID控制器需要对输入信号进行微分，而阶跃输入信号会导致PID的输出中包含脉冲输出，即时间无限短、值无限大的信号。
kd, kp, ki = 30, 200, 400
pid_ctrl = SYS([kd, kp, ki], [1, 1e-6])
feedback = (pid_ctrl *
plant).feedback( )
_, x2 = signal.step(feedback, T=t)
为了让PID控制器的输出在限定的范围之内，可以在反馈系统之前添加一个低通滤波器，一阶低通滤波器的传递函数为：
。添加低通滤波器之后，PID控制器的输入就是连续信号了，如图3-36所示。
图3-36 带低通滤波器的反馈控制系统框图
lp = SYS([1], [0.2, 1])
lp_feedback = lp *
(pid_ctrl *
plant).feedback( )
_, x3 = signal.step(lp_feedback, T=t)
由于PID控制器的传递函数的分子阶数高于分母阶数，因此无法使用lsim( )计算。我们可以把x_r当作系统的输入，把f当作输出，通过下面的方程计算从x_r到f的传递函数：
F(s)=(Xr (s)·LP(s)-F(s)·P(s))·C(s)
得到的传递函数为：
下面根据上面的公式计算带低通滤波器的控制系统中控制器的输出：
pid_out = (pid_ctrl *
lp) / (pid_ctrl *
plant + 1)
_, f3 = signal.step(pid_out, T=t)
图3-37显示了上述PI控制、PID控制以及带低通滤波的PID控制等系统中滑块的位移以及控制力。由于PID控制的控制力存在脉冲信号，因此无法在图中正确显示。由位移曲线可以看出低通+PID控制可以有效抑制过冲现象。
图3-37 滑块的位移以及控制力