ModelingToolkit - the concept of "state"

[BLI]

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?]

[ChrisRackauckas]

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

[BLI]

full_equations did the trick…

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

[BLI]

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

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

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

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?]:

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

[albheim]

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

[BLI]

Oops…

[BLI]

@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.]

[BLI]

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:

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?