Overdetermined ODEProblem after upgrading from MTK v8.75.0 to MTK v9.72.0

Specific Domains Modelling & Simulations
1

Hi all :waving_hand:

I’m currently upgrading a model from ModelingToolkit v8.75.0 (which worked fine), to v9.72.0.
I’m now running into an issue where my ODEProblem fails to initialize, even though my model contains capacitors (with dynamic equations) and valid initial conditions.

In MTK v8, I used to define a model using @component function style with internal ODESystem(...), and everything worked fine: simulation ran, voltage evolved, etc.

I switched to using @mtkmodel for components (Pin, Resistor, Capacitor, ResistorValve, etc.), and I’m composing the system like this:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations, Plots
@connector Pin begin
    v(t)
    i(t), [connect = Flow]
  end
@mtkmodel OnePort begin
    @components begin
        in = Pin()
        out = Pin()
    end
    @variables begin
        Δv(t)
        i(t)
    end
    @equations begin
        Δv ~ in.v - out.v
        0 ~ in.i + out.i
        i ~ in.i
    end
  end
 @mtkmodel Resistor begin
    @extend Δv, i = oneport = OnePort()
    @parameters begin
        R, [description = "Resistance"]
    end
    @equations begin
        Δv ~ i * R
    end
end
@mtkmodel Capacitor begin
    @extend Δv, i = oneport = OnePort()
    @parameters begin
        C, [description = "Capacitor"]
    end
    @equations begin
        D(Δv) ~ i / C
    end
end
@mtkmodel ResistorValve begin
    @extend Δv, i = oneport = OnePort()
    @parameters begin
        R, [description = "Resistance"]
    end
    @equations begin
        i ~ Δv / R * (Δv > 0)
    end
end
@mtkmodel ConstantVoltage begin
    @components begin
        node = Pin()
    end
    @parameters begin
        V, [description = "Voltage"]
    end
    @equations begin
        node.v ~ V
    end
end
params = Dict(
    # Resistance values
    :R1_v => 7980,
    :R3_v => 811.3,
    :R2_v => 119.7,
    :D1_v => 34e5,
    :D2_v => 698250,
    # Capacitor values
    :C1_v => 15e-6, # Aretrial compliance
    :C2_v => 34e-6, # Venous compliance
    # Voltage values
    :ext_v => 12635,
    :out_v => 8512,
    :in_v => 133000,
    :V1=> 33250,
    :V2=> 18620
)
@named R1 = Resistor(R=params[:R1_v])
@named C1 = Capacitor(C=params[:C1_v])
@named D1 = ResistorValve(R=params[:D1_v])
@named R2 = ResistorValve(R=params[:R2_v])
@named C2 = Capacitor(C=params[:C2_v])
@named R3 = Resistor(R=params[:R3_v])
@named D2 = Resistor(R=params[:D2_v])
@named in = ConstantVoltage(V=params[:in_v])
@named out = ConstantVoltage(V=params[:out_v])
@named ext = ConstantVoltage(V=params[:ext_v])
elements= [in, out, ext, R1, C1, D1, R2, C2, R3, D2]
assembly = [
    connect(in.node, R1.in, C1.in),
    connect(R1.out, R2.in, D1.in),
    connect(R2.out, C2.in, R3.in),
    connect(C1.out, D1.out, C2.out, D2.in, ext.node),
    connect(R3.out, D2.out, out.node),
];

Then I do:


@named _model = ODESystem(assembly, t)
@named model = compose(_model, elements)

sys = structural_simplify(model)
# → returns dC2.Δv(t)/dt = C2.i(t)/C2.C
# and      0 = C2.i(t) + R3.i(t) - R2.i(t)

I expose C2.Δv and R2.Δv to form the initial conditions:

u0 = [
    C2.Δv => 7209,
    R2.Δv => 1671 # values taken from the model on MTKv8 (steady state)
]

And then:

prob = ODEProblem(sys, u0, (0, 10))

When I try to solve the problem, I get this warning and failure:

┌ Warning: Initialization system is overdetermined. 1 equations for 0 unknowns...
Initialization status: OVERDETERMINED
Non-trivial mass matrix: true
retcode: InitialFailure

Inspecting the initialization equations shows:

equations(prob.f.initializeprob.f.sys)
# → returns 0 = -R1.i(t) + R2.i(t) + D1.i(t)

So from what I understand, the initialization system is overspecified, but I don’t get why — especially since this same structure worked in MTK v8.

:red_question_mark: My question:

Any tips on how to make this system “fully determined” in MTK v9 to avoid this initialization error?

I’d really appreciate any guidance or pointer :folded_hands:

Thanks a lot!

2

Since there is only one differential equation, you only need one initial condition.

That second condition is extra information. There’s only one condition to solve for but two conditions.

u0 = [
    C2.Δv => 7209,
]

That would fix it. You might want to do:

prob = ODEProblem(sys, u0, (0, 10), guesses = [R2.Δv => 1671])

basically, help the system by giving a guess for the algebraic value, but it’s not necessary and it may not be used.

MTK v8 used to not be smart enough to validate this. It would automatically treat the non-differential value R2.Δv => 1671 as a guess, but this confused users because sometimes they would see it change and say “but I told it I wanted this!”

So in MTK v9, we guarantee that if you say “this value is 1.2354”, then it will be that value. If it can’t be that value, we will error. We guarantee that all statements and equations will be respected or error if they cannot simultaneously. If you’re trying to pass information to help the Newton method guess a good starting point, that’s separted as guess so that it’s explicitly not an equation that is guaranteed. That’s the major breaking change of v9, force that all equations are respected. In the end, it gives much cleaner explanations of what each term means, but it does mean that some things that seem to have worked in v8 will have broke in v9 (though I stress “seemed” because if you check sol[1] you’ll see R2.Δv => 1671 didn’t hold, it just wasn’t clear that initial values on algebraic variables were interpreted as guesses).

3

Thanks a lot, Chris! That cleared things up perfectly — especially the distinction between initial conditions and guesses in MTK v9. Really appreciate your help!