Initialization question for the RC system with HeatingResistor

[wang890]

Question❓

The MTK code for initialization may change starting from September 2024. Before that, for example, in version 9.19.0, the following MRE code of RC system with HeatingResistor worked.

Although I went through the tutorial Parameter initialization by example and have done some tests (github page of the Jupyter file), I still don't know how to debug it according to the physical laws of RC. Could you help me, please? Thanks a lot!

Error: Initialization system is overdetermined. 1 equations for 0 unknowns.
isys3_prob.f.initializeprob.f.sys out:
0 =−resistor.v(t)+resistor.i(t)resistor.R(t), if deleted this equation in the model of HeatingResistor, Error:
ExtraVariablesSystemException: The system is unbalanced.
There are 25 highest order derivative variables and 24 equations.

If using the common Resistor model without heat_port, there is no error.

ENV

  julia 1.11.0
  [961ee093] ModelingToolkit v9.72.0
  [1dea7af3] OrdinaryDiffEq v6.93.0
  [91a5bcdd] Plots v1.40.11
  [0c5d862f] Symbolics v6.37.1

Jupyter test file of the MRE is on the github page
Here is the copy-paste MRE:

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkit: t_nounits as t, D_nounits as d

@connector Pin begin
    v(t) # Potential at the pin
    i(t), [connect = Flow] # Current flowing into the pin
end

@mtkmodel OnePort begin
    @components begin
        p = Pin()
        n = Pin()
    end
    @variables begin
        v(t) # Voltage drop of the two pins (= p.v - n.v)
        i(t) # Current flowing into the pin
    end
    @equations begin
        v ~ p.v - n.v
        0 ~ p.i + n.i
        i ~ p.i
    end
end

@connector HeatPort begin
    T(t)
    Q_flow(t), [connect = Flow]
end

@connector HeatPort_a begin
    @extend T, Q_flow = heatPort_a = HeatPort()
end

@connector HeatPort_b begin
    @extend T, Q_flow = heatPort_b = HeatPort()
end

@mtkmodel FixedTemperature begin
    @components begin
        port = HeatPort_b()
    end
    @parameters begin
        T, [description = "Fixed temperature boundary condition"]
    end
    @equations begin
        port.T ~ T
    end
end

# copy from ModelingToolkitStandardLibrary
@mtkmodel HeatingResistor begin
    @extend v, i = oneport = OnePort()
    @components begin
        heat_port = HeatPort()
    end
    @parameters begin
        R_ref = 1.0, [description = "Reference resistance"]
        T_ref = 300.15, [description = "Reference temperature"]
        alpha = 0, [description = "Temperature coefficient of resistance"]
    end
    @variables begin
        R(t)#, [guess = R_ref]
    end
    @equations begin
        R ~ R_ref * (1 + alpha * (heat_port.T - T_ref))
        heat_port.Q_flow ~ -v * i # -LossPower
        v ~ i * R # if delete this equation, ExtraVariablesSystemException: The system is unbalanced. 
        # There are 25 highest order derivative variables and 24 equations.
    end
end

@mtkmodel Ground begin
    @components begin
        p = Pin()
    end
    @equations begin
        p.v ~ 0
    end
end

@mtkmodel Capacitor begin
    @extend i,v = oneport = OnePort()
    @parameters begin
        C # Capacitance
    end
    @equations begin
        i ~ C * d(v)
    end
end

@mtkmodel ConstantVoltage begin
    @extend v, = oneport = OnePort()    
    @parameters begin
        V # Value of constant voltage
    end
    @equations begin
        v ~ V
    end
end

@mtkmodel RcHeating begin    
    @components begin
        resistor = HeatingResistor()           
        capacitor = Capacitor()
        source = ConstantVoltage()
        ground = Ground()
        fixedTemperature = FixedTemperature()
        
    end      
    @equations begin        
        connect(source.p, resistor.p)
        connect(resistor.n, capacitor.p)
        connect(capacitor.n, source.n, ground.p) 
        connect(resistor.heat_port, fixedTemperature.port)       
    end
end

@mtkbuild sys = RcHeating()
"""
Unknowns (2)
    capacitor₊v(t) ▪ Unassigned ▪ Voltage drop of the two pins (= p.v - n.v)
    resistor₊i(t) ▪ Unassigned ▪ Current flowing from pin p to pin n
Parameters (6)
    resistor₊alpha ▪ 0 ▪ Temperature coefficient of resistance
    fixedTemperature₊T ▪ Unassigned ▪ Fixed temperature at port
    source₊V ▪ Unassigned ▪ Value of constant voltage
    resistor₊R_ref ▪ 1.0 ▪ Reference resistance
    resistor₊T_ref ▪ 300.15 ▪ Reference temperature
    capacitor₊C ▪ Unassigned ▪ Capacitance
Equations (2)
    Differential(t)(capacitor₊v(t)) ~ capacitor₊i(t) / capacitor₊C
    0 ~ -resistor₊v(t) + resistor₊i(t)*resistor₊R(t)
"""

isys3 = generate_initializesystem(sys; 
            u0map = [sys.capacitor.v => 0.0, sys.resistor.i => 0.0],
            pmap = [sys.capacitor.C => 1.0, sys.source.V=>1.0,sys.fixedTemperature.T=>400.15],
            guesses = [sys.resistor.R=>sys.resistor.R_ref])

isys3_prob = ODEProblem(sys, [sys.capacitor.v => 0.0, sys.resistor.i => 0.0], (0.0, 10.0),
                [sys.capacitor.C => 1.0, sys.source.V=>1.0,sys.fixedTemperature.T=>400.15],
                guesses = [sys.resistor.R=>sys.resistor.R_ref])
                                
isys3_prob_initsys = isys3_prob.f.initializeprob.f.sys
# out: 0 =−resistor.v(t)+resistor.i(t)resistor.R(t)

if without guesses, the error is the same.

[ChrisRackauckas]

I feel like I answered this one somewhere else? It is overdetermined because u0map = [sys.capacitor.v => 0.0, sys.resistor.i => 0.0] that's too many initial conditions. Did you mean u0map = [sys.capacitor.v => 0.0]? Before we would drop the algebraic condition, see https://discourse.julialang.org/t/overdetermined-odeproblem-after-upgrading-from-mtk-v8-75-0-to-mtk-v9-72-0/128163/2 . We got better at validating it, and so we throw errors now when it (or any equation) is violated since sometime early last year. The answer is to make sure all set equations are correct, and I highly recommend fixing systems that are overdetermined since they are a major source of building unsatisfiable systems. If the system is unsatisfiable mathematically, there's simply nothing we can do.

[wang890]

I feel like I answered this one somewhere else? It is overdetermined because u0map = [sys.capacitor.v => 0.0, sys.resistor.i => 0.0] that's too many initial conditions. ... ...

Thanks a lot. I got it.