# MTK & Linearization

**URL:** <https://discourse.julialang.org/t/mtk-linearization/95803>\
**Category:** Modelling & Simulations\
**Created:** [March 9, 2023, 2:59pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803 "2023-03-09T14:59:35Z")\
**Posts on this page:** 20\
**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 9, 2023, 2:59pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/1 "2023-03-09T14:59:35Z")

</div>

### MTK documentation case

I’m looking at the linearization example in the MTK documentation:

```julia
using ModelingToolkit
@variables t x(t)=0 y(t)=0 u(t)=0 r(t)=0
@parameters kp = 1
D = Differential(t)

eqs = [u ~ kp * (r - y) # P controller
       D(x) ~ -x + u # First-order plant
       y ~ x] # Output equation

@named sys = ODESystem(eqs, t)

```

**Question 1** : statements like `@parameters kp = 1` and `@variables ... x(t)=0`, etc. – (i) is that a way to give _default_ values to the variables in question? (ii) Can these values be replaced when instantiating the numeric model via function `ODEProblem(..., par)`?

**Question 2** : the above code has 4 unknowns and 3 equations – is it a prerequisite for the following _linearize_ function that you have an “input” variable that has not been determined? In other words: if you specify nu input variables in the linearization, you have to have nu fewer equations than unknowns?

Next, you do the _linearization_:

```julia
matrices, simplified_sys = linearize(sys, [r], [y])

```

**Question 3** : (i) do you have to use the full model when you linearize? Or can you use the model found by `structural_simplify()`? (ii) OK – I tried to do `structural_simplify()` of this model, the simplification crashes because there are more unknowns than equations…

### Trying to linearize my own model…

I have tried to linearize my own model…

```julia
@variables t m(t) ṁ_i(t) ṁ_e(t) V(t) h(t)
@parameters ρ A K h_ς
Dt = Differential(t)

@register_symbolic ṁ(t)

eqs = [Dt(m) ~ ṁ_i-ṁ_e,
        m ~ ρ*V,
        V ~ A*h,
        ṁ_e ~ K*sqrt(h/h_ς),
        ṁ_i ~ ṁ(t)]

@named tank = ODESystem(eqs)

tank_simp = structural_simplify(tank)

```

If I try to linearize either `tank` or `tank_simp`, it fails in both cases; e.g.

```julia
linearize(tank, [ṁ], [h]; op = Dict(m=>2.4, ṁ=>2))

```

fails miserably.

**Question 4** : what is it that I do that is wrong?

---

<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 9, 2023, 3:08pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/2 "2023-03-09T15:08:52Z")

</div>

