求解隐式 ODE（微分代数方程 DAE）

Posted 2023-03-12

【中文标题】求解隐式 ODE（微分代数方程 DAE）

【英文标题】：Solve an implicit ODE (differential algebraic equation DAE)

【发布时间】：2014-06-28 00:27:27

【问题描述】：

我正在尝试使用 scipy 中的 odeint 解决二阶 ODE。我遇到的问题是函数隐式耦合到二阶项，如简化的 sn-p 所示（请忽略示例的假装物理）：

import numpy as np
from scipy.integrate import odeint

def integral(y,t,F_l,mass):
    dydt = np.zeros_like(y)
    x, v = y
    F_r =  (((1-a)/3)**2 + (2*(1+a)/3)**2) * v # 'a' implicit 
    a  = (F_l - F_r)/mass

    dydt = [v, a]
return dydt


y0 = [0,5]
time = np.linspace(0.,10.,21)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon,mass))

在这种情况下，我意识到可以对隐式变量进行代数求解，但是在我的实际场景中，F_r 和 a 的评估之间存在很多逻辑，并且代数操作失败。

我相信 DAE 可以使用 MATLAB 的 ode15i 函数来解决，但我会尽可能避免这种情况。

我的问题是 - 有没有办法解决 python 中的隐式 ODE 函数 (DAE)（最好是 scipy）？有没有更好的方法来解决上述问题？

作为最后的手段，从前一个时间步传递a 可能是可以接受的。如何在每个时间步之后将dydt[1] 传递回函数？

【问题讨论】：

你也可以看看DAE Tools。

【参考方案1】：

相当旧，但值得更新，因此它可能对任何偶然发现这个问题的人有用。目前在 python 中可用的包很少，可以解决隐式 ODE。 GEKKO (https://github.com/BYU-PRISM/GEKKO) 是软件包之一，专门针对混合整数、非线性优化问题进行动态优化，但也可以用作通用 DAE 求解器。

上述“假装物理”问题可以在 GEKKO 中解决如下。

m= GEKKO()
m.time = np.linspace(0,100,101)
F_l = m.Param(value=1000)
mass = m.Param(value =1000)
m.options.IMODE=4
m.options.NODES=3
F_r = m.Var(value=0)
x = m.Var(value=0)
v = m.Var(value=0,lb=0)
a = m.Var(value=5,lb=0)
m.Equation(x.dt() == v)
m.Equation(v.dt() == a)
m.Equation (F_r ==  (((1-a)/3)**2 + (2*(1+a)/3)**2 * v)) 
m.Equation (a == (1000 - F_l)/mass)
m.solve(disp=False)
plt.plot(x)

【讨论】：

您好，请问您是否知道是否可以用非常数系数解决这种二阶隐式问题？例如，采用等式 $f(x)y''(x) + g(x) y'(x)*y(x)$。快速浏览一下文档并没有给我一个明确的答案，我也没有找到涵盖此案例的示例。

【参考方案2】：

如果代数运算失败，您可以寻求约束的数值解，例如在每个时间步运行 fsolve：

import sys
from numpy import linspace
from scipy.integrate import odeint
from scipy.optimize import fsolve

y0 = [0, 5]
time = linspace(0., 10., 1000)
F_lon = 10.
mass = 1000.

def F_r(a, v):
    return (((1 - a) / 3) ** 2 + (2 * (1 + a) / 3) ** 2) * v

def constraint(a, v):
    return (F_lon - F_r(a, v)) / mass - a

def integral(y, _):
    v = y[1]
    a, _, ier, mesg = fsolve(constraint, 0, args=[v, ], full_output=True)
    if ier != 1:
        print "I coudn't solve the algebraic constraint, error:\n\n", mesg
        sys.stdout.flush()
    return [v, a]

dydt = odeint(integral, y0, time)

显然，这会减慢您的时间整合速度。始终检查 fsolve 是否找到了一个好的解决方案，并刷新输出，以便您可以在它发生时意识到它并停止模拟。

关于如何在前一个时间步“缓存”变量的值，您可以利用默认参数仅在函数定义时计算的事实，

from numpy import linspace
from scipy.integrate import odeint

#you can choose a better guess using fsolve instead of 0
def integral(y, _, F_l, M, cache=[0]):
    v, preva = y[1], cache[0]
    #use value for 'a' from the previous timestep
    F_r = (((1 - preva) / 3) ** 2 + (2 * (1 + preva) / 3) ** 2) * v 
    #calculate the new value
    a = (F_l - F_r) / M
    cache[0] = a
    return [v, a]

y0 = [0, 5]
time = linspace(0., 10., 1000)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon, mass))

请注意，为了使技巧起作用，cache 参数必须是可变的，这就是我使用列表的原因。如果您不熟悉默认参数的工作原理，请参阅 this 链接。

请注意，这两个代码不会产生相同的结果，并且您应该非常小心地使用前一个时间步的值，以确保数值稳定性和精度。第二个显然要快得多。

【讨论】：

嗨，谢谢 flebol，这里对这两种解决方案的基本分析给出了：%timeit odeint(integral1, y0, time)100 loops, best of 3: 9.03 ms per loop%timeit odeint(integral2, y0, time, args=(F_lon, mass))1000 loops, best of 3: 972 µs per loopintegral1是第一个示例。在这些条件下，两个示例之间的差异相当小（对于 1000 个时间步长，大约为 10^-7），但是，在第二个示例中更改某些值（例如，'F_lon=1000; mass=100'）方法失败。出于这个原因，我可以忍受时间惩罚。谢谢。

@Shaun_M，9.03ms 和 972micro s 之间的差异是 10 倍，第一个解决方案慢了 10 倍

抱歉不清楚。我的意思是解决方案之间的数值差异。我想说的是，在这种情况下，使用缓存积分的精度是可以接受的，但稳定性不是。谢谢

@Shaun_M 感谢您提供解决 DAE 的有趣方法。您的解决方案能否适应具有奇异“质量矩阵”的 ode 问题，例如 M*du/dt=f(u) ？