使用 ModelingToolkit.jl 消除守恒量

Question

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):

然后我找到了关于 DAE 索引缩减的教程，并认为这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目前是否有一个功能可以进行转换，或者是否需要采取不同的方法？

Answer 1

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

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)仍然是一样的。在这种特定情况下，也许有一天可以证明，给定微分方程组，守恒定律一定成立。但即使可以，那么格式将是您只给出微分方程，然后让它通过替换证明的守恒定律来删除方程：所以您仍然只给出两个状态系统的两个方程。