Python求解含嵌入导数的电池热模型ODE方程咨询

你的核心问题是：给定的方程是微分代数方程（DAE），而非普通ODE——它包含状态变量导数之间的代数约束，且你尝试的两种模型定义都未正确构建封闭的求解系统。以下是可行的解决方法：

问题分析 #

你给出的方程：

dTs/dt = Ts*a + Ta*b + (dTa/dt)*c + d

要同时求解Ts和Ta，系统需要两个独立方程。你的两种尝试均存在明显错误：

• 第一种尝试：状态向量包含Ts和Ta两个变量，但仅返回dTs/dt一个导数，长度不匹配，求解器无法处理。
• 第二种尝试：错误地将导数v1=dTs/dt、v2=dTa/dt作为状态变量，但未定义这些导数的变化率（即二阶导数），系统不封闭。

可行解决方案 #

方案1：补充缺失方程，构建封闭ODE系统 #

如果Ta是未知状态变量（比如冷却介质温度），你需要补充Ta的热平衡方程。假设Ta的微分方程为：

dTa/dt = m*Ts + n*Ta + p

将其代入原方程，转化为标准一阶ODE系统：

dTs/dt = (a + c*m)*Ts + (b + c*n)*Ta + (c*p + d)
dTa/dt = m*Ts + n*Ta + p

代码示例（适配离散时间步的系数）：

import numpy as np
from scipy.integrate import odeint

# 示例时间数组与系数（替换为你的实际数据）
t = np.linspace(0, 10, 100)
a = np.random.randn(100)
b = np.random.randn(100)
c = np.random.randn(100)
d = np.random.randn(100)
m = np.random.randn(100)  # Ta方程的系数
n = np.random.randn(100)
p = np.random.randn(100)

def Tmodel(z, t_idx):
    Ts, Ta = z
    # 获取当前时间步的系数
    ai, bi, ci, di = a[t_idx], b[t_idx], c[t_idx], d[t_idx]
    mi, ni, pi = m[t_idx], n[t_idx], p[t_idx]
    
    dTsdt = (ai + ci*mi)*Ts + (bi + ci*ni)*Ta + (ci*pi + di)
    dTadt = mi*Ts + ni*Ta + pi
    return [dTsdt, dTadt]

# 初始条件与求解
z0 = [0.0, 0.0]
Ts = np.zeros_like(t)
Ta = np.zeros_like(t)
Ts[0], Ta[0] = z0

for i in range(len(t)-1):
    tspan = [t[i], t[i+1]]
    z = odeint(Tmodel, z0, tspan, args=(i,))
    Ts[i+1] = z[-1][0]
    Ta[i+1] = z[-1][1]
    z0 = z[-1]

方案2：将DAE转化为隐式ODE（当c≠0时） #

若只有给定的一个方程，可将其整理为隐式ODE系统，使用支持隐式求解的solve_ivp（BDF/Radau方法）处理：

import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d

# 离散时间与系数（替换为你的数据）
t_data = np.linspace(0, 10, 100)
a_data = np.random.randn(100)
b_data = np.random.randn(100)
c_data = np.random.randn(100) + 2  # 确保c≠0，避免除以0
d_data = np.random.randn(100)

# 构建连续插值函数，用于任意时间点的系数查询
a = interp1d(t_data, a_data, kind='linear', fill_value="extrapolate")
b = interp1d(t_data, b_data, kind='linear', fill_value="extrapolate")
c = interp1d(t_data, c_data, kind='linear', fill_value="extrapolate")
d_val = interp1d(t_data, d_data, kind='linear', fill_value="extrapolate")

def implicit_ode(t, z, dzdt):
    Ts, Ta = z
    dTsdt, dTadt = dzdt
    # 原方程约束
    eq1 = dTsdt - a(t)*Ts - b(t)*Ta - c(t)*dTadt - d_val(t)
    # 从原方程解出dTa/dt的约束
    eq2 = dTadt - (dTsdt - a(t)*Ts - b(t)*Ta - d_val(t))/c(t)
    return [eq1, eq2]

# 初始条件与求解
z0 = [0.0, 0.0]
dzdt0 = [0.0, 0.0]  # 初始导数猜测值

sol = solve_ivp(implicit_ode, t_span=[t_data[0], t_data[-1]], 
                y0=z0, y0_dot=dzdt0, method='BDF', t_eval=t_data)

Ts_sol = sol.y[0]
Ta_sol = sol.y[1]

方案3：使用专门的DAE求解库 #

对于复杂DAE系统，可使用casadi或pyomo.dae等专业库，示例（casadi）：

import casadi as ca
import numpy as np

# 时间网格与系数
t = np.linspace(0, 10, 100)
n = len(t)
a_data = np.random.randn(n)
b_data = np.random.randn(n)
c_data = np.random.randn(n)+2
d_data = np.random.randn(n)

# 定义变量与插值系数
Ts = ca.MX.sym('Ts')
Ta = ca.MX.sym('Ta')
dTsdt = ca.MX.sym('dTsdt')
dTadt = ca.MX.sym('dTadt')

a = ca.interpolant('a', 'linear', [t], a_data)
b = ca.interpolant('b', 'linear', [t], b_data)
c = ca.interpolant('c', 'linear', [t], c_data)
d_val = ca.interpolant('d', 'linear', [t], d_data)

# 定义DAE残差
res = ca.vertcat(
    dTsdt - a(t)*Ts - b(t)*Ta - c(t)*dTadt - d_val(t),
    dTsdt - ca.diff(Ts, t),
    dTadt - ca.diff(Ta, t)
)

# 创建求解器并运行
dae = {'x': ca.vertcat(Ts, Ta), 'p': t, 'ode': res}
solver = ca.integrator('solver', 'idas', dae, {'tf': t[-1]})
sol = solver(x0=[0.0, 0.0])

Ts_sol = np.array(sol['xf'])[0]
Ta_sol = np.array(sol['xf'])[1]

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