Not an issue... but trying to understand MTK

[B-LIE]

I'm gradually learning to use MTK by converting some examples I have from Modelica to MTK. It is great fun (:grinning:). As part of this, I'd like to understand what I'm doing (:open_mouth:); I'm not reporting any issues -- it is more to get some clarification from more experienced users...

Below, I attempt to describe my understanding in sections 1. Variables and equations, 2. Model and index reduction. Observables, and 3. MTK vs. Modelica. I close with 4. Questions. Sorry for the lengthy/wordy/technical post.

Variables and equations

I define my independent variable t as well as dependent variables which I will refer to as unknowns, let us denote them $\upsilon$ (Greek "upsilon", for unknown).

@variables t υ(t)

I have "cheated" in that there can be many unknowns.

My model will also have a number of externally driving "force" functions $\phi$ (Greek "f", for force) which I define using the @register_symbolic macro:

@register_symbolic ϕ(t)

Again, I'm simplifying: I need one line of register_symbolic for each element in $\phi$.

I finally write my equations:

D = Differential(t)

eqs = [ D(υ) .~ F(υ, ϕ(t)) ;
            G(υ, ϕ(t)) .~ 0] 

where the equation is meant to be a vector of equations.

Model and Index reduction. Observables

With the above, I create the model as:

mod = ODESystem(eqs)

I then simplify the structure of the model, including index reduction to a DAE of index 0 (= ODE), or a DAE of index 1 using function structural_simplify:

mod_simp = structural_simplify(mod)

The resulting simplified model has the generic form: $$ \frac{d \varsigma}{d t} = f\left(\varsigma,\bar{\varsigma},\phi(t)\right) \ 0 = g\left(\varsigma, \bar{\varsigma}, \phi(t) \right) $$ If the result is an index 0 DAE (= ODE), then $\bar{\varsigma}$ doesn't exist, and $0=g(\cdot)$ vanishes. If the result is an index 1 DAE, then $\partial g/\partial \bar{\varsigma}$ is always invertible. For index 1 DAEs, an ODE solver with mass matrix is required.

In MTK, the state $\sigma$ (Greek "s", for state) is then defined as $\sigma = [\varsigma, \bar{\varsigma}]$.

Initial value for $\sigma$ must be specified by the user when converting the symbolic mod_simp to a numeric model using ODEProblem(), which is then solved.

If we go back to unknowns $\upsilon$, the state $\sigma$ is a "subset" of the unknowns $\upsilon$: in the process of building the model/structural simplifcation, MTK has "stripped off" parts of the model with some unknowns $\omega$ denoted observables (Greek "o", for observable); these observables can be expressed as $$ \omega = h(\sigma, \phi(t)) $$

The observables $\omega$ are not part of the model that is solved numerically. However, MTK keeps track of them, and they can easily be evaluated after the DAE has been solved.

MTK vs. Modelica

Modelica tools do something similar to what MTK does wrt. structural simplification, but normally under the hood.

In principle, it is sufficient to provide initial values for $\varsigma$ in order to solve the model $$ \frac{d \varsigma}{d t} = f\left(\varsigma,\bar{\varsigma},\phi(t)\right) \ 0 = g\left(\varsigma, \bar{\varsigma}, \phi(t) \right) $$ With value for $\varsigma(t=0)$ specified, the second equation (algebraic equation) can then be solved to find $\bar{\varsigma}(t=0)$, and then $\sigma(t=0)$ is known. This is why some fields within dynamic systems refer to $\varsigma$ as the state [in contrast to MTK, which refers to $\sigma$ as the state].

In practice, this is slightly more complicated in that the algebraic constraint may admit multiple solutions.

Modelica

In Modelica, it is sufficient to provide initial value for $\varsigma$. Then it is up to the solver to, if necessary, iterate on the equation $0 = g(\varsigma, \bar{\varsigma}, \phi(t))$ in order to find $\bar{\varsigma}$. However, the user may also provide an initial guess for $\bar{\varsigma}$. The user can then specify whether the initial value is meant to be enforced (attribute fixed = true), or meant to be an approximate guess (attribute fixed = false) -- attribute fixed can be used for both $\varsigma$ and for $\bar{\varsigma}$.

In practice, it may be necessary to provide initial guesses for $\bar{\varsigma}$ of the algebraic constraint, e.g., if it has multiple solutions.

In principle, I guess, all initial values can be specified as fixed = false, and then Modelica uses some least squares algorithm to make sure that the algebraic contraint $0 = g(\varsigma,\bar{\varsigma},\phi(t))$ is satisfied.

MTK

As far as I understand it, the user has to specify initial values for all of $\sigma = [\varsigma, \bar{\varsigma}]$.

What is not clear to me is whether $\sigma$ has to be specified so that the algebraic constraint $0 = g(\varsigma, \bar{\varsigma},\phi(t)) = g(\sigma, \phi(t))$ is satisfied, or whether it is sufficient with an approximate satisfaction of the algebraic constraint, and that the solvers has an initial phase where it iterates on $\sigma$ in order to satisfy the algebraic constraint.

Questions

1. Is there anything wrong in my "understanding" of the operation of MTK? [Probably...]
2. For index 1 DAEs, is it required to use an ODE solver that supports mass matrix? [...when I convert it to an ODESystem...]
3. Is there a list of ODE solvers with support for mass matrix?
4. The initial state ($\sigma$ in my notation...) -- is it required that the algebraic constraints ($0 = g(\sigma, \phi(t))$ are satisfied exactly?

[xtalax]

as for 3, you can find a list of DAE mass matrix solvers here https://docs.sciml.ai/DiffEqDocs/latest/solvers/dae_solve/

[B-LIE]

Ah. Thanks, @xtalax . I was looking under ODE solvers...

[ChrisRackauckas]

The initial state ( in my notation...) -- is it required that the algebraic constraints ( are satisfied exactly?

No, an initialization sequence will be ran if not. https://github.com/SciML/ModelingToolkit.jl/pull/2403 will improve this.

For index 1 DAEs, is it required to use an ODE solver that supports mass matrix? [...when I convert it to an ODESystem...]

Yes if you use ODEProblem. If you use DAEProblem then it needs a DAE solver. Alternative forms may be introduced in the future.

Is there anything wrong in my "understanding" of the operation of MTK? [Probably...]

Probably.