Use ModelingToolkit.jl to eliminate a conserved quantity

user_8791491 0
Viewed 156

ModelingToolkit.jl is such a great package that I frequently expect too much of it. For example, I often find myself with a model which boils down to the following:

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

Now, I know that I could get this down to one equation by substitution of the 0 ~ x + y - C. But in reality my systems are much larger, less trivial and programmatically generated, so I would like ModelingToolkit.jl to do it for me.

I have tried using structural_simplify, but the extra equation gets in the way:

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

Then I found the tutorial on DAE index reduction, and thought that dae_index_lowering might work for me:

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.

So the question is whether ModelingToolkit.jl currently has a feature which will do the transformation, or if a different approach is necessary?

julia modeling
1 Answers

The problem is that the system is unbalanced, i.e. there are more equations than there are states. In general it is impossible to prove that an overdetermined system of this sort is well-defined. Thus to solve it, you have to delete one of the equations. If you know the conservation law must hold true, then you can delete the second differential equation:

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)

And that will simplify down to a single equation. The problem is that in general it cannot prove that if it does delete that differential equation, that y(t) is still going to be the same. In this specific case, maybe it could one day prove that the conservation law must occur given the differential equation system. But even if it could, then the format would be for you to only give the differential equation and then let it remove equations by substituting proved conservation laws: so you would still only give two equations for the two state system.

user_15442030
Related