@mtkcompile / System(...): unknowns and required u0 change when I rename a component

Specific Domains Modelling & Simulations
,
1

Hi everyone,

I am working with a fairly large ModelingToolkit model built from components. I compile it with something like:

@mtkcompile sys = System(assembly, t; systems = elements)

where elements is the list of the different components previously defined with @named, and assembly contains the connect() calls used to build the full system.

Then I provide u0 and guesses (as symbolic mappings) to build an ODEProblem and solve it.

I am observing a confusing behavior regarding the system unknowns.

In my elements list, I have a component named θ_dynamic. If I rename that component to θ (changing both the @named declaration and the corresponding entry in the elements list), the set of unknowns in sys changes, and therefore the variables I need to provide in u0 / guesses also change.

For example, with the name θ_dynamic, I am asked for something like R1.q (the flow through the R1 resistance). But if I only rename that component to θ, I am instead asked for R2.q.

:confused: What confuses me is:

  • I expected u0 to refer to the actual state variables, so I did not expect it to change just because I renamed one component.
  • More generally, I would like the initialization to be consistent and not depend on component names or on the order of components in elements.

I also noticed another strange behavior: I have the “same” model both in a notebook and in a script, but I do not obtain exactly the same system of equations or the same unknowns.

:red_question_mark: So my questions are:

  1. Is it expected that renaming a subsystem can change the compiled unknowns?
  2. Could this change be related to equation sorting, tearing, or alias elimination during mtkcompile?
  3. What is the best practice to make initialization robust so that it does not depend on unstable unknown ordering?

Thank you for your help !

2

When you say the set of unknowns in sys changes, do you mean the number of unknowns, or just the order? If the number of unknowns is changing, this means the number of equations is changing and therefore the tearing and equation simplification would be somehow different even though all that changed was a name. That definitely should not be happening.

Regarding a different results from the same model in a script and notebook, I would check that the version of MTK is the same in both environments. MTK is changing rapidly and it’s very possible you have a different versions loaded.

3

Thank you for your reply!

When I say that the set of unknowns in sys changes, I do not mean that the number of unknowns changes — it remains the same. However, one unknown is replaced by another. For example, with the original name θ_dynamic, I am asked for R1.q, but if I rename the component to θ, the system instead asks for R2.q.

So the number of unknowns stays identical, but which variable appears as an unknown changes when I rename the θ element.

For your suggestion on the MTK version, I checked it and both environments are using the same ModelingToolkit version, and they are running in the same Julia environment.

4

Can you provide the full list of unknowns for both cases? If it’s simply a change of order, providing a symbolic map for u0 should work regardless of the order. If you can provide a full MWE that would be really helpful, but at least the full list of unknowns for both cases would be very useful.

6

Here is a MWE with the change of name that make a change in the unknows

using ModelingToolkit
using DifferentialEquations
using SciCompDSL
using ModelingToolkit: t_nounits as t, D_nounits as D
using SciMLStructures: Tunable, canonicalize, replace
using PreallocationTools
using Plots

@connector Pin begin
  p(t)
  q(t), [connect = Flow]
end

@mtkmodel OnePort begin
  @components begin
      in = Pin()
      out = Pin()
  end
  @variables begin
      Δp(t)
      q(t)
  end
  @equations begin
      Δp ~ in.p - out.p
      0 ~ in.q + out.q
      q ~ in.q
  end
end

@mtkmodel ConstantPressure begin
  @components begin
      node = Pin()
  end
  @parameters begin
      P
  end
  @equations begin
      node.p ~ P
  end
end


@mtkmodel ForcedFlow begin
    @components begin
        in = Pin()
        out = Pin()
    end

    @parameters begin
        Q
    end

    @equations begin
        out.q ~ - Q 
        in.q ~ Q
    end
end

@mtkmodel Resistor begin
  @extend Δp, q = oneport = OnePort()
  @parameters begin
      R
  end
  @equations begin
      Δp ~ q * R
  end
end

@mtkmodel ResistorValve begin
    @extend Δp, q = oneport = OnePort()

    @parameters begin
        R
    end

    @equations begin
        q ~ Δp / R * (Δp > 0) 
    end
end

@mtkmodel Capacitor_1 begin
    @components begin
        in = Pin()
        out = Pin()
    end

    @parameters begin 
        C
    end

    @variables begin 
        V(t)
    end

    @equations begin 
        out.q ~ - in.q 
        D(V) ~ in.q 
        in.q ~ (D(in.p)-D(out.p)) * C
    end
end

