使用ModelingToolkit.jl消去一个守恒量

Use ModelingToolkit.jl to eliminate a conserved quantity

ModelingToolkit.jl 是一个非常棒的包,我经常对它抱有太多期望。例如,我经常发现自己的模型归结为以下内容:

@variables t x(t) y(t)
@parameters a b C
d = Differential(t)
eqs = [
    d(x) ~ a * y - b * x,
    d(y) ~ b * x - a * y,
    0 ~ x + y - C
]

@named sys = ODESystem(eqs)

现在, 知道我可以通过代入 0 ~ x + y - C 将其归结为一个方程式。但实际上我的系统要大得多,不那么琐碎并且是通过编程生成的,所以我想 ModelingToolkit.jl 为我做这件事。

我已经尝试使用structural_simplify,但是额外的等式妨碍了:

julia> structural_simplify(sys)
ERROR: ExtraEquationsSystemException: The system is unbalanced. There are 2 highest order derivative variables and 3 equations.
More equations than variables, here are the potential extra equation(s):

然后我找到了 tutorial on DAE index reduction,并认为 dae_index_lowering 可能适合我:

julia> dae_index_lowering(sys)
ERROR: maxiters=8000 reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiters`. Try to increase the maxiters by `pantelides(sys::ODESystem; maxiters=1_000_000)` if your system has an incredibly high index and it is truly extremely large.

所以问题是 ModelingToolkit.jl 目前是否有一个功能可以进行转换,或者是否需要不同的方法?

问题在于系统是不平衡的,即方程的数量多于状态的数量。一般来说,不可能证明这种超定系统是明确定义的。因此,要解决它,您必须删除其中一个方程式。如果你知道守恒定律必须成立,那么你可以删除第二个微分方程:

using ModelingToolkit
@variables t x(t) y(t)
@parameters a b C
d = Differential(t)
eqs = [
    d(x) ~ a * y - b * x,
    0 ~ x + y - C
]

@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)

这将简化为一个方程式。问题是,一般来说,它不能证明如果删除那个微分方程,y(t) 仍然是一样的。在这种特定情况下,也许有一天它可以证明在给定微分方程系统的情况下守恒定律必须出现。但即使可以,那么格式将是你只给出微分方程,然后通过代入已证明的守恒定律让它去掉方程:所以你仍然只给出两个状态系统的两个方程。