MTK & Linearization

[BLI]

MTK documentation case

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

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:

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…

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

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

fails miserably.

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

[baggepinnen]

@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 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

@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

How to Linearize a ModelingToolkit model

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

[BLI]

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

[baggepinnen]

@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

  (; 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 to be exact), with an example demonstrating exactly this problem and the solution.

[BLI]

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

[BLI]

@baggepinnen

yes, but you can make use of analysis points 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, 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.,

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:

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

instead.

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

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?

[baggepinnen]

@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

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 connections 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:

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

[BLI]

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…

# 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:

λ = - (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?

[baggepinnen]

@BLI

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

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 λ

[BLI]

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?

[baggepinnen]

@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

add symbolic linearization function

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.

[BLI]

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

[baggepinnen]

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.

[BLI]

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.

[baggepinnen]

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

[BLI]

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

[baggepinnen]

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

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

[baggepinnen]

Here’s how you can use the new linearize_symbolic

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

(jl_n6mpyc) pkg> add ModelingToolkit#master

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_ς*ρ)))

[BLI]

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.

# 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).

[BLI]

I tried to do

(@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

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?