@mtkmodel Capacitor_2 begin
    @components begin
        in = Pin()
        out = Pin()
    end

    @parameters begin 
        k
        p0
    end

    @variables begin 
        V(t)
        C(t)
    end

    @equations begin 
        out.q ~ - in.q 
        D(V) ~ in.q 
        in.q * (in.p - out.p - p0) * k ~ (D(in.p)-D(out.p))
        C ~ 1/(k*(in.p - out.p -p0))
    end
end

# -------------------------
# Numerical values
vals = Dict{Symbol, Float64}()
vals[:Pa_v] = 100.0
vals[:Psinus_v] = 6.34
vals[:Pic_v] = 9.5
vals[:P0_v] = 3.2
vals[:Prod_v] = 0.006
vals[:Ra_v] = 6.0
vals[:Rc_v] = 0.08
vals[:Rv_v] = 0.77
vals[:Ro_v] = 525.0
vals[:Ca_v] = 0.02
vals[:E_v] = 0.35*1.3

# Define elements
@named Pa = ConstantPressure(P = vals[:Pa_v]);
@named Psinus = ConstantPressure(P = vals[:Psinus_v]);
@named Pic = ConstantPressure(P = vals[:Pic_v]);
@named Prod = ForcedFlow(Q = vals[:Prod_v]);
@named Ra = Resistor(R = vals[:Ra_v]);
@named Rc = Resistor(R = vals[:Rc_v]);
@named Rv = Resistor(R = vals[:Rv_v]);
@named Ro = ResistorValve(R = vals[:Ro_v]);

@named Ca = Capacitor_1(C = vals[:Ca_v]);
@named Cv = Capacitor_2(k = vals[:E_v], p0 = vals[:P0_v]);

# Define the elements used in the steady-state model
elements_steady_state = [Pa, Psinus, Pic, Prod, Ra, Rc, Rv, Ro, Ca, Cv];

# Describe how components are connected
assembly_steady_state = [
    connect(Pa.node, Ra.in, Ca.in);
    connect(Ra.out, Rc.in, Prod.in);
    connect(Rc.out, Cv.in, Rv.in);
    connect(Ca.out, Prod.out, Cv.out, Ro.in, Pic.node);
    connect(Rv.out, Ro.out, Psinus.node);
];

# Declare the system as an ODESystem
@mtkcompile sys_steady_state = System(assembly_steady_state,t; systems = elements_steady_state)

We have
Model sys_steady_state:
Equations (4):
4 standard: see equations(sys_steady_state)
Unknowns (4): see unknowns(sys_steady_state)
Cv₊in₊q(t)
Cv₊V(t)
Cv₊in₊p(t)
Ca₊V(t)
Parameters (11): see parameters(sys_steady_state)
Pa₊P
Psinus₊P
Pic₊P
Prod₊Q
⋮
Observed (46): see observed(sys_steady_state)

Then when doing this

# Define elements
@named Pa_steady_state = ConstantPressure(P = vals[:Pa_v]);
@named Psinus = ConstantPressure(P = vals[:Psinus_v]);
@named Pic = ConstantPressure(P = vals[:Pic_v]);
@named Prod = ForcedFlow(Q = vals[:Prod_v]);
@named Ra = Resistor(R = vals[:Ra_v]);
@named Rc = Resistor(R = vals[:Rc_v]);
@named Rv = Resistor(R = vals[:Rv_v]);
@named Ro = ResistorValve(R = vals[:Ro_v]);
@named Ca = Capacitor_1(C = vals[:Ca_v]);
@named Cv = Capacitor_2(k = vals[:E_v], p0 = vals[:P0_v]);

# Define the elements used in the steady-state model
elements_steady_state = [Pa_steady_state, Psinus, Pic, Prod, Ra, Rc, Rv, Ro, Ca, Cv];

# Describe how components are connected
assembly_steady_state = [
    connect(Pa_steady_state.node, Ra.in, Ca.in);
    connect(Ra.out, Rc.in, Prod.in);
    connect(Rc.out, Cv.in, Rv.in);
    connect(Ca.out, Prod.out, Cv.out, Ro.in, Pic.node);
    connect(Rv.out, Ro.out, Psinus.node);
];

# Declare the system as an ODESystem
@mtkcompile sys_steady_state = System(assembly_steady_state,t; systems = elements_steady_state)

I have this :
Model sys_steady_state:
Equations (4):
4 standard: see equations(sys_steady_state)
Unknowns (4): see unknowns(sys_steady_state)
Rc₊q(t)
Cv₊V(t)
Cv₊in₊p(t)
Ca₊V(t)
Parameters (11): see parameters(sys_steady_state)
Pa_steady_state₊P
Psinus₊P
Pic₊P
Prod₊Q
⋮
Observed (46): see observed(sys_steady_state)

So I have a change from Cv.in.q(t) to Rc.q(t) in the unknows.