Inconsistent unknown selection in MTK with identical systems (@mtkcompile)

Specific Domains Modelling & Simulations
1

I encountered an inconsistent behavior with ModelingToolkit.jl regarding the selected unknowns after structural simplification / compilation.

I built a minimal working example (MWE) where simply changing the name of one component changes the resulting unknowns of the compiled system, even though:

  • the equations are identical,
  • the assembly is identical,
  • the Julia version is the same,
  • and the package versions are the same.

For example, in the first case:

@named Pa = ConstantPressure(P = vals[:Pa_v]);

the compiled system gives:

Unknowns:
Cv₊in₊q(t)
Cv₊V(t)
Cv₊in₊p(t)
Ca₊V(t)

But when I only rename the component:

@named Pa_steady_state = ConstantPressure(P = vals[:Pa_v]);

the system becomes:

Unknowns:
Rc₊q(t)
Cv₊V(t)
Cv₊in₊p(t)
Ca₊V(t)

So Cv₊in₊q(t) is replaced by Rc₊q(t) even though the physical model and assembly are unchanged.

What is even more surprising is that the selected unknowns / guesses can also differ across two different computers that have:

  • the same Julia version,
  • the same package versions,
  • and supposedly the same environment.

This makes the behavior appear non-deterministic or at least dependent on internal ordering / hashing.

Is this expected behavior from MTK structural simplification / alias elimination, or could this indicate an issue in the compilation pipeline?

MWE below:

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)


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

There seem to have been a lot of changes with how state variables are selected over the past few updates of MTK, from circa 11.18, when I started, to 11.26.

I have been staying up to date with the latest releases and often had to change which states to put in my initial guesses even though everything else remained the same.

Have you seen if the problem is fixed in very latest version? 11.26.5 as of today?