> [@BLI](#):
>
> (i) is that a way to give _default_ values to the variables in question?

yes

> [@BLI](#):
>
> (ii) Can these values be replaced when instantiating the numeric model via function `ODEProblem(..., par)`?

yes

> [@BLI](#):
>
> is it a prerequisite for the following _linearize_ function that you have an “input” variable that has not been determined?

yes, but you can make use of [analysis points](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/) to analyze fully connected models as well.

> [@BLI](#):
>
> do you have to use the full model when you linearize? Or can you use the model found by `structural_simplify()`?

You have to use the full model, at least for now.

> [@BLI](#):
>
> `[ṁ]`

This is a function, you can only linearize w.r.t. a variable. In the example below, I added default value 1 to all variables that were unspecified

```julia

@variables t m(t) ṁ_i(t) ṁ_e(t) V(t) h(t)
@parameters ρ=1 A=1 K=1 h_ς=1
Dt = Differential(t)

@variables ṁinp(t)=1

eqs = [Dt(m) ~ ṁ_i-ṁ_e,
        m ~ ρ*V,
        V ~ A*h,
        ṁ_e ~ K*sqrt(h/h_ς),
        ṁ_i ~ ṁinp]

@named tank = ODESystem(eqs)

linearize(tank, [ṁinp], [h]; op = Dict(m=>2.4, ṁ=>2))

julia> matrices
(A = [-0.3227486121839514;;], B = [1.0;;], C = [1.0;;], D = [-0.0;;])

```

I have a video going through some aspects of linearization of MTK models here

[![](https://global.discourse-cdn.com/julialang/original/3X/e/f/efc5408321efbc8e0ff1414297073a7ec6bac349.jpeg "How to Linearize a ModelingToolkit model") ](https://www.youtube.com/watch?v=-XOux-2XDGI&t=21s)

unfortunately, the font size is really small, so you need a large monitor 😉

---

<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 9, 2023, 4:39pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/3 "2023-03-09T16:39:01Z")

</div>

> [@baggepinnen](#):
>
> > [@BLI](#):
> >
> > (i) is that a way to give _default_ values to the variables in question?
> 
> yes

OK – since linearization is performed on the _symbolic_ model, before parameters are added in the `ODEProblem` function, I guess it is _required_ to specify default values to be able to `linearize` the model?

> [@baggepinnen](#):
>
> > [@BLI](#):
> >
> > is it a prerequisite for the following _linearize_ function that you have an “input” variable that has not been determined?
> 
> yes, but you can make use of [analysis points](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/) to analyze fully connected models as well.

My _interpretation_ of the little I read about an “analysis point” is that it is a way to add an “artificial input” at a given point in the system, so that this added input gives the needed extra undefined input??

> [@baggepinnen](#):
>
> > [@BLI](#):
> >
> > do you have to use the full model when you linearize? Or can you use the model found by `structural_simplify()`?
> 
> You have to use the full model, at least for now.

Makes sense, since I don’t think `structural_simplify()` works unless you have the same number of equations and unknowns.

Hm. I need to have a “balanced”/complete model that can be simulated in order to find the steady states with given constant input. But then I need to change the model and make at least one _variable_ undefined in order to be able to linearize the model.

Unless I can:

- temporarily add an (input) variable at an analysis point in the linearization routine, or
- “un-specify” a variable in the linearization routine.

Anyways – I understand a little bit more now.

---

<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 9, 2023, 4:55pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/4 "2023-03-09T16:55:00Z")

</div>

> [@BLI](#):
>
> I guess it is _required_ to specify default values to be able to `linearize` the model?

You can also pass the keyword argument `op` (for operating point) to `linearize`. From the docstring

```julia
  (; A, B, C, D), simplified_sys = linearize(sys, inputs, outputs; t=0.0, op = Dict(), allow_input_derivatives = false, kwargs...)
  (; A, B, C, D) = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false)

  Return a NamedTuple with the matrices of a linear statespace representation on the form

  \begin{aligned}
ẋ &= Ax + Bu\\
y &= Cx + Du
\end{aligned}

  The first signature automatically calls linearization_function internally, while the second signature expects
  the outputs of linearization_function as input.

  op denotes the operating point around which to linearize. If none is provided, the default values of sys are
  used.

```

> [@BLI](#):
>
> My _interpretation_ of the little I read about an “analysis point” is that it is a way to add an “artificial input” at a given point in the system, so that this added input gives the needed extra undefined input??

An analysis point behaves according to what problem you ask it to solve. If you just run `structural_simplify`, it’s like they weren’t there at all. If you ask for a sensitivity function, then an artificial input is added, if you ask for a loop opening, the connection through the analysis point is broken.

> [@BLI](#):
>
> Makes sense, since I don’t think `structural_simplify()` works unless you have the same number of equations and unknowns.

Correct, thsi is part of why analysis points exist, to make this easier so that you can linearize the same model as you would have called simplification on.

> [@BLI](#):
>
> But then I need to change the model and make at least one _variable_ undefined in order to be able to linearize the model.

Yes, _unless_ you use analysis points, which solve exactly this problem. This is described in the video (at [3:51](https://youtu.be/-XOux-2XDGI?t=231) to be exact), with an example demonstrating exactly this problem and the solution.

---

<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 9, 2023, 5:04pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/5 "2023-03-09T17:04:02Z")

</div>

With a FullHD resolution, 13.5inch screen, I think I have to look at it when I get access to a bigger screen 😉 .

---

<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 10, 2023, 8:13am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/6 "2023-03-10T08:13:22Z")

</div>

> [@baggepinnen](#):
>
> yes, but you can make use of [analysis points](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/) to analyze fully connected models as well.

Hm. Your link points me to the documentation of `ModelingToolkitStandardLibrary`, where it is stated:

“An analysis point can be created explicitly using the constructor [`AnalysisPoint`](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/#ModelingToolkitStandardLibrary.Blocks.AnalysisPoint-Tuple%7BAny,%20Any%7D), or automatically when connecting two causal components using `connect`”

**Question** : in what package is constructor `AnalysisPoint` located? I’m trying to use it, but am told `UndefVarError: AnalysisPoint not defined`.

In my first attempt of “Building component-based, hierarchical models” (not the one in the documentation), I successfully defined connections by, e.g.,

```julia
connections = [tank_U.ṁ_e ~ tank_L.ṁ_i,
                tank_L.ṁ_e ~ tank_B.ṁ_i,
                tank_B.ṁ_e ~ pipe.ṁ,
                pipe.ṁ ~ tank_U.ṁ_i,
                pipe.h_b ~ tank_B.h]

```

I assume I could have used:

```julia
connections = [connect(tank_U.ṁ_e, tank_L.ṁ_i), ...]

```

instead.

And that I could have inserted Analysis Points, e.g., by:

```julia
connections = [connect(tank_U.ṁ_e, :u, tank_L.ṁ_i), ...]

```

**Question** : Is there a difference between setting `a ~ b` and `connect(a,b)`? Would it be possible to do `a ~ :u ~ b`?

---

<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 10, 2023, 8:44am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/7 "2023-03-10T08:44:28Z")

</div>

> [@BLI](#):
>
> Hm. Your link points me to the documentation of `ModelingToolkitStandardLibrary`

Yes, analysis points are defined in the standard library, the first paragraph on the link sums it up

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

The reason is that analysis points are causal in nature, and MTK does not have a concept of causal modeling, only the components in `ModelingToolkitStandardLibrary.Blocks` have this causal nature, and analysis points thus only make sense in `connection`s between Blocks form this standard library, such as `RealInput` and `RealOutput`.

A `connect` statement is translated to one or several equations depending on what connectors are included in the connection. A connection between two electrical pins, for example, implies Kirchoff’s current and voltage laws being introduced by `connect`.

> [@BLI](#):
>
> I assume I could have used:
> 
> ```julia
> connections = [connect(tank_U.ṁ_e, tank_L.ṁ_i), ...]
> 
> ```
> 
> instead.

No, you cannot `connect` variables, only connectors. A connector is a special `ODESystem` defined using the `@connector` annotation. See examples here

> **[Acausal Component-Based Modeling · ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/dev/tutorials/acausal_components/)**
>
> Documentation for ModelingToolkit.jl.

---

<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 11, 2023, 8:11pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/8 "2023-03-11T20:11:49Z")

</div>

The `Analysis Point` idea looks really handy. However, it requires some deeper knowledge of `ModelingToolkit` than I need in my course.

I figured out a way to “sort of” get around it… Don’t know how smart the way is, though. Here is a very simple example…

```julia
# Packages... perhaps I don't need all
using ModelingToolkit
using ModelingToolkit: inputs, outputs
using ModelingToolkitStandardLibrary
using DifferentialEquations
using Plots, LaTeXStrings
using ControlSystems, ControlSystemsMTK
using LinearAlgebra

# External driving function
ṁ_const(t) = 2 # Used for finding steady state, etc.
ṁ_step(t) = t < 1 ? 2 : 1.5 # Could be used to inject a step change

ṁ(t) = ṁ_const(t)  
@register_symbolic ṁ(t)

# Creating a "partial" tank model (a water tank...) where external driving function is excluded
@parameters ρ=1 A=5 K=5 h_ς=3
@variables t m(t)=1.5*ρ*A ṁ_i(t) ṁ_e(t) V(t) h(t)
Dt = Differential(t)

eqs_p = [Dt(m) ~ ṁ_i-ṁ_e,
        m ~ ρ*V,
        V ~ A*h,
        ṁ_e ~ K*sqrt(h/h_ς)]

@named tank_p = ODESystem(eqs_p) # Partial tank model, i.e., without inputs

# Creating an "input" model
eqs_i = [ṁ_i ~ ṁ(t)]
@named tank_i = ODESystem(eqs_i)

# Extending the partial model with the input model, and simplifying it:
tank = extend(tank_i, tank_p)
tank_simp = structural_simplify(tank)

# Quering what is the state; the result is the mass `m`:
states(tank_simp)

# Finding steady state with the given constant input:
tspan = (0.0,1e3)
prob = ODEProblem(tank_simp, [], tspan)
sol = solve(prob)
m_ss = sol[end][1]

# Linearizing the model using MTK
mats, tank_ = linearize(tank_p, [ṁ_i], [h], op=Dict(m=>m_ss, ṁ_i=>ṁ(0)))
sys = ss(mats...)

# Alternatively, linearizing the model using ControlSystemsMTK
named_ss(tank_p, [ṁ_i], [h], op=Dict(m=>m_ss, ṁ_i=>ṁ(0)))

```

OK - I think this approach is useful for _simple_ models; clearly it is not as “smart” as the `Analysis Point` idea.

**Question** : suppose I have computed the eigenvalue _theoretically_/by hand, and have found that it is:

```julia
λ = - (K/√(h_ς))^2/(2*ρ*A*ṁ(0))

```

which is a symbolic expression.

… is there a way to “evaluate” this expression numerically, by using the specified default values for the `parameters`?

---

<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, 2023, 12:47pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/9 "2023-03-12T12:47:58Z")

</div>

> [@BLI](#):
>
> … is there a way to “evaluate” this expression numerically, by using the specified default values for the `parameters`?

```julia
julia> ModelingToolkit.defaults(tank_simp)
Dict{Any, Any} with 5 entries:
  h_ς => 3
  K => 5
  m(t) => 1.5A*ρ
  A => 5
  ρ => 1

julia> λ = - (K/√(h_ς))^2/(2*ρ*A*ṁ(0))
(-((K / sqrt(h_ς))^2)) / (4A*ρ)

julia> substitute(λ, ModelingToolkit.defaults(tank_simp))
-0.4166666666666668

julia> substitute(λ, ModelingToolkit.defaults(tank_simp)) |> typeof
Num

julia> substitute(λ, ModelingToolkit.defaults(tank_simp)) |> Symbolics.value
-0.4166666666666668

julia> substitute(λ, ModelingToolkit.defaults(tank_simp)) |> Symbolics.value |> typeof
Float64

```

You could also use `build_function` to create a function that takes the parameters as inputs and evaluates `λ`

---

<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, 2023, 8:52pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/10 "2023-03-14T20:52:08Z")

</div>

Hm… I create matrices `A, B, C, D` by `calculate_jacobian` in MTK. If I take the result when I have substituted numeric values into the matrices, and pipe it through `Symbolics.value`, the `eltype` of the matrices is still `Num`.

Because of this (?), I doing `ControlSystems.ss(A, B, C, D)` fails.

- Is there a way to convert the `eltype` of these matrices into something that `ss()` will accept?

---

<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 15, 2023, 5:23am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/11 "2023-03-15T05:23:32Z")

</div>

> [@BLI](#):
>
> I create matrices `A, B, C, D` by `calculate_jacobian`

Are you sure you mean `calculate_jacobian`? This function does not do exactly this, and if it did in your case, it’s likely due to luck and not something that can be relied upon always working.

I added a function to compute the `A,B,C,D` matrices with symbolic differentiation

> <https://github.com/SciML/ModelingToolkit.jl/pull/2119>
>
> For simple systems, it's possible to perform linearization symbolically, this ca…n be super handy at times.
> 
> This PR adds the function \`linearize\_symbolic\` that behaves very similar to \`linearize\`, but returns symbolic matrices instead. For now, it has some known limitations, it can not handle \`ifelse\` statements appearing in the dynamics function. This can likely be fixed in the future by either
> \- allowing code gen to output code with \`ifelse\` statements intact rather than inserting branches
> \- Not relying on \`generate\_function\` and instead generalizing the internals of \`generate\_function\` sufficiently to support this use case.

but I’m not sure if the approach I’ve taken is the best, it has some known limitations stated in the PR. It does work for simple systems though, see the tests added in the PR.

---

<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 15, 2023, 8:35am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/12 "2023-03-15T08:35:14Z")

</div>

> [@baggepinnen](#):
>
> Are you sure you mean `calculate_jacobian`?

Yes, `calculate_jacobian`. Of course, I added quite a few other lines to deconstruct the Jacobian and find matrices `A`, `B`, `C`, `D`.

However, currently, my function is “clumsy” and based on some restrictive assumptions. And – since I’m not 100% sure as to how `compute_jacobian` treats the equations I have (does it manipulate them? Does it change order? Etc.), I may switch to the Symbolics `jacobian` algorithm, which takes a vector of terms and a vector of variables as inputs, and thus is more transparent about what it does.

---

<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 15, 2023, 9:48am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/13 "2023-03-15T09:48:03Z")

</div>

I reiterate my advice to not rely on `calculate_jacobian` for this, even with additional post processing. It is not hard to construct an example system on which it fails and does not do what you expect. For example, even simple component-based systems with `connect` statements is likely to produce garbage Jacobians with this approach. Also, you can not guarantee that the signal you consider an output `y = Cx + Du` will actually be present in the equations that are used to compute the jacobian, it might have been moved to the observed variables, furthermore, the order of the state variables is not guaranteed to be stable between MTK versions, neither is the selection of which variables are included in the state. This is why I submitted the PR linked above, it ensures that you get correct input and output jacobians for the specified inputs and outputs.

To summarize, `calculate_jacobian` does not know anything about inputs or outputs, it only knows about variables. Variables that you consider inputs and outputs may thus be removed from the equations of the system and placed in the observed equations (even inputs), and you may thus fail to compute any of `B,C,D`.

The `linearize_symbolic` function in the PR performs the required `io_preprocessing` in order to properly mark and handle inputs and outputs.

---

<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 15, 2023, 10:45am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/14 "2023-03-15T10:45:18Z")

</div>

Ah. I’m not as ambitious as this… I use a “partial” model (underdetermined), i.e., a model where the number of unknown dependent variables [`n_v`] is larger than the number of equations [`n_eq`], and I assume that the number of inputs `n_u` must satisfy `n_u = n_v - n_eq`. Also, I assume that the resulting DAE is index 1 (which is a strong limitation, but probably ok for my use).

Essentially, I do the following:

1. I create a model `mod` from the equations `eq` as `mod = ODESystem(eq)`.
2. I extract unknown dependent variables by `states(mod)`; I need this to extract matrices.
3. I create the Jacobian by `calculate_jacobian(mod)`.

My _uncertainty_ with `calculate_jacobian` is the following:

A. Are the _states_ I find by `states(mod)` equal to the dependent variables in the original equations `eq`? Or does the constructor `ODESystem` _in itself_ manipulate the model/extract _observed_ variables, etc.?

- NOTE: I do not issue any `structural_simplify` operation or any other manipulation on the model, but it is possible that the constructor `ODESystem` does something in itself.

B. Are the _equations_ found by `equations(mod)` identical to the original equations `eq`?

- Again, note that I have not issued any method to manipulate the system, but it is possible that the constructor `ODESystem` does this on its own.

If the answer to either A or B is _NO_, then I can not use `calculate_jacobian` _for my use_.

I can get around this by:  
I. Directly using `eq` instead of `equations(mod)`; that is no problem.  
II. Find a way to automate which _variables_ are involved in the equations `eq` (which `states(mod)` kind of does – unless `ODESystem` changes the model, strips off `observed`, etc…).

So far, I have been hoping that since, in my way of doing it, I have an underdetermined system (more unknown variables than equations), `ODESystem` doesn’t do much… I may be wrong.

Anyways, I can get around using MTK for this particular problem if I can automate finding the variables in `eq`.

---

<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 15, 2023, 10:49am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/15 "2023-03-15T10:49:00Z")

</div>

What’s the reason for not using `linearize_symbolic` from the PR?

---

<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 15, 2023, 10:57am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/16 "2023-03-15T10:57:59Z")

</div>

Two reasons:

1. I want to learn something and understand how things work,
2. I have barely ever used Git, so I don’t even know how to install anything “from the PR”.

---

<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 15, 2023, 1:44pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/17 "2023-03-15T13:44:48Z")

</div>

> [@BLI](#):
>
> Are the _states_ I find by `states(mod)` equal to the dependent variables in the original equations `eq`?

The `states(mod)` will include also input variables, so not only dependent variables.

> [@BLI](#):
>
> Are the _equations_ found by `equations(mod)` identical to the original equations `eq`?

based on the implementaiton of `equations`, it looks like they are identical, except that they are namespaced appropriately for composed systems.

```julia
function equations(sys::AbstractSystem)
    eqs = get_eqs(sys)
    systems = get_systems(sys)
    if isempty(systems)
        return eqs
    else
        eqs = Equation[eqs;
                       reduce(vcat,
                              namespace_equations.(get_systems(sys));
                              init = Equation[])]
        return eqs
    end
end

```

---

<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 15, 2023, 1:49pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/18 "2023-03-15T13:49:48Z")

</div>

Here’s how you can use the new `linearize_symbolic`

```julia
(@v1.8) pkg> activate --temp

(jl_n6mpyc) pkg> add ModelingToolkit#master

```

```julia
julia> using ModelingToolkit

julia> @variables t m(t) ṁ_i(t) ṁ_e(t) V(t) h(t);

julia> @parameters ρ=1 A=1 K=1 h_ς=1;

julia> Dt = Differential(t);

julia> @variables ṁinp(t)=1;

julia> eqs = [Dt(m) ~ ṁ_i-ṁ_e,
               m ~ ρ*V,
               V ~ A*h,
               ṁ_e ~ K*sqrt(h/h_ς),
               ṁ_i ~ ṁinp];

julia> @named tank = ODESystem(eqs);

julia> matrices, _ = ModelingToolkit.linearize_symbolic(tank, [ṁinp], [h])
((A = Num[(-K) / (2A*h_ς*ρ*sqrt(m(t) / (A*h_ς*ρ)));;], B = Num[1;;], C = Num[1 / (A*ρ);;], D = Num[0;;]), ODESystem(0x0000000000000000, Equation[Differential(t)(m(t)) ~ ṁ_i(t) - ṁ_e(t)], t, Any[m(t)], Any[ρ, A, h_ς, K, ṁinp(t)], nothing, Dict{Any, Any}(:m => m(t), :ρ => ρ, :A => A, :h_ς => h_ς, :K => K, :V => V(t), :h => h(t), :ṁinp => ṁinp(t), :ṁ_i => ṁ_i(t), :ṁ_e => ṁ_e(t)…), Any[], Equation[ṁ_i(t) ~ ṁinp(t), V(t) ~ m(t) / ρ, h(t) ~ V(t) / A, ṁ_e(t) ~ K*sqrt(h(t) / h_ς)], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :tank, ODESystem[], Dict{Any, Any}(ṁinp(t) => 1, h_ς => 1, K => 1, A => 1, ρ => 1), nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[ModelingToolkit.SymbolicContinuousCallback(Equation[], Equation[])], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, TearingState of ODESystem, ModelingToolkit.Substitutions(Equation[ṁ_i(t) ~ ṁinp(t), V(t) ~ m(t) / ρ, h(t) ~ V(t) / A, ṁ_e(t) ~ K*sqrt(h(t) / h_ς)], [Int64[], Int64[], [2], [3, 2]], nothing), false, nothing, nothing))

julia> matrices.A
1×1 Matrix{Num}:
 (-K) / (2A*h_ς*ρ*sqrt(m(t) / (A*h_ς*ρ)))

```

---

<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 15, 2023, 5:32pm UTC](https://discourse.julialang.org/t/mtk-linearization/95803/19 "2023-03-15T17:32:05Z")

</div>

OK – here is an “academic” example with my “linearize” function `linearize_partial_model`.

- I assume that the model is an Index 1 DAE, so not for any problem
- Base the algorithm on a “partial”/underdetermined model
- “Trusts” that `states(mod)` and `equations(mod)` have not changed the variables and equations through the `ODESystem` constructor.

```julia
# Packages -- not sure if everything is used
using ModelingToolkit
using ModelingToolkit: inputs, outputs
using DifferentialEquations
using Plots, LaTeXStrings
using ControlSystems, ControlSystemsMTK
using LinearAlgebra
using Symbolics

# My linearization function
function linearize_partial_model(mod_p, u_mod, y_mod; x_opt = [])
    #
    # This function assumes that `mod_p` is an index 1 DAE 
    #
    # mod_p: ModelingToolkit partial model, i.e., with inputs undefined
    # u_mod: vector of variables which are inputs in `states(mod_p)`
    # y_mod: vector of variables in `state(mod_p)` which are outputs (inputs not allowed)
    # x_opt: optional state vector ξ if other than differential variables x
    #
    # Unknown dependent variables v 
    v_mod = states(mod_p)
    # Finding vector of states and vector of equations
    
    eqs = equations(mod_p)
    n_v = length(v_mod)
    n_eqs = length(eqs)
    n_y = length(y_mod)
    #
    # Finding indices of differential and algebraic equations
    # -- extracting lhs and rhs of equations
    eqs_lhs = [eq.lhs for eq in eqs]
    eqs_rhs = [eq.rhs for eq in eqs]
    # -- vector of differential of unknown variables
    dvdt_mod = Differential(t).(v_mod)
    # -- assuming lhs of equations has only one occurrence of
    # -- differential of a given variable
    id = [findfirst(isequal.(eq,dvdt_mod)) for eq in eqs_lhs]
    # -- finding indices of differential and algebraic equations
    idxs_eq_d = id[id .!== nothing]
    idxs_eq_a = setdiff(1:n_eqs,idxs_eq_d)
    #
    # Finding vector field for DAE
    eq_field = [eqs_rhs[idxs_eq_d]; eqs_rhs[idxs_eq_a].-eqs_lhs[idxs_eq_a]]
    # Finding Jacobian
    J = Symbolics.jacobian(eq_field,v_mod)
    #
    # Finding indices for differential variables, inputs,
    # algebraic variables, and outputs
    # -- differential variables
    id = [findfirst(isequal.(Differential(t)(var),eqs_lhs[idxs_eq_d])) for var in v_mod]
    idxs_x = id[id .!== nothing]
    # -- inputs
    idxs_u = [findmax(isequal.(u,v_mod))[2] for u in u_mod]
    # -- finding algebraic variables
    idxs_z = setdiff(1:n_v,idxs_x,idxs_u)
    # -- outputs
    idxs_y = [findmax(isequal.(y,v_mod))[2] for y in y_mod]
    # -- checking that outputs are not inputs
    if !isdisjoint(idxs_y,idxs_u)
        println("Outputs can not be inputs.")
        println("Outputs: $(y_mod)")
        println("Inputs: $(u_mod)")
    end
    #
    # Deconstructing Jacobian
    J_d_x = J[idxs_eq_d,idxs_x]
    J_d_z = J[idxs_eq_d,idxs_z]
    J_d_u = J[idxs_eq_d,idxs_u]
    J_a_x = J[idxs_eq_a,idxs_x]
    J_a_z = J[idxs_eq_a,idxs_z]
    J_a_u = J[idxs_eq_a,idxs_u]
    #
    # LTI assuming states are differential variables
    # -- A and B matrices
    AA = simplify.(J_d_x - J_d_z*(J_a_z\J_a_x))
    BB = simplify.(J_d_u - J_d_z*(J_a_z\J_a_u))
    # --Jacobian matrix for outputs; contains ones and zeros
    J_y = zeros(n_y,n_v)
    for i in 1:n_y
        J_y[i,:] = isequal.(y_mod[i],v_mod)
    end
    J_y_x = J_y[:,idxs_x]
    J_y_z = J_y[:,idxs_z]
    J_y_u = J_y[:,idxs_u]
    # -- C and D matrices
    CC = simplify.(J_y_x - J_y_z*(J_a_z\J_a_x))
    DD = simplify.(J_y_u - J_y_z*(J_a_z\J_a_u))
    #
    # Checking whether optional states are requested
    n_x = length(idxs_x)
    n_ξ = length(x_opt)
    #
    if n_ξ == 0
        return AA, BB, CC, DD, v_mod[idxs_x]
    elseif n_ξ != n_x
        println("Incorrect size of optional state")
        return AA, BB, CC, DD, v_mod[idxs_x]
    else
        # Finding transformation matrix TT from original to new state
        # -- indices for new state
        idxs_ξ = [findmax(isequal.(x,v_mod))[2] for x in x_opt]
        # Jacobian matrix for new state; contains ones and zeros
        J_ξ = zeros(n_ξ,n_v)
        for i in 1:n_ξ
            J_ξ[i,:] = isequal.(x_opt[i],v_mod)
        end
        J_ξ_x = J_ξ[:,idxs_x]
        J_ξ_z = J_ξ[:,idxs_z]
        J_ξ_u = J_ξ[:,idxs_u]
        #
        TT = simplify.(J_ξ_x - J_ξ_z*(J_a_z\J_a_x))
        TT_u = simplify.(J_ξ_u - J_ξ_z*(J_a_z\J_a_u))
        # -- not allowed to have inputs in ξ, i.e., must require TT_u ≡ 0
        if sum(TT_u) != 0
            println("Optional states lead to input derivatives")
            return AA, BB, CC, DD, v_mod[idxs_x]
        end
        AA_ = simplify.(TT*AA*inv(TT))
        BB_ = simplify.(TT*BB)
        CC_ = simplify.(CC*inv(TT))
        DD_ = DD
        #
        return AA_, BB_, CC_, DD_, x_opt
    end
end

# External forcing functions
u_ṁ_step(t) = t<2 ? 1 : 1.5
u_u_step(t) = t<15 ? 1 : 0.8
u_T_i_step(t) = t<10 ? 288.15 : 293.15
u_Q̇_step(t) = t<5 ? 0 : 1e5

u_ṁ_const(t) = 1
u_u_const(t) = 1
u_T_i_const(t) = 288.15
u_Q̇_const(t) = 0

# Parameters, variables, registered functions
@parameters ρ=1e3 A=1.5e-2 K_e=2 h_ς=3e-1 Ĥ_o=0 ĉ_p=4.2e3 T_o=298.15 p_o=1.01e5 p_a=1.01e5 g=9.81

@variables t h(t)=1.5e-1 m(t)=ρ*h*A ṁ_i(t) ṁ_e(t) p(t)=1.01e5 T(t)=293.15 
@variables Ĥ(t)=Ĥ_o + ĉ_p*(T - T_o) + (p - p_o)/ρ
@variables U(t)= m*Ĥ - p*m/ρ Ḣ_i(t) Ḣ_e(t) Ẇ_v(t) Q̇(t)
@variables V(t) u(t) 
@variables Ĥ_i(t) T_i(t) Q̇_(t)

D = Differential(t)

@register_symbolic u_ṁ(t) 
@register_symbolic u_u(t)
@register_symbolic u_T_i(t)
@register_symbolic u_Q̇_(t)

# DAE Model: partial for linearization, and complete for simulation
# DAE Formulation
eqs_p = [D(m) ~ ṁ_i - ṁ_e,
        D(U) ~ Ḣ_i - Ḣ_e + Ẇ_v + Q̇,
        m ~ ρ*V,
        V ~ A*h,
        ṁ_e ~ K_e*u*sqrt(h/h_ς),
        U ~ m*Ĥ - p*V,
        Ĥ ~ Ĥ_o + ĉ_p*(T-T_o) + (p-p_o)/ρ,
        Ĥ_i ~ Ĥ_o + ĉ_p*(T_i-T_o) + (p_a-p_o)/ρ,
        Ḣ_i ~ ṁ_i*Ĥ_i,
        Ḣ_e ~ ṁ_e*Ĥ,
        Ẇ_v ~ -p*(ṁ_i - ṁ_e)/ρ,
        Q̇ ~ Q̇_*h,
        p ~ p_a + ρ*g*h]

@named sys_p = ODESystem(eqs_p)

eqs_i = [ṁ_i ~ u_ṁ(t), u ~ u_u(t), T_i ~ u_T_i(t), Q̇_ ~ u_Q̇_(t)]
@named sys_i = ODESystem(eqs_i)

sys = extend(sys_i, sys_p)

sys_simp = structural_simplify(sys)

states(sys_simp)

# ODE formulation of system
# ODE Formulation
eqs_ode_p = [D(h) ~ (ṁ_i - K_e*u*sqrt(h/h_ς))/(ρ*A),
            D(T) ~ (ṁ_i*(ĉ_p*(T_i - T) - g*h) + Q̇_*h)/(ρ*A*h*ĉ_p)]

@named sys_ode_p = ODESystem(eqs_ode_p)

sys_ode = extend(sys_i, sys_ode_p)

sys_ode_simp = structural_simplify(sys_ode)

states(sys_ode_simp)

# "Partial" equations -- these are assumed to be the same as `eqs_p`:
equations(sys_p)

# "States" of partial model -- these are assumed to be the same as the variables used in `eqs_p`
v_sys = states(sys_p)

# Specifying vectors of inputs and outputs for the linearization
u_sys = [ṁ_i,u,T_i,Q̇_]
y_sys = [h,T]

# Linearization -- uses differential variables (here: `[m, U]`) as states in result
AA, BB, CC, DD, state = linearize_partial_model(sys_p, u_sys,y_sys)
state

# Linearization -- with specification of states to be `x_opt = [h,T]` instead of default
AA, BB, CC, DD, state = linearize_partial_model(sys_p, u_sys,y_sys; x_opt=[h,T])
state

# Resulting matrices
AA
BB
CC
DD

# Running numeric model (complete model) to find steady states
# -- specifying driving functions

u_ṁ(t) = u_ṁ_const(t)
u_u(t) = u_u_const(t)
u_T_i(t) = u_T_i_const(t)
u_Q̇_(t) = u_Q̇_const(t)

# -- finding steady states
tspan = (0.0,1e3)

prob = ODEProblem(sys_simp, [], tspan)
sol = solve(prob)
name_s = states(sys)
s_ss = sol(tspan[2]; idxs = name_s)
m_ss = sol(tspan[2]; idxs = m)
U_ss = sol(tspan[2]; idxs = U)
Ĥ_ss = sol(tspan[2]; idxs = Ĥ)
h_ss = sol(tspan[2]; idxs = h)
T_ss = sol(tspan[2]; idxs = T)

# -- evaluating numerical values of AA, BB, CC, DD
default_vals = ModelingToolkit.defaults(sys_p)

linsys_symb_A = substitute(AA, Dict(name_s .=> s_ss)) |> x -> substitute(x, default_vals) |> Symbolics.value

linsys_symb_B = substitute(BB, Dict(name_s .=> s_ss)) |> x -> substitute(x, default_vals)

linsys_symb_C = substitute(CC, Dict(name_s .=> s_ss)) |> x -> substitute(x, default_vals)

linsys_symb_D = substitute(DD, Dict(name_s .=> s_ss)) |> x -> substitute(x, default_vals)

# Linearizing the model using MTK
sys_in = [ṁ_i, u, T_i, Q̇_]
sys_out = [h, T]
op_0 = Dict(m=>m_ss, U=>U_ss, Ĥ=>Ĥ_ss, ṁ_i=>u_ṁ(0), u=>u_u(0), T_i=>u_T_i(0), Q̇_=>u_Q̇_(0))
op_ode_0 = Dict(h=>h_ss, T=>T_ss, ṁ_i=>u_ṁ(0), u=>u_u(0), T_i=>u_T_i(0), Q̇_=>u_Q̇_(0))

mats, sys_ = linearize(sys_p, sys_in, sys_out, op=op_0)
linsys_mtk = ss(mats...)

# Linearizing the model with ControlSystemsMTK
linsys_csmtk = named_ss(sys_p, sys_in, sys_out; op=op_0)

```

**Some comments** :

1. MTK finds a linear model with 3 states based on `sys_p`: `m`, `U`, and `Ĥ`.
2. ControlSystemsMTK finds exactly the same 3 state model as MTK in this case.
3. Both MTK and ControlSystemsMTK finds the same 2 state model for the `sys_ode_p` model, i.e., the ODE formulation.
4. Function `linearize_partial_model` finds 2 states based on `sys_p` both when the default states are used (= differential variables `[m, U]`), and when optional states `[h, T]` are given; the models are symbolic. When substituting in numerical values, the version with state `[h, T]` is identical to the one found from MTK and ControlSystemsMTK based on `sys_ode_p`.
5. Function `linearize_partial_model` finds 2 states based on `sys_ode_p`. When the default states are used (`[h, T]`), the symbolic model appears to be the same as the one found based on the DAE `sys_p` with optional states (`[h, T]`) [there are some terms that have not been simplified away, so they _appear_ to be slightly different].
6. The numeric evaluation of the matrices found from `linearize_partial_model` are of type `Num`, and I have not been able to insert them into ControlSystem’s `ss` function.
7. Function `linearize_partial_model` is extremely poorly tested, so I’m sure there are some bugs in it. Also, it is relatively limited in scope.
8. I’m sure `linearize_partial_model` can be made to run more efficiently…
9. I have assumed that `ODESystem()` does not change the order of equations, remove equations, or remove states. If this assumption is false, some minor modifications must be introduced.
10. A slight surprise when it comes to symbolic linear algebra… `A\B` works, but `B/A` does not work. In the latter case, I had to do `B*inv(A)`.

---

<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 20, 2023, 11:52am UTC](https://discourse.julialang.org/t/mtk-linearization/95803/20 "2023-03-20T11:52:43Z")

</div>

I tried to do

```julia
(@v1.8) pkg> activate --temp

(jl_n6mpyc) pkg> add ModelingToolkit#master

```

in my terminal window. I run my code in a Jupyter notebook of VSCode. Does the “temp” activation work when I switch to run it in the Jupyter notebook?

- I’m told that `linearize_symbolic` doesn`t exist.

Next, I issued commands

```julia
using Pkg
Pkg.activate --temp

```

in the Jupyter notebook, but then I’m told that `--` is an invalid operator.

Seems like you rund it in the command window?

[Next page](https://discourse.julialang.org/t/mtk-linearization/95803.md?page=2)