# ModelingToolkit - the concept of "state"
**URL:** https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304
**Category:** Modelling & Simulations
**Created:** [October 26, 2022, 2:37pm UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304 "2022-10-26T14:37:06Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [October 26, 2022, 2:37pm UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/1 "2022-10-26T14:37:06Z")
I’m a little confused about the use of the word “state” (not only in ModelingToolkit). In classic dynamic systems theory, the _state_ is the minimum information necessary to describe current time of a dynamic system, so that the model evolution (with model parameters) is uniquely described when the _state_ + _future inputs_ are specified.
Macroscopic models based on balance laws can often be described by _semi-explicit DAEs_ of form:
dx/dt = f(x,z,u;t)
0 = g(x,z,u;t)
Sometimes, an output is added:
y = h(x,z,u;t)
Here, the unknown variables are x and z (and possibly y), while u is a known input.
It may be convenient to classify “x” as a differential (or: differentiated?) variable, and z as an algebraic variable.
Some documentation (not only MTK) suggests that (x,z) is the state. But this is not so if we use the classic definition of “state” as the minimum necessary information.
In DAE literature, it is common to denote (x,z) as the _descriptor_ of the system (my books are still in boxes after changing office…, so I can’t find references right now); alternatively, the _unknown variables_.
For examples of DAEs, see, e.g.,
- Brenan, Campbell, Petzold (1987). Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations (see Amazon)
- Ascher, Petzold (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, SIAM
Ascher and Petzold, in Example 9.7 gives a stylized model of a 2D pendulum, constrained to move on a circle:
q1’ = v1
q2’ = v2
v1’ = -L.q1
v2’ = -L.q2-g
0 = q1^2 + q2^2 - 1
where L(t) is unknown together with q1, q2 (positions) and v1, v2 (velocities). What is the state of this system? Is it the variables (q1, q2, v1, v2, L)? Or something else? The state should probably be the “differential variables” after index reduction.
Using MTK, the model is:
and function `structural_simplify` leads to:
Inspecting the resulting simplified model from `structural_simplify`, some equations must be missing. Yes, we can compute q1 when q2 is known (fourth equation, with an assumption of sign), but we cannot simultaneously compute q2\_tt and L from the third equation. And what about v1? In Ascher and Petzold, they propose both:
0 = q1.v1 + q2.v2
-L = q2.g - v1^2 - v2^2
The first of these two additional equations gives us v1, while the second gives us L. Finally, we can compute q2\_tt from the third equation provided by MTK.
Both of these additional equations are necessary, but `structural_simplify` doesn’t spit them out.
In **summary** , it seems to be sufficient to specify initial values for q2 and v2, thus the system has 2 states - based on the classical definition of state from (dynamic) systems theory. It follows that the state is (q2,v2), or any invertible transformation of these two. With the key point being that the system has 2 states.
The original system had 5 unknown variables (q1,v1,q2,v2,L) [sometimes denoted the _descriptor_], but index reduction added a sixth unknown q2\_tt.
**Summary** of _summary_: I’m not claiming that the classical definition of “state” in dynamic systems theory is the only valid definition. Maybe other fields use the concept of state differently? Most likely, many users of MTK (and SciML tools) will have a background in classical dynamic systems theory, so it might be useful to keep this in mind in documentation.
[And… am I doing something wrong wrt. `structural_simplify` in that 2 equations appear to be missing?]
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 26, 2022, 3:48pm UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/2 "2022-10-26T15:48:23Z")
> [@BLI](#):
>
> Some documentation (not only MTK) suggests that (x,z) is the state. But this is not so if we use the classic definition of “state” as the minimum necessary information.
The splitting of `x` vs `z` is not necessarily possible in higher index nonlinear scenarios. We thus say that it’s MTK’s job to figure out which of the variables are observed values, states, and it’s MTK’s job to determine which ones are what and whether to treat certain objects as diffferential or algebraic (which is really just a distinction for initialization).
> [@BLI](#):
>
> Both of these additional equations are necessary, but `structural_simplify` doesn’t spit them out.
`full_equations(sys)`. Those equations are not necessary because they are eable to be easily eliminated.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [October 27, 2022, 7:08am UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/3 "2022-10-27T07:08:28Z")
`full_equations` did the trick…
```julia
@parameters g
@variables t, q1(t), q2(t), v1(t), v2(t), λ(t)
D = Differential(t)
eqs = [D(q1) ~ v1,
D(q2) ~ v2,
D(v1) ~ -λ*q1,
D(v2) ~ -λ*q2 - g,
0 ~ q1^2 + q2^2 - 1]
@named sys = ODESystem(eqs)
sys = structural_simplify(sys)
full_equations(sys)
```
leads to:
The resulting model has two differential variables and has been index-reduced. So from classical dynamic systems theory, the system has _two states_ = required initial conditions, and this is also true for the original system with 4 differential equations and one algebraic constraint.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [October 27, 2022, 12:35pm UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/4 "2022-10-27T12:35:15Z")
[**fixed bug** detected by @albheim]
OK— perhaps the thread becomes too long. But I have pondered more on this 2D “pendulum” model from Ascher and Petzold… to recapture:
```julia
@parameters g
@variables t, q₁(t), q₂(t), v₁(t), v₂(t), λ(t)
D = Differential(t)
eqs = [D(q₁) ~ v₁,
D(q₂) ~ v₂,
D(v₁) ~ -λ*q₁,
D(v₂) ~ -λ*q₂ - g,
0 ~ q₁^2 + q₂^2 - 1]
@named sys = ODESystem(eqs)
```
is the encoding of the model:
If I do `structural_simplify` and show the `full_equations`,
```julia
sys_simp = structural_simplify(sys)
full_equations(sys_simp)
```
After some “tedious” manipulation and by eliminating \lambda, I find that — with states q\_2 and v\_2, the model can be re-formulated as:
It is obvious that this system has two states, and it should be trivial to solve it with any ODE solver… except we will have problems for q\_1 =0… Perhaps a solution could be to replace q\_1 ^2 \rightarrow q\_1^2 + \epsilon. If we need \lambda, we can find it from:
If I try to solve the above “manipulated” system as an ODE by coding my own function for the differential equations, I get
```julia
function pendulum!(du,u,p,t)
g = 9.81
q2 = u[1]
v2 = u[2]
q1_2 = 1-q2^2
du[1] = v2
du[2] = -q2*v2^2/q1_2 - g*q1_2
end
u0 = [sqrt(0.5); -1.0]
tspan = (0.0,10.0)
prob = ODEProblem(pendulum!,u0,tspan)
sol = solve(prob)
```
Perhaps I need to tweak tolerances? Specify another solver?
For comparison, here is the similar Modelica code [Sundials solver?]:
```julia
model pendulum
constant Real g = 9.81;
parameter Real q20 = sqrt(0.5);
parameter Real v20 = -1;
//
Real q2(start = q20, fixed = true);
Real v2(start = v20, fixed = true);
Real q12;
Real q1;
equation
der(q2) = v2;
der(v2) = -q2*v2^2/q12 - g*q12;
q12 = 1 - q2^2;
q1^2 = q12;
end pendulum;
```
which gives the solution without problems:
[except that Modelica chooses the wrong root `q1` (purple) from `q1^2` (red) at times…]
If I try to solve the problem using the MTK model, it seems like I need to specify 4 initial values (although there clearly is only 2 independent initial conditions, i.e., 2 states…)
---
### Author:  [@albheim](https://discourse.julialang.org/u/albheim)
#### Post date: [October 27, 2022, 4:21pm UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/5 "2022-10-27T16:21:06Z")
> [@BLI](#):
>
> `du[u2] = -q2*v2^2/q1_2 - g*q1_2`
I think `2` instead of `u2` here for the crash, though fixing it there is still the instability making for an aborted solve.
If you only give two initial values, you can’t really decide q1 though? I guess the simulator will decide positive or negative initial value for you, and I guess that does not really matter, if you wanted it the other way you can just negate the solution.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [October 27, 2022, 4:22pm UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/6 "2022-10-27T16:22:15Z")
Oops…
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [October 28, 2022, 7:06am UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/7 "2022-10-28T07:06:42Z")
> [@albheim](#):
>
> If you only give two initial values, you can’t really decide q1 though?
Yeah, that is probably true. I do think the solution is “unphysical” considering it a pendulum, though. `q1` should not bounce back from `q1 = 0` coming from above. I may have to use callback and use the velocity `v1` (not computed above) to determine the physically correct root.
Anyways, my main “concerns” are:
- What is the minimal initial information/initial values to simulate the system? Answer: 2, which means there are 2 _states_.
- This traditional meaning of the word _state_ for dynamical systems is at odds with the way the word is used in MTK (e.g., in function `states`).
- My MTK version of the model gives an error message if I try to specify the 2 required initial values. It seems like I need to specify one value for each of the 4 “states”. [In my view it would be better if MTK talked about _unknowns_ or the _descriptor_ instead of introducing a new meaning to the word “state”].
- Modelica has no problem solving the model when I specify 2 initial values.
- Does MTK write the equations of this example into a DAE with mass matrix, and then requires two initial states and two initial guesses for the algebraic variables?? Is that why 4 “states” are required?
- In Modelica, with 2 states, it suffices with 2 initial values (syntax: `q2(start = q20,fixed = true)` which implies that the initial value _has_ to be adhered to) and no initial guesses for the algebraic variables — for the above “pendulum” case. In more complex cases (i.e., more complicated algebraic equations), initial guesses which are considered as suggestions can be used (syntax: `q12(start = q120,fixed = false)`) — then Modelica solvers will use some least squares choice for compatible initial guesses.
- So … if my “guess” that MTK writes the above model in (singular) mass matrix form is correct, and that the “states” are the unknowns, and that it is required to specify initial values and initial guesses for all unknowns — the question is: must the initial guesses completely satisfy the algebraic constraints? [Not necessary in Modelica.]
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [October 28, 2022, 8:07am UTC](https://discourse.julialang.org/t/modelingtoolkit-the-concept-of-state/89304/8 "2022-10-28T08:07:53Z")
Hm… if I implement the `structural_simplify` model from MTK in Modelica, the result is identical to what I get with my own further simplification (simulation results above) — except that it crashes after 8-9 seconds.
However, if I implement the _original_ model as given in Ascher and Petzold:
```julia
model pend3
constant Real g = 9.81;
parameter Real q20 = sqrt(0.5);
parameter Real v20 = -1;
//
Real q2(start = q20, fixed = true);
Real v2(start = v20, fixed = true);
Real q1;
Real v1;
Real lambda;
equation
der(q1) = v1;
der(q2) = v2;
der(v1) = -q1*lambda;
der(v2) = -g - q2*lambda;
0 = -1 + q1^2 + q2^2;
end pend3;
```
I get the physically correct result (also for `q1`):
Observe that this time, `q1` (green curve) passes through zero and can become negative.
NOTE: I have used OpenModelica, which also does index reduction, and uses an ODE solver (by default). Perhaps a different index reduction algorithm?