如何在scipy.solve_ivp(lsoda法)中获取误差的时间导数

针对你在物理系统仿真中PID控制环的误差导数获取问题，以下提供两种符合需求的解决方案，均保留lsoda求解器的使用：

一、解析法计算误差导数（精准无近似） #

误差var_error = pos - pos_control，其中pos_control是vel的函数，根据链式法则可直接推导误差的时间导数，无需额外引入状态变量：

• d(pos)/dt就是当前状态中的vel
• d(pos_control)/dt = d(pos_control)/d(vel) * d(vel)/dt，其中d(vel)/dt是加速度acc，d(pos_control)/d(vel)可对pos_control = exp(0.25*vel)求导得到0.25*exp(0.25*vel)

修改后的代码如下：

import numpy as np
from scipy.integrate import solve_ivp

def rhs_func(t, rhs_integrated):
    pos = rhs_integrated[0]
    vel = rhs_integrated[1]
    var_error_integrated = rhs_integrated[2]

    # 控制增益
    i_gain = 0.1
    p_gain = 0.1
    d_gain = 0.1
    
    pos_control = np.exp(0.25*vel)
    var_error = pos - pos_control
    
    # 解析计算误差导数
    d_pos_control_dvel = 0.25 * np.exp(0.25*vel)
    acc = 10 - 0.5*vel  # 先计算加速度，因为误差导数需要用到它
    var_error_derivative = vel - d_pos_control_dvel * acc
    
    # 更新速度（加入D项）
    new_vel = vel - p_gain*var_error + i_gain*var_error_integrated + d_gain*var_error_derivative
    
    return np.hstack((new_vel, acc, var_error))

# 初始状态：pos=10, vel=0, 误差积分=0
sol = solve_ivp(rhs_func, [0, 10], [10, 0, 0], t_eval=np.arange(1,11))
print(sol)

二、数值近似法（模拟实际控制系统的固定dt差分） #

要模拟实际控制中导数=(最新值-旧值)/dt的近似方式，需在求解器内部记录历史采样点的误差和时间。由于rhs_func是无状态的，可通过闭包保存历史变量，同时设置固定采样间隔dt，仅当求解器步长跨越采样点时更新导数：

import numpy as np
from scipy.integrate import solve_ivp

def create_rhs_func(dt):
    # 闭包存储历史状态：上一次采样的时间和误差
    prev_t = None
    prev_error = None
    var_error_derivative = 0.0  # 初始导数设为0

    def rhs_func(t, rhs_integrated):
        nonlocal prev_t, prev_error, var_error_derivative
        pos = rhs_integrated[0]
        vel = rhs_integrated[1]
        var_error_integrated = rhs_integrated[2]

        # 控制增益
        i_gain = 0.1
        p_gain = 0.1
        d_gain = 0.1
        
        pos_control = np.exp(0.25*vel)
        var_error = pos - pos_control
        
        # 当时间到达采样间隔时更新导数
        if prev_t is not None and t - prev_t >= dt - 1e-6:  # 加小容差避免浮点误差
            var_error_derivative = (var_error - prev_error) / dt
            prev_t = t
            prev_error = var_error
        elif prev_t is None:
            # 初始化历史值
            prev_t = t
            prev_error = var_error
        
        # 更新速度
        new_vel = vel - p_gain*var_error + i_gain*var_error_integrated + d_gain*var_error_derivative
        acc = 10 - 0.5*new_vel
        
        return np.hstack((new_vel, acc, var_error))
    
    return rhs_func

# 设置实际控制的采样间隔dt=0.1
rhs = create_rhs_func(dt=0.1)
sol = solve_ivp(rhs, [0, 10], [10, 0, 0], t_eval=np.arange(1,11))
print(sol)

注意事项 #

• 解析法完全依赖系统模型的导数推导，误差最小，但需要确保pos_control的导数可解析求解
• 数值近似法更贴近实际控制器的实现，但需注意：lsoda是自适应步长求解器，可能会跳过部分采样点，可通过设置max_step=dt强制求解器步长不超过采样间隔，让差分更准确
• 若要评估近似方法的误差，可将两种方法的仿真结果对比，分析输出差异

内容的提问来源于stack exchange，提问作者ZiePeregrine