MTK & linearization with analysis points

Specific Domains Modelling & Simulations
1

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:

image
image1332Ă—920 86.3 KB

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?

2

image
image919Ă—141 14.3 KB

_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
3

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

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?

5

OK. If I change the operating point description to:

op = Dict(v => sol_rx_s[v][end] for v in unknowns(sys_rx_s)[[1;6:8]])

thus only including the temperature and the mole numbers, I avoid the overdetermined system error. But when I attempt to linearize, I get another error message:

julia> matrices, _ = linearize(_sys_rx_s, inputs, outputs; op=op)
UndefVarError: `pb_e_sâ‚Šaâ‚Šx` not defined in `ModelingToolkitBase`
Suggestion: check for spelling errors or missing imports.

Stacktrace:
  [1] macro expansion
    @ C:\Users\Bernt\.julia\packages\Symbolics\29fkz\src\build_function.jl:143 [inlined]
  [2] macro expansion
6

Don’t know if I am making progress or not… I have looked at the “mtkcompiled” symbolic model equations:

image
image996Ă—819 70.6 KB

with unknowns:
image
image623Ă—558 29.5 KB

My understanding (I may be wrong…) is that mtkcompile reduces the original model _sys_rx_s to an index 1 DAE (sys_rx_s). Here, the reduced model has 4 “states” (i.e., differential equations that require initial values) and 4 algebraic variables that need “guess” values to ensure that the nonlinear solver is able to solve for the algebraic variables.

In principle, it should be relatively easy to linearize this model; in addition to the equations, I should just need operating values for the “states” and good guess values for the algebraic variables, + the operating values of the inputs. I have done that in the past for index 1 DAEs using Symbolics.jl. It should not be necessary to differentiate signals, unless those show up from the index reduction procedure. Here, there is the added difficulty that some functions in the model are not “traceable”; these have been registered symbolically.

Since I have added analysis points for two inputs, I assume that MTK.linearize knows that the operating value of these are zero, so I do not provide them.

I have tried to provide operating point values for the “states” (or rather: for 3 of the states + the temperature; that should in principle work), and guess values for the algebraic values:

op = Dict(v => sol_rx_s[v][end] for v in unknowns(sys_rx_s)[[4;6:8]])
op_guess = Dict(v => sol_rx_s[v][end] for v in unknowns(sys_rx_s)[[1:3;5]])

Next, I provide inputs and outputs:

inputs  = [_sys_rx_s.Tc_ap, _sys_rx_s.Twi_ap]
outputs = [rx_s.T, rx_s.xi]

Here, the “inputs” are the analysis points.

Then I try to do the linearization based on the composed model, i.e., before “mtkcompile”:

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

which crashes after 55 seconds on my lunar lake laptop, with error message:

UndefVarError: `pb_e_sâ‚Šaâ‚Šx` not defined in `ModelingToolkitBase`
Suggestion: check for spelling errors or missing imports.

Stacktrace:
  [1] macro expansion
    @ C:\Users\Bernt\.julia\packages\Symbolics\29fkz\src\build_function.jl:143 [inlined]
  [2] macro expansion

My model may be tricky to linearize as I think it is rather stiff (I include compressibility of liquid).

A handful of questions:

  1. Is my understanding above correct?
  2. If I instead of specifying temperature T as an operating point, choose to specify U (internal energy), I get warning about over-specified set of equations. In my reactor component function I specify initial value of the temperature, and not of the internal energy
# Algebraic Variables
        T(t) = T0,          [unit = u"K", description="Temperature"]

Does MTK remember this when I try to linearize the model? If so, if I also specify U as an operating point value, this could explain that the system then has too many equations…

  1. Is it considered bad practice to specify initial values in component functions? I realize that this may cause problems if I in addition specify initial values in the ODEProblem function. Is it better to specify guess values for values in a component function?
7

… Here is what Google AI Studio thinks about the key question on how to specify initial value for the simulation vs. operating point for linearization. I don’t know whether this is accurate, though:

The way ODEProblem handles initial values is the exact same mechanism linearize uses for the operating point.

To answer your question directly: Yes, it will try to use the default value of T from the component, but it treats it as a “Guess” (soft constraint) rather than a “Fixed Value” (hard constraint).

Here is exactly what happens inside the solver when you provide U in the u0 of an ODEProblem but leave out T:

1. The Initialization Hierarchy

