# ModelingToolkit: Index reduction and linear algebraic loops
**URL:** https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit
**Created:** [October 6, 2021, 3:41pm UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312 "2021-10-06T15:41:50Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@AnHeuermann](https://discourse.julialang.org/u/AnHeuermann)
#### Post date: [October 6, 2021, 3:41pm UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312/1 "2021-10-06T15:41:50Z")
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:
```julia
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](https://www.mic-journal.no/ABS/MIC-2020-4-1.asp/) 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](https://mtk.sciml.ai/stable/mtkitize_tutorials/modelingtoolkitize_index_reduction/#Transforming-Higher-Order-DAEs-to-Index-1-DAEs)
```julia
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
```julia
┌ 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:
```julia
@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.
---
### Author:  [@SCher](https://discourse.julialang.org/u/SCher)
#### Post date: [October 6, 2021, 4:56pm UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312/2 "2021-10-06T16:56:15Z")
> [@AnHeuermann](#):
>
> ```julia
>
> ```
>
> during solving of the problem and no solution.
Try specifying different initial conditions:
```julia
u0 = [x => 1.0,
y => 3.0,
a => 0.0]
p = []
tspan = (0.0,1.0)
prob = ODEProblem(sys_simp, u0, tspan, p)
```
```julia
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())
```julia
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 6, 2021, 5:40pm UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312/3 "2021-10-06T17:40:30Z")
Yeah, it looks like a problem with the consistent initialization algorithm `BrownFullBasicInit` in OrdinaryDiffEq. If you do `@show solve(prob, Rodas4())[1]` you see:
```julia
(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.
---
### Author:  [@AnHeuermann](https://discourse.julialang.org/u/AnHeuermann)
#### Post date: [October 7, 2021, 10:54am UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312/4 "2021-10-07T10:54:46Z")
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
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).
---
### Author:  [@SCher](https://discourse.julialang.org/u/SCher)
#### Post date: [October 7, 2021, 2:37pm UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312/5 "2021-10-07T14:37:30Z")
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)
```julia
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)
```julia
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)
```julia
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.
---
### Author:  [@AnHeuermann](https://discourse.julialang.org/u/AnHeuermann)
#### Post date: [October 7, 2021, 3:32pm UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312/6 "2021-10-07T15:32:22Z")
I see how transforming the DAE `f(x(t),dx(t),y(t),u(t),p,t)` into two sets of equations
```julia
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](https://github.com/SciML/ModelingToolkit.jl/issues/988) 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](https://www.researchgate.net/publication/337672027_On_the_choice_of_initial_guesses_for_the_Newton-Raphson_algorithm)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 7, 2021, 9:57pm UTC](https://discourse.julialang.org/t/modelingtoolkit-index-reduction-and-linear-algebraic-loops/69312/7 "2021-10-07T21:57:02Z")
Yeah, this looks like it needs dummy derivative stuff.