MTK & `@mtkbuild` -- structure of result?

[B-LIE]

Question on @mtkbuild model structure❓

Consider a DAE model described by @mtkmodel in the form:

$$ \frac{dx}{dt} = f(x,z,u) $$ $$ 0 = g(x,z,u) $$

Here, $x$ is the "differential variable", $z$ the "algebraic variable", and $u$ an input which is set equal to a register_symbolic function. Assume that the user is inexperienced, and may have created some differential equations by differentiating simple algebraic equations, e.g., differentiating $x_1 = x_2 + 1$.

Next, macro @mtkbuild is applied, leading to a transformed model, say, in the form:

$$ \frac{d \xi}{dt} = \phi (\xi,\zeta, u) $$ $$ 0 = \gamma(\xi, \zeta, u) $$

Question 1: Is this transformed model in Index 1 DAE form? In other words, is $\partial \gamma/ \partial \zeta$ non-singular?

Question 2: Will in, in general, $\dim \xi \le \dim x$, or is that uncertain? [For sure, if "fake" differential equations from differentiating $x_1 = x_2 + 1$ are removed and made into algebraic equations, $\dim \xi$ should be smaller than $\dim x$. But perhaps new differential equations are introduced in the index reduction process?]

In addition to the above model, some observed relations are included. These are of form:

$$ y = h(\xi, \zeta, u) $$

Question 3: Will the set of variables ${\mathcal S}(x,z,u) \subseteq {\mathcal S}(\xi, \zeta, y, u)$?

[AayushSabharwal]

1. Yes, simplified models are reduced to index-1 form.
2. I think that for the form you've proposed, this should be true.
3. Yes.

[B-LIE]

Follow-up related to linear sub-constraints.

Given the transformed Index 1 DAE form:

$$ \frac{d\xi}{dt} = \phi (\xi, \zeta, u) $$ $$ 0 = \gamma(\xi, \zeta, u) $$

I guess $\xi$ can be considered the "state" (i.e., the minimal information that has to be provided together with future inputs $u$ to solve the model).

• My experience (solving DAEs, not linearization) is that instead of providing initial values for $\xi$, provided that there is a linear/invertible mapping between some observed quantities $y$ (e.g., subset of $y$) and the same number of $\xi$ (subset of $\xi$), then MTK is able to use these known observed values and find the compatible $\xi$, and the initialization is accepted.

Question 4: Is this correct and general?

In principle, the algebraic quantities $\zeta$ don't need to be provided, since in principle they can be inferred from the algebraic constraints $0 = \gamma(\xi, \zeta, u)$ when $\xi$ and $u$ are known.

However, there is a chance that the algebraic constraints have non-unique solutions.

Question 5: Is it correct to say that if some of the algebraic constraints have unique solutions given $\xi$ and $u$, it is not necessary to provide "guess" values for the algebraic variables $\zeta$. However, if some algebraic constraints have multiple valid solutions, then "guess" values for them are required. A guess value could in principle be in the form of a range of valid values, or alternatively a guess point value which is such that the nonlinear solver converges to the correct value $\zeta$?

Question 6: Would the above principles also be valid for providing operational point in the case of linearization?

---

OK -- hope you guys can clarify these questions. Understanding them will make it easier (for me) to understand how to correctly initialize a system.

[AayushSabharwal]

4. Yes. If you have k differential variables $\xi$ and instead provide initial conditions for k observed variables $y$ (or a mix of differential and observed variables, totaling to k) and there is a mapping between the values you provided initial conditions for and $\xi$, MTK can figure out what $\xi$ is and from there figure out $\zeta$. Note that the mapping doesn't have to be linear/invertible. MTK is capable of doing a nonlinear solve to figure out $\xi$. This is where guess values come into play. E.g. if x ~ y^2 is an observed equation where y $\in\xi$, and x => 4.0 is an initial condition, then y could be $\pm$ 2. If y has a positive guess, it'll be 2 and if it is negative then -2.
5. Not necessarily - the requirement for a guess depends on whether MTK can figure out how to symbolically invert the mapping. For example, it currently can't symbolically solve a linear system of equations. So if x + y ~ 4, x - y ~ 7 (somehow) end up being initial conditions, it can't symbolically solve for x and y and instead will use a nonlinear solver
6. Yes, these are valid for linearization