When you call solve(prob), the first thing Julia does is solve the Initialization Problem. It needs to find a consistent set of values for every variable such that all algebraic equations (like your Energy Balance) are satisfied at t=0.

  • Fixed Variables: Since you specified U = U_{ss} and n = n_{ss} in the u0 of the ODEProblem, the solver “locks” these values. They are non-negotiable.
  • Free Variables: Since T is not in your u0, the solver marks T as an unknown that needs to be found.
  • The Guess: The solver looks at the reactor component and sees @variables T(t) = T0. It takes T_0 and says: “I will use this as my starting guess for the Newton-Raphson iteration.”

2. The Root-Finding Process

The solver then attempts to solve:

0 = f(T) = U_{fixed} - \text{thermo\_energy}(T, n_{fixed})
  1. It starts at T = T_0.
  2. It calculates the residual (U_{fixed} - f(T_0)).
  3. It adjusts T iteratively until the residual is zero.

3. The Danger (Why it might crash)

This is where your PC-SAFT thermodynamics make things difficult.

  • If your steady-state U corresponds to 85°C, but your component default T_0 is 20°C, the solver starts the iteration very far from the solution.
  • Because PC-SAFT is highly non-linear, the derivative \partial U / \partial T at 20°C might point the solver in the wrong direction, or lead it into a region where the EoS returns NaN (e.g., a phase change boundary).
  • Result: The initialization fails, and you get a Linearization failed to find consistent initial conditions or Dt <= dtmin error.

4. What happens if you specify both U and T?

If you put both U and T in the u0 of the ODEProblem:

  • MTK sees you have provided 2 values for a relationship that only has 1 degree of freedom (U \leftrightarrow T).
  • It will usually give you the “Overdetermined system” warning you’ve been seeing.
  • It then uses a Least Squares approach to find a compromise between your U and your T. If they are perfectly consistent (from your steady-state solution), there is no problem. If they are slightly inconsistent, it will shift them both a tiny bit to satisfy the EoS.

Summary for your Workflow:

If you want to start a simulation (or a linearization) from steady state using U as the primary input:

  1. Fixed Values (op or u0): Provide U and n.
  2. Guesses (guesses or u0 defaults): Provide the steady-state T.

Crucial Advice: The reason you were getting errors is that when you leave T to be handled by the “component default,” the solver is starting its search at the wrong temperature.

Does this explain why providing only U caused a conflict with the component’s T_0? The solver was trying to reconcile your “New U” with your “Old T Guess.”

AI Studio is somewhat off, since I actually specified guess values.

8

This is not your fault, it’s a bug in MTK.

1 Like
9

OK. I hope this is on the todo list. It is kind of crucial that it works (in my case) since x is a vector.

Here is what Google AI Studio says about it – and I don’t know if I trust the AI solution…

The pb_e_s₊a₊x Error (The “UndefVarError”)

This error is a Symbolics.jl code-generation bug.
The variable pb_e_s₊a₊x is a vector (mass fractions at the port). When linearize tries to build the Jacobian, it creates a “Runtime Generated Function” (RGF).

  • The compiler sees the vector-variable pb_e_sâ‚Šaâ‚Šx.
  • It fails to “unroll” this vector into individual scalar symbols (x_1, x_2, x_3).
  • It looks for a single variable named exactly pb_e_sâ‚Šaâ‚Šx in the ModelingToolkitBase namespace and fails.

How to bypass this immediately:
You are linearizing _sys_rx_s (unsimplified). This system contains hundreds of “copy” equations for mass fraction vectors. The linearizer is getting lost in the vector-symbolic forest.

The Solution: You must linearize the simplified system to prune these redundant vector identities.

# 1. Simplify the system first
sys_simple = structural_simplify(_sys_rx_s)

# 2. Re-map your OP to the simplified unknowns
# Use the solution object to find the values of whatever unknowns survived simplification
op_simple = Dict(v => sol_rx_s[v][end] for v in unknowns(sys_simple))

# 3. Reference inputs/outputs from the simplified system
inputs  = [sys_simple.Tc_ap, sys_simple.Twi_ap]
outputs = [rx_s.T, rx_s.xi]

# 4. Linearize
matrices, _ = linearize(sys_simple, inputs, outputs; op=op_simple)

Why this specific workflow is “Perfection”:

  1. Removes Vectors: structural_simplify turns vector-to-vector port connections into scalar-to-scalar equations. This eliminates the pb_e_s₊a₊x error because the “vector” variable no longer exists in the system’s core equations.
  2. Solves Overdetermination: By simplifying, you reduce the 182 equations to the 8 core equations you identified (U, n, T, \dots). Your op will now match the number of unknowns exactly.
  3. Speed: The 54-second runtime will likely drop to under 5 seconds.

I tested the AI “solution”, and it did not work.

10

Please do not post this kind of faulty output from Ai tools. Next iteration of Ai models will learn from this and perpetuate this incorrect information.

4 Likes