3.5.3 ode类
使用odeint( )可以很方便地计算微分方程组的数值解，只需调用一次odeint( )就能计算出一组时间点上的系统状态。但是有时我们希望一次只前进一个时间片段，从而对求解对象进行更精确的控制。在下面的例子中，我们通过ode类模拟如图3-26所示的弹簧系统每隔1ms周期的系统状态，并通过PID控制器控制滑块的位置。
图3-26 质量-弹簧-阻尼系统
该系统的微分方程为：
。其中x为滑块的位移，x ̈为位移对时间的二次导数，即滑块的加速度，x ̇为滑块的速度，m为滑块的质量，b为阻尼系数，k为弹簧的系数，F为外部施加于滑块的控制力。这是一个二次微分方程，为了使用ode对系统求解，需要将其转换成如下一阶微分方程组：
其中x为滑块的位移，u为滑块的速度。这两个变量构成了系统的状态，它们对时间的导数可以通过这两个方程直接算出。
def mass_spring_damper(xu, t, m, k, b, F):
x, u = xu.tolist( )
dx = u
du = (F - k*
x - b*
u)/m
return dx, du
下面使用odeint( )对该系统进行求解，初值为滑块在位移0.0处，起始速度为0，外部控制力恒为1.0。如图3-27所示，系统经过约两秒钟，最终停在了位移0.05米处。
图3-27 滑块的速度和位移曲线
m, b, k, F = 1.0, 10.0, 20.0, 1.0
init_status = 0.0, 0.0
args = m, k, b, F
t = np.arange(0, 2, 0.01)
result = odeint(mass_spring_damper, init_status, t, args)
我们希望通过控制外力F，使得滑块更迅速地停在位移1.0处，这时可以使用PID控制器进行控制。在介绍PID控制器之前，首先让我们用ode类重写odeint( )模拟的部分。ode类和odeint( )一样，也需要一个计算各个状态的导数的函数。
❶这里使用MassSpringDamper类的方法f( )计算状态点处的导数。注意该方法的参数顺序和odeint( )所需的函数mass_spring_damper( )不同，第一个参数为时间，第二个参数为系统状态。并且该方法必须返回一个列表来表示各个状态的导数，不能返回元组。这里使用MassSpringDamper类将系统的各个参数M、k、b以及F等包装在对象内部。
❷创建ode对象之后，通过set_integrator( )设置积分器相关的参数。它的第一个参数为积分器的算法，其后的关键字参数设置该积分器算法的各个参数。关于各个参数的具体含义请读者参阅ode类的文档说明。然后调用set_initial_value( )设置系统的初始状态和初始时间。
❸在while循环中，以dt为间隔对系统进行积分求解。ode对象的属性r.t为当前的模拟时间，调用r.integrate(r.t + dt)计算r.t + dt处的状态。系统的状态保存在r.y中。我们用两个列表t和result分别保存模拟时间和系统状态。
由allclose( )的比较结果可知使用ode的结果与odeint( )的完全相同。
from scipy.integrate import ode
class MassSpringDamper(object): ❶
def __init__(self, m, k, b, F):
self.m, self.k, self.b, self.F = m, k, b, F
def f(self, t, xu):
x, u = xu.tolist( )
dx = u
du = (self.F - self.k*
x - self.b*
u)/self.m
return [dx, du]
system = MassSpringDamper(m=m, k=k, b=b, F=F)
init_status = 0.0, 0.0
dt = 0.01
r = ode(system.f) ❷
r.set_integrator('vode', method='bdf')
r.set_initial_value(init_status, 0)
t = [ ]
result2 = [init_status]
while r.successful( ) and r.t + dt < 2: ❸
r.integrate(r.t + dt)
t.append(r.t)
result2.append(r.y)
result2 = np.array(result2)
np.allclose(result, result2)
True
由于在while循环中我们逐点对系统的状态进行模拟，因此可以在其中增加控制代码，改变作用于滑块上的力，从而使滑块更迅速地停在目标位置。这里我们将采用控制理论中最常用的控制方法：PID控制。采用PID控制器的系统框图如图3-28所示。
图3-28 PID控制系统
图3-28中，Setpoint为控制的目标状态，在本例中为滑块的目标位移，控制目标系统Process（即质量-弹簧-阻尼系统）的输出Output为系统的输出，在本例中为滑块的实际位移，二者的误差Error作为PID系统的输入。PID控制器由三个独立的部分构成：
●P：比例项，输出和误差成正比。
●I：积分项，输出和误差的积分成正比。
●D：微分项，输出和误差的微分成正比。
三项之和作为PID控制器的输出，在本例中控制器的输出为作用在滑块之上的外力F。系统框图中使用积分和微分符号表示积分项和微分项，在实际的控制系统中，我们以一定的时间间隔计算控制器的输出，这时积分项可以用累加变量self.status表示，而微分项可以用当前的误差以及前次的误差self.last_error之间的差进行计算。
下面是PID控制器的程序，我们使用类把kp、ki、kd、dt等参数、累加变量status以及上一次的误差last_error等包装起来。
class PID(object):
def __init__(self, kp, ki, kd, dt):
self.kp, self.ki, self.kd, self.dt = kp, ki, kd, dt
self.last_error = None
self.status = 0.0
def update(self, error):
p = self.kp *
error
i = self.ki *
self.status
if self.last_error is None:
d = 0.0
else:
d = self.kd *
(error - self.last_error) / self.dt
self.status += error *
self.dt
self.last_error = error
return p + i + d
下面的程序使用PID控制器对系统进行控制，让滑块更快地停在位移1.0处，为了后续调用优化工具自动搜索合适的PID参数，这里将整个系统模拟用函数pid_control_system( )封装起来，函数的参数为PID控制器的三个参数。
程序的基本构造与前面的无控制的程序相同，只是增加了PID控制器方面的运算：❶计算目标位置1.0与当前位置之间的误差，❷使用该误差更新PID控制器，获得控制器的输出F，❸更新目标系统中的控制力。
由程序的输出可知，由于PID控制器的控制，系统在两秒之内就已经停在了位移1.0处，如图3-29所示。
图3-29 使用PID控制器让滑块停在位移1.0处
def pid_control_system(kp, ki, kd, dt, target=1.0):
system = MassSpringDamper(m=m, k=k, b=b, F=0.0)
pid = PID(kp, ki, kd, dt)
init_status = 0.0, 0.0
r = ode(system.f)
r.set_integrator('vode', method='bdf')
r.set_initial_value(init_status, 0)
t = [0]
result = [init_status]
F_arr = [0]
while r.successful( ) and r.t + dt < 2.0:
r.integrate(r.t + dt)
err = target - r.y[0] ❶
F = pid.update(err) ❷
system.F = F ❸
t.append(r.t)
result.append(r.y)
F_arr.append(F)
result = np.array(result)
t = np.array(t)
F_arr = np.array(F_arr)
return t, F_arr, result
t, F_arr, result = pid_control_system(50.0, 100.0, 10.0, 0.001)
print u"控制力的终值:", F_arr[-1]
控制力的终值: 19.9434046839
通过调节PID控制器的三个参数可以获得最佳的控制效果，这里我们使用前面介绍过的optimize库中的函数自动寻找最优的PID参数。为了使用最优化函数，需要编写一个对控制结果进行评价的函数。由于我们的目标是让滑块尽快地停在位移1.0处，因此可以用前两秒钟滑块位移与目标位移差的绝对值之和作为控制结果的评价，该值越小表示控制得越好。为了让最优化运行得快一些，这里将控制器的时间间隔改为0.01秒。
%%time
from scipy import optimize
def eval_func(k):
kp, ki, kd = k
t, F_arr, result = pid_control_system(kp, ki, kd, 0.01)
return np.sum(np.abs(result[:, 0] - 1.0))
kwargs = {"method":"L-BFGS-B",
"bounds":[(10, 200), (10, 100), (1, 100)],
"options":{"approx_grad":True}}
opt_k = optimize.basinhopping(eval_func, (10, 10, 10),
niter=10,
minimizer_kwargs=kwargs)
print opt_k.x
[ 199.81255771 100. 15.20382074]
Wall time: 1min 15s
下面使用优化器的输出作为PID的参数对系统进行模拟，可以看到控制开始0.5秒之后滑块已经基本上稳定在了位移1.0处，如图3-30所示。
图3-30 优化PID的参数以降低控制响应时间
kp, ki, kd = opt_k.x
t, F_arr, result = pid_control_system(kp, ki, kd, 0.01)
idx = np.argmin(np.abs(t - 0.5))
x, u = result[idx]
print "t={}, x={:g}, u={:g}".format(t[idx], x, u)
t=0.5, x=0.979592, u=0.00828481