ModelingToolkit: Index reduction and linear algebraic loops

Specific Domains Modelling & Simulations
1

Hey all,

I’m trying out ModelingToolkit.jl and want to investigate what is currently possible and what not and compare it to the OpenModelica Compiler Backend.

Let’s say we have the following DAE:

der(x) = y + a
y = 3*x
der(y) = a*x
b = a/5

which should be simple but interesting example because it has:

  • Differential index 2 (Maximum number of differentiations of all equations needed to get a solvable ODE).
  • A linear loop.
  • Algebraic variables.

To get an executable simulation from this OpenModelica will do:

  • Matching and sorting
  • Index reduction with dummy-states

See paper The OpenModelica Integrated Environment for Modeling, Simulation, and Model-Based Development chapter 4 on how the OpenModelica Backend is working.

I’m trying to solve the example from above with ModelingToolkit.jl (v6.5.2) but wasn’t successful yet, so I guess I’m missing something. I’m following the documentation and the pendulum example: Automated Index Reduction of DAEs · ModelingToolkit.jl

using ModelingToolkit
using OrdinaryDiffEq

@parameters t
@variables x(t) y(t) a(t) b(t)
der = Differential(t)

eqs = [der(x) ~ y + a,
       0 ~ 3*x - y,
       der(y) ~ a*x,
       0 ~ a/5 - b]

@named sys = ODESystem(eqs,t)

sys = dae_index_lowering(sys)
sys = structural_simplify(sys)

u0 = [y => 1.0,
      x => 3.0,
      a => 0.0]
p = []

tspan = (0.0,1.0)
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Rodas4())

but get

┌ Warning: Instability detected. Aborting
└ @ SciMLBase C:\Users\...\.julia\packages\SciMLBase\n3U0M\src\integrator_interface.jl:351

during solving of the problem and no solution.

I think the problem is in the ODESystem I’m declaring right at the beginning:

@named sys = ODESystem(eqs,t)

In fact the system I want to describe is not an ODE but an DAE.
Also usually I would not need to give start values for all x,y,a and b. Only one (x or y) is a state, so only one of these needs a start value.

I’m not sure if that is the problem or that we (should) have a linear equation system where we need a linear solver to solve several equations at once.

1 Like
2

Try specifying different initial conditions:

u0 = [x => 1.0,
      y => 3.0,
      a => 0.0]
p = []

tspan = (0.0,1.0)
prob = ODEProblem(sys_simp, u0, tspan, p)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 3-element Vector{Float64}:
 1.0
 3.0
 0.0

sol = solve(prob, Rodas4())

retcode: Success
Interpolation: specialized 3rd order "free" stiffness-aware interpolation
t: 17-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.03479403880815428
 0.06638720509605027
 0.10919245871227796
 0.16354007662252074
 0.23195113377060073
 0.31589925312034817
 0.41802992132379013
 0.5409179851394389
 0.6879618323552728
 0.8629971165185748
 1.0
u: 17-element Vector{Vector{Float64}}:
 [1.0, 3.0, -4.5]
 [0.9999985000028127, 2.999995500008438, -4.4999898750265785]
 [0.9999835003403045, 2.9999505010209133, -4.499888628215853]
 [0.9998335346445447, 2.999500603933634, -4.498876452366555]
 [0.9983369632571392, 2.995010889771416, -4.488783828438379]
 [0.9836726314301416, 2.9510178942904224, -4.390682692451393]
 [0.9509782787859907, 2.8529348363579694, -4.177019806867765]
 [0.9112823663520662, 2.7338470990561956, -3.92659142890958]
 [0.8634629624516273, 2.5903888873548793, -3.6372714698271382]
 [0.8107083942206305, 2.432125182661888, -3.33275521142727]
 [0.7541425065360393, 2.2624275196081145, -3.02213238193428]
 [0.6960505334773193, 2.0881516004319547, -2.719004903299619]
 [0.6377736294194034, 1.913320888258207, -2.4298929030455176]
 [0.580686254432733, 1.7420587632981954, -2.1601871937913346]
 [0.5256735953773559, 1.5770207861320642, -1.9120585066743794]
 [0.4734288967326366, 1.420286690197906, -1.68641802306022]
 [0.4397881499773854, 1.3193644499321526, -1.5460020112815398]

Looks like the problem was inconsistent initial conditions (x and y switched?). For DAE systems consistent initial conditions are important.

2 Likes
3

