有没有更好的方法来使用 Scipy 的 solve_ivp() 在 Python 中模拟 PID 控制? [英] Is there a better way to simulate PID control in Python with Scipy's solve_ivp()?

问题描述

我正在解决家庭作业问题.我正在尝试使用 Scipy 的 integrate.solve_ivp() 函数在 Python 中模拟 PID 控制.

我的方法是在函数的右侧运行 PID 代码，使用全局变量并在每个时间步的末尾将它们附加到全局矩阵，如下所示:

solution =integrate.solve_ivp(rhs, tspan, init, t_eval=teval)

这是我的代码:

def rhs(dt, init):全局 old_time、omega0dot、rhs_t、omega0dotmat时间步长 = dt - old_timeold_time = dt# 解压初始x = 初始化 [0]y = 初始化 [1]z = 初始化[2]xdot = init[3]ydot = 初始化 [4]zdot = init[5]阿尔法=初始化[6]测试版 = 初始化 [7]伽玛 = 初始化 [8]alphadot = init[9]betadot = init[10]伽玛点 = init[11]# 求解方程(xddot, yddot, zddot, alphaddot, betaddot, gammaddot) = 动力学(k_d, k_m, x, y, z, xdot, ydot, zdot, alpha, beta, gamma, alphadot, betadot, gammadot, omega0dot)# 控制系统z_des = 10err_z = z_des - zzPID = (1*err_z) + 悬停omega0dot = zPIDrhs_t.append(dt)omega0dotmat.append(omega0dot)返回 [xdot, ydot, zdot, xddot, yddot, zddot, alphadot, betadot, gammadot, alphaddot, betaddot, gammaddot]

全局变量在这个函数之外被初始化.您可能会注意到，我专门尝试模拟四旋翼飞行器，其中四旋翼飞行器的线性和角运动取决于 omega0dot，它代表转子速度，我试图用 PID 控制.

我的难点在于integrate.solve_ivp() 的时间步长.PID 控制的积分和微分部分都依赖于时间步长，但是 solve_ivp() 函数的时间步长是可变的，有时甚至似乎在时间上倒退，有时没有时间步长(即dt <= 0).

我想知道是否有更好的方法来进行这种 PID 控制，或者我是否在解释 solve_ivp() 中的 dt 术语是错误的.

解决方案

我们来看一个更简单的系统，无处不在的带阻尼弹簧

 y'' + c*y' + k*y = u(t)

其中 u(t) 可以表示电磁体施加的力(这立即提出了通过引入更现实的电压和磁场关系来使系统更复杂的方法).

现在在 PID 控制器中，我们有参考输出的误差 e = yr - y 和

 u(t) = kD*e'(t) + kP*e(t) + kI*integral(e(t))

为了使用 ODE 求解器处理这个问题，我们立即看到需要使用新组件扩展状态 E(t) 其中 E'(t)=e(t).下一个困难是实现不一定可微的表达式的导数.这可以通过使用非标准的一阶实现来完全避免区分这个表达式来实现(标准是使用 [y,y',E] 作为状态).>

本质上，将公式中的所有派生表达式以它们的集成形式收集为

 v(t)=y'(t)+c*y-kD*e(t).

然后回到导数得到一阶系统

 v'(t) = y''(t) + c*y'(t) - kD*e'(t)= kP*e(t) + kI*E(t) - k*y(t)y'(t) = v(t) - c*y(t) + kD*e(t)E'(t) = e(t)

这现在允许将受控系统实现为 ODE 系统，而无需涉及全局内存或类似的技巧

def odePID(t,u,params):c,k,kD,kP,kI = 参数y,v,E = ue = yr(t)-y返回 [ v-c*y+kD*e, kP*e+kI*E-k*y, e ]

您应该能够在更复杂的模型中使用一阶系统的类似转换.

I am working on a homework problem. I'm trying to simulate a PID control in Python with Scipy's integrate.solve_ivp() function.

My method is to run the PID code within the right-hand-side of the function, using global variables and appending them to a global matrix at the end of each timestep, like so:

solution = integrate.solve_ivp(rhs, tspan, init, t_eval=teval)

Here is my code:

def rhs(dt, init):

    global old_time, omega0dot, rhs_t, omega0dotmat
    timestep = dt - old_time
    old_time = dt

    # UNPACK INITIAL
    x = init[0]
    y = init[1]
    z = init[2]
    xdot = init[3]
    ydot = init[4]
    zdot = init[5]
    alpha = init[6]
    beta = init[7]
    gamma = init[8]
    alphadot = init[9]
    betadot = init[10]
    gammadot = init[11]

    # SOLVE EQUATIONS
    (xddot, yddot, zddot, alphaddot, betaddot, gammaddot) = dynamics(k_d, k_m, x, y, z, xdot, ydot, zdot, alpha, beta, gamma, alphadot, betadot, gammadot, omega0dot)

    # CONTROL SYSTEMS
    z_des = 10

    err_z = z_des - z

    zPID = (1*err_z) + hover

    omega0dot = zPID

    rhs_t.append(dt)
    omega0dotmat.append(omega0dot)

    return [xdot, ydot, zdot, xddot, yddot, zddot, alphadot, betadot, gammadot, alphaddot, betaddot, gammaddot]

The global variables are initialized outside this function. You might notice that I am specifically trying to simulate a quadcopter, where the linear and angular motion of the quadrotor are dependent on omega0dot, which represents rotor velocity and which I am trying to control with PID.

My difficulty is with the timestep of integrate.solve_ivp(). Both the integral and derivative part of the PID control rely on the timestep, but the solve_ivp() function has a variable time step and seems to even go backwards in time sometimes, and sometimes makes no timestep (i.e. dt <= 0).

I was wondering if there was a better way to go about this PID control, or if maybe I'm interpreting the dt term in solve_ivp() wrong.

解决方案

Let's look at a more simple system, the ubiquitous spring with dampening

    y'' + c*y' + k*y = u(t)

where u(t) could represent the force exerted by an electro-magnet (which immediately suggests ways to make the system more complicated by introducing a more realistic relation of voltage and magnetic field).

Now in a PID controller we have the error to a reference output e = yr - y and

    u(t) = kD*e'(t) + kP*e(t) + kI*integral(e(t))

To treat this with an ODE solver we immediately see that the state needs to be extended with a new component E(t) where E'(t)=e(t). The next difficulty is to implement the derivative of an expression that is not necessarily differentiable. This is achievable by avoiding to differentiate this expression at all by using a non-standard first-order implementation (where the standard would be to use [y,y',E] as state).

Essentially, collect all derived expressions in the formula in their integrated form as

    v(t)=y'(t)+c*y-kD*e(t). 

Then going back to the derivative one gets the first order system

    v'(t) = y''(t) + c*y'(t) - kD*e'(t) 
          = kP*e(t) + kI*E(t) - k*y(t)

    y'(t) = v(t) - c*y(t) + kD*e(t)
    E'(t) = e(t)

This now allows to implement the controlled system as ODE system without tricks involving a global memory or similar

def odePID(t,u,params):
    c,k,kD,kP,kI = params
    y,v,E = u
    e = yr(t)-y
    return [ v-c*y+kD*e, kP*e+kI*E-k*y, e ]

You should be able to use similar transformations of the first order system in your more complex model.

这篇关于有没有更好的方法来使用 Scipy 的 solve_ivp() 在 Python 中模拟 PID 控制?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！