MTK & linearization with analysis points

I am playing around with MTK v11 (v11.14) and trying to understand connectors, components + signals. My code now work, and I can specify signals. I use a chemical reactor (liquid phase) as an example, where the temperature evolves as:

The model involves comparison of (i) a Clapeyron.jl EoS vs. a crude model (IS), and (ii) a model with compressible liquid (solid) vs. one with specified reactor pressure (dash).

So, now I want to try to find a linear approximation of my model. As often in the past, I have some problem in linearizing the model, and I don’t know whether this is due to updates in MTK or not.

Here are some snapshots of how I (a) included the analysis points, e.g.,

connect(T_src_w.y, :Twi_ap, pb_H2O_s.T_set),

inserts the analysis point :Twi_ap between signal T_src_w.y and a component port pb_H2O_s.T_set ,

next, (b) how I compose the model:

@named __sys_rx_s = System(eqs_rx_s, t)
@named _sys_rx_s = compose(__sys_rx_s, [p_src_w, T_src_w, p_src_eo, T_src_eo,
            u_v_src_w, u_v_src_eo, u_v_src_e, T_src_c, p_src_e, T_src_e,
            pb_H2O_s, pb_EO_s, pb_e_s, va_H2O_s, va_EO_s, va_e_s, rx_s])
sys_rx_s = mtkcompile(_sys_rx_s)

For simulation, I create and solve the numeric problem:

tspan = (0.0, 100).*60

prob_rx_s = ODEProblem(sys_rx_s, [], tspan; use_scc=false)
sol_rx_s = solve(prob_rx_s,Rodas5P(); reltol=1e-7, abstol=1e-7)

[The use_scc keyword is to get around a limitation in NonlinearSolve, I think: the case is rather stiff due to including liquid compression.]

Simulating for 100 minutes leads to steady state conditions.

Now, for the linearization part, which should have the form:

linearize(model, inputs, outputs; op=op)

But which “model” should I use? sys_rx_s? _sys_rx_s? structural_simplify(_sys_rx_s)?

Suppose my outputs are [sys_rx_s.T, sys_rx_s.xi]… (or something else?), what should be the inputs for my two analysis points :Twi_ap and :Tc_ap? Should I set inputs to [:Twi_ap, :Tc_ap], or something else?

And how do I specify the operating point?

_sys_rx_s.Twi_ap, if this analysis point is at the top level. :Twi_ap probably also works but we are moving away from the :x kind of symbol in favor of standard MTK-style access.

The same way as you do for ODEProblem(sys_rx_s, [], tspan; use_scc=false), which here is an empty array. Use a Dict(var => value, ...)

1 Like

OK – my mtkcompiled sys_rx_s has 4 algebraic + 4 differential equations. The differential equations are for internal energy U, and the number of moles n[1:3] in the mixture. Of the algebraic equations, one is for temperature, and 3 appear to be related to three signals.

When I do:

prob_rx_s = ODEProblem(sys_rx_s, [], tspan; use_scc=false)

the empty bracket ([]) could be 4 pairs of initial value specifications, e.g.,

[rx_s.T=>50+273.15, rx_s.n[1]=>10, rx_s.n[2]=>3 rx_s.n[3]=>5.0]

[Note: this is not a dictionary.] Next, I would need to specify 4 “guess values”, which are specified with keyword guesses:

guesses = [rx_s.U => 0.0, va_e_s.T_i=>20+273.15, 
    va_H2O_s.T_i=>60+273.15, va_EO_s.T_i=>60+273.15]

But in my case, I have specified these values in the components.


So – how do I specify the operating point op?

  • With a dictionary (instead of vector of pairs)… OK.
  • Which variables do I need to specify for the uncompiled system _sys_rx_s? The same as for the mtkcompiled system sys_rx_s? Do I make a distinction between initial values and guesses?
  • Since I am using analysis points, I assume I do not need to specify the inputs, since the analysis points essentially are additional perturbations to the input signals with nominal value zero, and MTK knows this?

If I do the following (the solution sol_rx_s has been run to steady state):

op = Dict(v => sol_rx_s[v][end] for v in unknowns(sys_rx_s))
inputs  = [_sys_rx_s.Tc_ap, _sys_rx_s.Twi_ap]
outputs = [rx_s.T, rx_s.xi]

matrices, _ = linearize(_sys_rx_s, inputs, outputs; op=op)

I get the error message:

Warning: Initialization system is overdetermined. 7 equations for 3 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here.
│ 
│ 
│ To suppress this warning, pass `warn_initialize_determined = false`. To turn this warning into an error, pass `fully_determined = true`.
â”” @ ModelingToolkit C:\Users\Bernt\.julia\packages\ModelingToolkit\4uiLn\src\initialization.jl:22
TypeError: in typeassert, expected DataType, got Type{Union{Nothing, Real}}

Stacktrace:
  [1] safe_eltype
    @ C:\Users\Bernt\.julia\packages\SymbolicUtils\lqSte\src\methods.jl:34 [inlined]
  [2] promote_symtype
...

You can pass a keyword argument guesses as well if that’s needed.

There is currently no good way of using an existing solution as the operating point, an open issue about exists

The initialization system is likely over determined because your approach specifies the algebraic variables as well. I’m not sure why it errors though, @cryptic.ax?