Yeah, it looks like a problem with the consistent initialization algorithm BrownFullBasicInit in OrdinaryDiffEq. If you do @show solve(prob, Rodas4())[1] you see:

(solve(prob, Rodas4()))[1] = [3.0, 1.0, -Inf]

Oops, for some reason the initializer diverged, and I can see it diverges in one step too which is odd. If you can open an issue on OrdinaryDiffEq.jl we can follow up there. @SCher 's answer is correct but we might as well fix this too.

1 Like
4

Indeed I used an incorrect initialization. With @SCher u0 the simulation works. Thanks!

I’ll open an issue on SciML/OrdinaryDiffEq.jl about this. I would like a warning that will tell the user that the initial system could not be solved because it is inconsistent due to 3 != 0.3333 = start(x).

What values of the initialization are needed and which ones are optional?

For this example I know, that x is a state and y is a dummy-state (or the other way around) and that a and b are algebraic variables wich are inside a linear loop.
So I would expect that one start value for x or y would be enough. When a and b are inside a non-linear loop I would need to set start values for them as well, so the non-linear iterative solver (e.g. Newton) can converge to the solution.

How can I get this information from MTK?

julia> states(sys)
4-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 x(t)
 y(t)
 a(t)
 b(t)

Calling states(sys) will give me all state varaibles x and algebraic variables y. But the information I’m looking for is what are the unknowns z=(der(x),y) of my system, so I know which variables I’ll need to give a meaningful start value.

Of course this information is only available after matching, sorting and index reduction have figured out what variables are the states and which variables are solved inside a non-linear strong component (non-linear loop).

5

You can do dae_index_lowering() and structural_simplify() on your system, which will create a minimal and most efficient system for numerical simulation, partitioning the original system into dynamic and algebraic equations.

sys_low = dae_index_lowering(sys)

Differential(t)(x(t)) ~ a(t) + y(t)
Differential(t)(y(t)) ~ a(t)*x(t)
0 ~ (1//5)*a(t) - b(t)
0 ~ 3a(t) + 3y(t) - a(t)*x(t)

sys_simp = structural_simplify(sys_low)

Differential(t)(x(t)) ~ a(t) + y(t)
Differential(t)(y(t)) ~ a(t)*x(t)
0 ~ 3a(t) + 3y(t) - a(t)*x(t)

states(sys_simp)

3-element Vector{Any}:
 x(t)
 y(t)
 a(t)

For large nonlinear algebraic systems of equations, I guess under the hood structural_simplify performs optimal tearing of equations (after BLT transformation) to figure out which variables to use for most efficient numerical formulation.
Coming from OpenModelica and Dymola background, I think selection of the right variables to initialize is more of an art than exact science though arguably Dymola gives better info to users about tearing variables.
ModelingToolkit tearing() function will probably give more info about what’s going on in the background, however I haven’t see much documentation on this. Hopefully will come soon.
@ChrisRackauckas and @YingboMa should be able to provide more info.

6

I see how transforming the DAE f(x(t),dx(t),y(t),u(t),p,t) into two sets of equations

dx(t) = h(x(t),u(t),p,t)
y(t) = k(x(t),u(t),p,t)

is done by dae_index_lowering and structural_simplify.
For notation: x are states, dx are state derivatives, y are algebraic variables, u are input variables, p are paramters/constants, t is time.
I think in the docs of MTK states are interpreted as (x,y) (so with algebraic variables). An that is what states(sys) is returning.

But currently dae_index_lowering() won’t find that there are dummy-states in my example above and therefore only variable y would need a start value, but expects me to give start values for x, y and b.
To get the correct start values I would need to solve the initial system myself or specify in some way which start value I care about and which start values should be find automatically.

But that can be solved. Collection the unknowns (dx,y) is already done during Pantelides, so it should be possible to return x from the states(sys) or some similar function.
State selection · Issue #988 · SciML/ModelingToolkit.jl · GitHub is talking about adding the dumy-derivative method, so that should work.

Coming from OpenModelica and Dymola background, I think selection of the right variables to initialize is more of an art than exact science though arguably Dymola gives better info to users about tearing variables.

I would say getting the variables that need a start value is possible to get from a tool. These are the states (without dummy-states) and non-linear iteration (tearing) variables.
Getting the value right for these variables on the other hand is no easy task. There are a few ideas around, see e.g. this pre-print: https://www.researchgate.net/publication/337672027_On_the_choice_of_initial_guesses_for_the_Newton-Raphson_algorithm

7

Yeah, this looks like it needs dummy derivative stuff.