[AayushSabharwal]

I accidentally sent the message early. Edited it now with the full information.

[B-LIE]

Great. To me, it is useful to structure the understanding like the above. If you have ideas on how to make this initialization process simpler for the user, I look forward to see that :-). Some "tricky" points are (from my experience):

• Setting default values is useful (from my perspective), but if I leave some code for a while or use the code of others, it may be difficult to remember/find which variables in code have been given defaults. So a simple way to figure out which values have defaults and set these to nothing would be welcome.
• Because the translation process introduce new variable, say Vd_d_t in my case, I don't know how to "address" that variable. Attempting to specify hp_ns_oL.Vd_d_t => 0 doesn't work. Also, D(hp_ns_oL.Vd_d) => 0 fails. It is not clear to me how I should refer to such a variable.
• Regarding linearization, one would normally always do that at a point where the time derivatives are zero; if not the matrices will be time varying. If the system is asymptotically stable, the simplest way to find the operating point is probably to simulate for a long time. That should provide good initial guesses for the algebraic variables! (But how to address variables such as Vd_d_t?). Only for oscillatory or unstable systems would one use a nonlinear solver to find the equilibrium point, I'd think.

[AayushSabharwal]

So a simple way to figure out which values have defaults and set these to nothing would be welcome.

defaults(sys) returns a dictionary mapping variables to their default values. If sys is completed, then this dictionary can be mutated to update the defaults.

Because the translation process introduce new variable, say Vd_d_t in my case, I don't know how to "address" that variable. Attempting to specify hp_ns_oL.Vd_d_t => 0 doesn't work. Also, D(hp_ns_oL.Vd_d) => 0 fails.

D(hp_ns_oL.Vd_d) => 0 should work. Do you have an example I can look into?

[B-LIE]

OK -- here is what happens:

So -- it seems like Vd_d_t and Dt(Vd_d) are two different entities. [I have imported D_nounits as Dt because I use D to signify diameter`.

[AayushSabharwal]

Ah yeah, so D(Vd_d) is translated under the hood into Vd_d_t, but if you have both in the map it'll give priority to the latter.

[B-LIE]

Ah yeah, so D(Vd_d) is translated under the hood into Vd_d_t, but if you have both in the map it'll give priority to the latter.

Two questions:

1. "priority to the latter" -- does that mean "priority to the first" (of D(Vd_d) and Vd_d_t), or "to the latter" (of V_d_t(t) and Differential(t)(V_d(t)) ?
2. If I had imported D_nounits as D instead, would I then have got D(hp_ns_oL.Vd_d) -> Vd_d_t?

[AayushSabharwal]

1. It'll give priority to Vd_d_t.
2. What you import it as doesn't matter. The object is the same, just the binding to it in your REPL changed.

[B-LIE]

defaults(sys) returns a dictionary mapping variables to their default values. If sys is completed, then this dictionary can be mutated to update the defaults.

OK. If I do defaults(hp_ns_oL), I get a dictionary of 116 entries. Not the easiest to digest.

[AayushSabharwal]

Yeah it's all the defaults you provide + some extras that MTK needs

[B-LIE]

1. It'll give priority to Vd_d_t.

OK. So if I only do:

guesses(hp_ns_LoL)[Dt(hp_ns_oL.Vd_d)] = sol_ss(tspan_ss[2], idxs=Dt(hp_ns_oL.Vd_d))

resulting in

Dict{Any, Any} with 3 entries:
  w(t)                     => 55.9432
  p_d1(t)                  => (H_i + H_p + H_r)*g*rho
  Differential(t)(Vd_d(t)) => 8.44046e-15

then Differential(t)(Vd_d(t)) will be interpreted as Vd_d_t? And if I in addition set a guess value of Vd_d_t, then that will take presedense?

[AayushSabharwal]

Yes. Differential(t)(Vd_d(t)) is the application of the operator Differential(t) to the symbolic variable Vd_d. It represents the differential of Vd_d with respect to t. During simplification of the model, the dummy derivatives algorithm determines that some differential variables (in this case D(Vd_d)) are actually algebraic variables in disguise and introduces a new variable Vd_d_t to solve this. D(Vd_d) and Vd_d_t are different symbolic variables, but represent the same quantity.