Python: 一个带积分项的常微分方程怎么解

Question

我不知道之前是否有人在 SO 中问过这个问题，我会继续 post 在这里，我试图用 PID 控制器求解一个简单的系统，我的微分方程组是下面给出。我基本上是在尝试编写非常基本的 PID 算法。我的控件 u 的结构取决于误差项的导数和积分。我对导数项没有任何问题，它是在我的代码中产生问题的积分项。当我在开始时分配 s=0 时，问题就出现了并按照下面代码中的描述在我的函数中使用它。有没有办法绕过它？我尝试分配 s 并告诉为全局变量，但它没有解决我的问题。简而言之，我正在做的是 - 我每次都添加状态 x1 并乘以 dt（由 t-told 表示）。

请帮我解决这个问题，PFA 下面附上我的代码。

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
plt.style.use('bmh')

t0=0
y0=[0.1,0.2]
kp,kd,ki=2,0.5,0.8
s,told=0,0
def pid(t,Y):
    x1,x2=Y[0],Y[1]
    e=x1-1
    de=x2
    s=(x1+s)
    integral=s*(t-told)
    told=t
    #ie=
    u=kp*e+kd*de+ki*integral
    x1dot=x2
    x2dot=u-5*x1-2*x2
    return[x1dot,x2dot]

solver=ode(pid).set_integrator('dopri5',rtol=1e-6,method='bdf',nsteps=1e5,max_step=1e-3)
solver.set_initial_value(y0,t0)
t1=10
dt=5e-3
sol = [ [yy] for yy in y0 ]
t=[t0]
while solver.successful() and solver.t<t1:
    solver.integrate(solver.t+dt)
    for k in range(2): sol[k].append(solver.y[k]);
    t.append(solver.t)
    print(len(sol[0]))
    print(len(t))
x1=np.array(sol[0])
x2=np.array(sol[1])
e=x1-1
de=x2
u=kp*e+kd*de
for k in range(2):
    if k==0:
        plt.subplot(2,1,k+1)
        plt.plot(t,sol[k],label='x1')
        plt.plot(t,sol[k+1],label='x2')
        plt.legend(loc='lower right')
    else:
        plt.subplot(2,1,k+1)
        plt.plot(t,u)
plt.show()

Answer 1

首先你需要在pid函数中包含"s"变量。

' def pid(s, t, Y): ... '

Answer 2

我现在能看到的最简单的解决方案是创建一个 class，其中 s 和 told 作为此 class 的属性：

class PIDSolver:
    def __init__(self)
        self.t0=0
        self.y0=[0.1,0.2]
        self.kp,self.kd,self.ki=2,0.5,0.8
        self.s,self.told=0,0

    def pid(t,Y):
        x1,x2=Y[0],Y[1]
        e=x1-1
        de=x2
        self.s=(x1+self.s)
        integral=self.s*(t-self.told)
        self.told=t
        #ie=
        u=self.kp*e+self.kd*de+self.ki*integral
        x1dot=x2
        x2dot=u-5*x1-2*x2
        return[x1dot,x2dot]

对于你问题的第一部分。在您的解决方案的下一部分使用 pidsolver = PIDSolver()。

Answer 3

我自己解决了这个问题，方法是使用 set_f_params() 方法并在 itz 参数中传递一个列表。此外，我在 pid() 中传递了第三个参数，即 pid(t,Y,arg)。最后我分配了 s,told=arg[0],arg[1]。

Answer 4

您对求解器及其访问的时间步长做出了不合理的假设。通过你对积分的破解，即使它在数学上是合理的（它应该看起来像 integral = integral + e*(t-told)，它给出了阶数为 1 的积分方法），你可以减少任何积分方法的阶数，可能降到 1，如果你是幸运只订购了 2.

实现这个系统的数学上正确的方法是为e的积分引入第三个变量x3，即x3的导数是e.正确的 order 1 系统必须是维度 3 可以理解为（消除 e）您的系统有 3 个 differentiation/integration 操作。这样你的系统就变成了

def pid(t,Y):
    x1, x2, x3 =Y
    e=x1-1
    x1dot = x2
    edot = x1dot
    x3dot = e
    u=kp*e+kd*edot+ki*x3
    x2dot=u-5*x1-2*x2
    return[x1dot, x2dot, x3dot]

请注意，不需要全局动态变量，只有常量（也可以作为参数传递，任何看起来更有效或可读的东西）。

现在您还需要 x3 的初始值，从系统中看不到集成变量必须是什么，您的代码似乎建议 0.

Python: 一个带积分项的常微分方程怎么解

Python: How to solve an ordinary differential equation with integral term in it

python

numerical-integration

ode

differential-equations