# MTK & linearization with analysis points

**URL:** <https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165>\
**Category:** Modelling & Simulations\
**Created:** [March 12, 2026, 3:37pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165 "2026-03-12T15:37:40Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 12, 2026, 3:37pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/1 "2026-03-12T15:37:41Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/b/a/baf3bc3e1b5f967c3ab0b5c3c37ae51bb52b2fc4.jpeg)

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

```julia-auto
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:

```julia-auto
@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:

```julia-auto
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:

```julia-auto
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?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 12, 2026, 3:59pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/2 "2026-03-12T15:59:08Z")

</div>

> [@BLI](#):
>
> But which “model” should I use? `sys_rx_s`? `_sys_rx_s`?

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/8/d81e95efc3e83897ffc7beb84a1fdf58945594e6.png)

> [@BLI](#):
>
> what should be the _inputs_ for my two analysis points `:Twi_ap` and `:Tc_ap`?

`_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.

> [@BLI](#):
>
> And how do I specify the operating point?

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, ...)`

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 12, 2026, 5:00pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/3 "2026-03-12T17:00:06Z")

</div>

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:

```julia-auto
prob_rx_s = ODEProblem(sys_rx_s, [], tspan; use_scc=false)

```

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

```julia-auto
[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`:

```julia-auto
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):

```julia-auto
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:

```julia-auto
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
...

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 12, 2026, 5:24pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/4 "2026-03-12T17:24:12Z")

</div>

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

> <https://github.com/SciML/ModelingToolkit.jl/issues/4159>
>
> This feature exists in ControlSystemsMTK but it has proven pretty much impossibl…e to maintain, so I'm asking for it to exist in MTK directly instead. 
> 
> It would be very useful if one could pass a solution object as the operating point to \`linearize\`, and the operating point would be automatically constructed from the solution and the specified time point. It should also be possible to provide a vector of time points in order to linearize along a trajectory.

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?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 12, 2026, 6:50pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/5 "2026-03-12T18:50:29Z")

</div>

OK. If I _change_ the operating point description to:

```julia-auto
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
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

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 13, 2026, 5:46pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/6 "2026-03-13T17:46:52Z")

</div>

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

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/5/4558febe8549af8472199ce1f1db1a5634b2dd80.png)  
with unknowns:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/3/b39cc6a0383e6842ee2166bfbc66acdc88a69f90.png)

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:

```julia-auto
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:

```julia-auto
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”:

```julia
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:

```julia-auto
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

```julia
# 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?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 13, 2026, 8:11pm UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/7 "2026-03-13T20:11:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 14, 2026, 6:44am UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/8 "2026-03-14T06:44:37Z")

</div>

> [@BLI](#):
>
> which crashes after 55 seconds on my lunar lake laptop, with error message:
> 
> ```julia-auto
> 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
> 
> ```

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

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 14, 2026, 9:19am UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/9 "2026-03-14T09:19:36Z")

</div>

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.

```julia
# 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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 14, 2026, 11:21am UTC](https://discourse.julialang.org/t/mtk-linearization-with-analysis-points/136165/10 "2026-03-14T11:21:08Z")

</div>

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.