Fix up initializesystem for hierarchical models

[ChrisRackauckas]

The first version only worked on non-hierarchical models. Now it handles all models. Here's some example usage. Here's an example model:

using ModelingToolkit, DifferentialEquations

@parameters t
D = Differential(t)

@connector Port begin
    p(t)
    dm(t)=0, [connect = Flow]
end

@connector Flange begin
    dx(t)=0
    f(t), [connect = Flow]
end

# Components ----
@mtkmodel Orifice begin
    @parameters begin
        Cₒ=2.7
        Aₒ=0.00094
        ρ₀=1000
        p′=0
    end
    @variables begin
        dm(t)=0
        p₁(t)=p′
        p₂(t)=p′
    end
    @components begin
        port₁ = Port(p=p′)
        port₂ = Port(p=p′)
    end
    begin
        u = dm/(ρ₀*Aₒ)
    end
    @equations begin
        dm ~ +port₁.dm
        dm ~ -port₂.dm
        p₁ ~ port₁.p
        p₂ ~ port₂.p

        p₁ - p₂ ~ (1/2)*ρ₀*u^2*Cₒ
    end
end

@mtkmodel Volume begin
    @parameters begin
        A=0.1
        ρ₀=1000
        β=2e9
        direction=+1
        p′
        x′
    end
    @variables begin
        p(t)=p′
        x(t)=x′
        dm(t)=0
        f(t)=p′ * A
        dx(t)=0
        r(t), [guess = 1000]
        dr(t), [guess = 1000]
    end
    @components begin
        port = Port(p=p′)
        flange = Flange(f=-p′ * A * direction)
    end
    @equations begin
        D(x) ~ dx
        D(r) ~ dr

        p ~ +port.p
        dm ~ +port.dm # mass is entering
        f ~ -flange.f * direction # force is leaving
        dx ~ flange.dx * direction

        r ~ ρ₀*(1 + p/β)
        dm ~ (r*dx*A) + (dr*x*A)
        f ~ p * A
    end
end

@mtkmodel Mass begin
    @parameters begin
        m = 100
        f′
    end
    @variables begin
        f(t)=f′
        x(t)=0
        dx(t)=0
        ẍ(t)=f′/m
    end
    @components begin
        flange = Flange(f=f′)
    end
    @equations begin
        D(x) ~ dx
        D(dx) ~ ẍ

        f ~ flange.f
        dx ~ flange.dx

        m*ẍ ~ f
    end
end

@mtkmodel Actuator begin
    @parameters begin
        p₁′
        p₂′
    end
    begin #constants
        x′=0.5
        A=0.1
    end
    @components begin
        port₁ = Port(p=p₁′)
        port₂ = Port(p=p₂′)
        vol₁ = Volume(p′=p₁′, x′=x′,  direction=-1)
        vol₂ = Volume(p′=p₂′, x′=x′,  direction=+1)
        mass = Mass(f′=(p₂′ - p₁′)*A)
        flange = Flange(f=0)
    end
    @equations begin
        connect(port₁, vol₁.port)
        connect(port₂, vol₂.port)
        connect(vol₁.flange, vol₂.flange, mass.flange, flange)
    end
end

@mtkmodel Source begin
    @parameters begin
        p′
    end
    @components begin
        port = Port(p=p′)
    end
    @equations begin
        port.p ~ p′
    end
end

@mtkmodel Damper begin
    @parameters begin
        c = 1000
    end
    @components begin
        flange = Flange(f=0)
    end
    @equations begin
        flange.f ~ c*flange.dx
    end
end

@mtkmodel System begin
    @components begin
        res₁ = Orifice(p′=300e5)
        res₂ = Orifice(p′=0)
        act = Actuator(p₁′=300e5, p₂′=0)
        src = Source(p′=300e5)
        snk = Source(p′=0)
        dmp = Damper()
    end
    @equations begin
        connect(src.port, res₁.port₁)
        connect(res₁.port₂, act.port₁)
        connect(act.port₂, res₂.port₁)
        connect(res₂.port₂, snk.port)
        connect(dmp.flange, act.flange)
    end
end

@mtkbuild sys = System()

It's having troubles initializing, so what is the system it's trying to initialize?

initsys = initializesystem(sys)

That gives it symbolically. We have 98 equations for 54 variables, shesh that's overloaded.

julia> initprob = NonlinearProblem(initsys,[])
ERROR: ArgumentError: Equations (98), states (54), and initial conditions (54) are of different lengths. To allow a different number of equations than states use kwarg check_length=false.
Stacktrace:
 [1] check_eqs_u0(eqs::Vector{…}, dvs::Vector{…}, u0::Vector{…}; check_length::Bool, kwargs::@Kwargs{})
   @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\abstractsystem.jl:1871
 [2] process_NonlinearProblem(constructor::Type, sys::NonlinearSystem, u0map::Vector{…}, parammap::SciMLBase.NullParameters; version::Nothing, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, kwargs::@Kwargs{…})
   @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:339
 [3] process_NonlinearProblem
   @ C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:323 [inlined]
 [4] (NonlinearProblem{…})(sys::NonlinearSystem, u0map::Vector{…}, parammap::SciMLBase.NullParameters; check_length::Bool, kwargs::@Kwargs{})
   @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:368
 [5] NonlinearProblem (repeats 2 times)
   @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:365 [inlined]
 [6] #NonlinearProblem#698
   @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:362 [inlined]
 [7] NonlinearProblem(sys::NonlinearSystem, args::Vector{Any})
   @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:361
 [8] top-level scope
   @ REPL[1]:1
Some type information was truncated. Use `show(err)` to see complete types.

So we can't build a NonlinearProblem. But we can do this:

initprob = NonlinearProblem(initsys,[], check_length=false, checkbounds = true)
lsqprob = NonlinearLeastSquaresProblem(NonlinearFunction(initprob.f.f, resid_prototype = zeros(98)), initprob.u0, initprob.p)
initsol = solve(lsqprob, GaussNewton(), reltol = 1e-12, abstol = 1e-12)
retcode: Success
u: 54-element Vector{Float64}:
      0.5
   1015.0
      0.5
   1000.0
      0.0
     -1.208682818218319e-28
      ⋮
      5.048709793414476e-28
     -3.52499105861307e-29
     -1.5146129380243427e-28
     -3.0e6
 -30000.0
     -3.0e6

And now we can make our initial conditions match that:

allinit = states(initsys) .=> initsol
prob = ODEProblem(sys, allinit, (0,0.1))

and solve:

sol = solve(prob)
tmp = [0.0, 0.0, -0.0, -0.0, -0.0, -0.0, 8.123585990009498e-54, 1.1599794698483246e-51, 1.5375323209568473e-26, -2.914685306392616e-26]
┌ Warning: Instability detected. Aborting
└ @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\3Mujd\src\integrator_interface.jl:602

Here is modified the solver so tmp is the resid vector inside of the initializer. You can see it worked because the initialization errors are all zero. This is as opposed to with the default initial conditions:

sysguesses = [ModelingToolkit.getguess(st) for st in states(sys)]
hasaguess = findall(!isnothing, sysguesses)
sysguesses = Dict(states(sys)[hasaguess] .=> sysguesses[hasaguess])
prob = ODEProblem(sys, sysguesses, (0,0.1))
sol = solve(prob)

julia> sol = solve(prob)
tmp = [-0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.0e7, 0.0, 50.0, 50.0]
nlsol.resid = [-3.0e7, 0.0, 50.0, 50.0]
nlsol.retcode = SciMLBase.ReturnCode.MaxIters
retcode: InitialFailure
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.5, 1000.0, 0.5, 1000.0, 0.0, 0.0, 0.0, 0.0, 1000.0, 1000.0]

The model still fails since the initial mass is zero but the derivative of the mass is negative, and thus we have negative mass right. Unless this is for an interdimensional space ship that means this model is incorrect, so that's why it has this issue, but it serves as a fun case to play with.

Todo:

• Make the guess part of the system and use ModelingToolkit.guesses(sys) in the initailization piece
• Allow for passing geuss in ODESystem
• Put a nicer syntax on building a NonlinearLeastSquaresProblem, since that is a nice utility to have around
• Integrate this into ODEProblem so that the initialization automatically applies this nonlinear system

Fixes #2151 fixes #2103 fixes #2000 fixes #1836 fixes #1565

[YingboMa]

@ChrisRackauckas

initsys = initializesystem(sys)
ss = structural_simplify(initsys, fully_determined = false)
initprob = NonlinearProblem(ss, [], check_length=false, checkbounds = true)
lsqprob = NonlinearLeastSquaresProblem(NonlinearFunction(initprob.f.f, resid_prototype = zeros(length(equations(ss)))), initprob.u0, initprob.p)
initsol = solve(lsqprob, GaussNewton(), reltol = 1e-12, abstol = 1e-12)

now works. But

allinit = states(initsys) .=> initsol[states(initsys)]

doesn't because the observed function isn't hooked in

julia> initsol.prob.f.observed
DEFAULT_OBSERVED_NO_TIME (generic function with 1 method)

[bradcarman]

I'm putting together a dead simple case to understand how to use and apply the initializesystem. This is a simple mass-damper model. I want to specify the initial velocity of the mass and have the remaining unset states determined, for example the initial velocity of the damper.dx should match the mass.dx. Here is the model

using ModelingToolkit

@parameters t
D = Differential(t)

@connector Flange begin
    dx(t), [guess = 0]
    f(t), [guess = 0, connect=Flow]
end

@mtkmodel Mass begin
    @parameters begin
        m = 100
    end
    @variables begin
        dx(t)
        f(t)=0
    end
    @components begin
        flange = Flange()
    end
    @equations begin
        # connectors
        flange.dx ~ dx
        flange.f ~ -f
        
        # physics
        f ~ m*D(dx)
    end
end

@mtkmodel Damper begin
    @parameters begin
        d = 1
    end
    @variables begin
        dx(t), [guess = 0]
        f(t), [guess = 0]
    end
    @components begin
        flange = Flange()
    end
    @equations begin
        # connectors
        flange.dx ~ dx
        flange.f ~ -f
        
        # physics
        f ~ d*dx
    end
end

@mtkmodel System begin
    @components begin
        mass = Mass(;dx=100,m=10)
        damper = Damper(;d=1)
    end
    @equations begin
        connect(mass.flange, damper.flange)
    end
end

@named odesys = System()
equations(odesys)

@mtkbuild sys = System()
isys = initializesystem(sys)

This gives the overdetermined system with 9 equations and 8 states. If I look at the equations I see

julia> eqs = equations(isys)
9-element Vector{Equation}:
 0 ~ 100 - mass₊dx(t)
 0 ~ -mass₊f(t)
 0 ~ mass₊dx(t) - mass₊flange₊dx(t)
 0 ~ mass₊dx(t) - damper₊dx(t)
 0 ~ -damper₊f(t) + damper₊d*damper₊dx(t)
 0 ~ -damper₊flange₊dx(t) + damper₊dx(t)
 0 ~ -damper₊f(t) - mass₊f(t)
 0 ~ -damper₊f(t) - damper₊flange₊f(t)
 0 ~ -mass₊f(t) - mass₊flange₊f(t)

It seems the extra equation is: 0 ~ mass₊dx(t) - damper₊dx(t). I'm not sure where this equation is coming from, it is not directly defined in the system (it can only be inferred).

[MarcBerliner]

allinit = states(initsys) .=> initsol[states(initsys)]

doesn't [work] because the observed function isn't hooked in

julia> initsol.prob.f.observed
DEFAULT_OBSERVED_NO_TIME (generic function with 1 method)

@YingboMa this works if we replace resid_prototype in the NonlinearFunction and keep all the other fields

using Setfield: @set!
initsys = initializesystem(sys)

ss = structural_simplify(initsys, fully_determined = false)
initprob = NonlinearProblem(ss, [], check_length=false, checkbounds = true)

nl_f = initprob.f
@set! nl_f.resid_prototype = zeros(length(equations(ss)))
lsqprob = NonlinearLeastSquaresProblem(nl_f, initprob.u0, initprob.p)

initsol = solve(lsqprob, GaussNewton(), reltol = 1e-12, abstol = 1e-12)

allinit = states(initsys) .=> initsol[states(initsys)] # 54-element Vector{Pair{...}}

[YingboMa]

Ah, of course.

[ChrisRackauckas]

It seems the extra equation is: 0 ~ mass₊dx(t) - damper₊dx(t). I'm not sure where this equation is coming from, it is not directly defined in the system (it can only be inferred).

That's just one of the connection equations. It needs to be there.

[ChrisRackauckas]

Okay, this got derailed because I was being fed incorrect information in the examples 😅 . I was waiting for a bit of time to really dive in and find out where extra equations were coming from, turns out there are none and @bradcarman 's example is just overdetermined. It properly found all of the stated equations, and if you want a not over-determined system you could do something like:

using ModelingToolkit

@parameters t
D = Differential(t)

@connector Flange begin
    dx(t), [guess = 0]
    f(t), [guess = 0, connect=Flow]
end

@mtkmodel Mass begin
    @parameters begin
        m = 100
    end
    @variables begin
        dx(t)
        f(t), [guess = 0]
    end
    @components begin
        flange = Flange()
    end
    @equations begin
        # connectors
        flange.dx ~ dx
        flange.f ~ -f
        
        # physics
        f ~ m*D(dx)
    end
end

@mtkmodel Damper begin
    @parameters begin
        d = 1
    end
    @variables begin
        dx(t), [guess = 0]
        f(t), [guess = 0]
    end
    @components begin
        flange = Flange()
    end
    @equations begin
        # connectors
        flange.dx ~ dx
        flange.f ~ -f
        
        # physics
        f ~ d*dx
    end
end

@mtkmodel System begin
    @components begin
        mass = Mass(;dx=100,m=10)
        damper = Damper(;d=1)
    end
    @equations begin
        connect(mass.flange, damper.flange)
    end
end

@named odesys = System()
equations(odesys)

@mtkbuild sys = System()
isys = initializesystem(sys) |> structural_simplify

initprob = NonlinearProblem(isys,[], check_length=false, checkbounds = true)
initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12)

and this all works. So 🎉 we're good.

[ChrisRackauckas]

I think the only issue left now is that NonlinearLeastSquaresProblem is having issue with observed since it thinks it's time dependent:

using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

@connector Port begin
    p(t)
    dm(t) = 0, [connect = Flow]
end

@connector Flange begin
    dx(t) = 0
    f(t), [connect = Flow]
end

# Components ----
@mtkmodel Orifice begin
    @parameters begin
        Cₒ = 2.7
        Aₒ = 0.00094
        ρ₀ = 1000
        p′ = 0
    end
    @variables begin
        dm(t) = 0
        p₁(t) = p′
        p₂(t) = p′
    end
    @components begin
        port₁ = Port(p = p′)
        port₂ = Port(p = p′)
    end
    begin
        u = dm / (ρ₀ * Aₒ)
    end
    @equations begin
        dm ~ +port₁.dm
        dm ~ -port₂.dm
        p₁ ~ port₁.p
        p₂ ~ port₂.p

        p₁ - p₂ ~ (1 / 2) * ρ₀ * u^2 * Cₒ
    end
end

@mtkmodel Volume begin
    @parameters begin
        A = 0.1
        ρ₀ = 1000
        β = 2e9
        direction = +1
        p′
        x′
    end
    @variables begin
        p(t) = p′
        x(t) = x′
        dm(t) = 0
        f(t) = p′ * A
        dx(t) = 0
        r(t), [guess = 1000]
        dr(t), [guess = 1000]
    end
    @components begin
        port = Port(p = p′)
        flange = Flange(f = -p′ * A * direction)
    end
    @equations begin
        D(x) ~ dx
        D(r) ~ dr

        p ~ +port.p
        dm ~ +port.dm # mass is entering
        f ~ -flange.f * direction # force is leaving
        dx ~ flange.dx * direction

        r ~ ρ₀ * (1 + p / β)
        dm ~ (r * dx * A) + (dr * x * A)
        f ~ p * A
    end
end

@mtkmodel Mass begin
    @parameters begin
        m = 100
        f′
    end
    @variables begin
        f(t) = f′
        x(t) = 0
        dx(t) = 0
        ẍ(t) = f′ / m
    end
    @components begin
        flange = Flange(f = f′)
    end
    @equations begin
        D(x) ~ dx
        D(dx) ~ ẍ

        f ~ flange.f
        dx ~ flange.dx

        m * ẍ ~ f
    end
end

@mtkmodel Actuator begin
    @parameters begin
        p₁′
        p₂′
    end
    begin #constants
        x′ = 0.5
        A = 0.1
    end
    @components begin
        port₁ = Port(p = p₁′)
        port₂ = Port(p = p₂′)
        vol₁ = Volume(p′ = p₁′, x′ = x′, direction = -1)
        vol₂ = Volume(p′ = p₂′, x′ = x′, direction = +1)
        mass = Mass(f′ = (p₂′ - p₁′) * A)
        flange = Flange(f = 0)
    end
    @equations begin
        connect(port₁, vol₁.port)
        connect(port₂, vol₂.port)
        connect(vol₁.flange, vol₂.flange, mass.flange, flange)
    end
end

@mtkmodel Source begin
    @parameters begin
        p′
    end
    @components begin
        port = Port(p = p′)
    end
    @equations begin
        port.p ~ p′
    end
end

@mtkmodel Damper begin
    @parameters begin
        c = 1000
    end
    @components begin
        flange = Flange(f = 0)
    end
    @equations begin
        flange.f ~ c * flange.dx
    end
end

@mtkmodel System begin
    @components begin
        res₁ = Orifice(p′ = 300e5)
        res₂ = Orifice(p′ = 0)
        act = Actuator(p₁′ = 300e5, p₂′ = 0)
        src = Source(p′ = 300e5)
        snk = Source(p′ = 0)
        dmp = Damper()
    end
    @equations begin
        connect(src.port, res₁.port₁)
        connect(res₁.port₂, act.port₁)
        connect(act.port₂, res₂.port₁)
        connect(res₂.port₂, snk.port)
        connect(dmp.flange, act.flange)
    end
end

@mtkbuild sys = System()

initprob = ModelingToolkit.InitializationProblem(sys)
sol = solve(initprob)
getter = ModelingToolkit.getu(initprob, unknowns(sys))
getter(sol)

[ChrisRackauckas]

Okay, initialization stuff is now complete but requires SciML/OrdinaryDiffEq.jl#2151

[bradcarman]

As you point out, there is error in the algebraic constraints. We can see this by running

julia> eq = full_equations(sys)[7]
       ModelingToolkit.fixpoint_sub(eq.rhs, defs)
-2.2724270820617676e-7

However, the model is not "wrong", this is a floating point error issue, which we can see this if we do manual substitutions rather than using fixpoint_sub:

julia> eq = full_equations(sys)[7]
       defs = ModelingToolkit.defaults(sys)
       defs[sys.res₁.ṁ]
       eq = substitute(eq, sys.res₁.ṁ => defs[sys.res₁.ṁ])
       defs[sys.act.vol₁.r]
       eq = simplify(substitute(eq, sys.act.vol₁.r => defs[sys.act.vol₁.r]))
       defs[sys.act.vol₁.p′]
       eq = substitute(eq, sys.act.vol₁.p′ => defs[sys.act.vol₁.p′])
       defs[sys.act.p₁′]
       eq = substitute(eq, sys.act.p₁′ => defs[sys.act.p₁′])
       defs[sys.src.p′]
       eq = substitute(eq, sys.src.p′ => defs[sys.src.p′])
0 ~ 0.0

The floating point error is originating from the density variable sys.act.vol₁.r, which is proven by the fact that fixedpoint_sub gives:

julia> ModelingToolkit.fixpoint_sub(defs[sys.act.vol₁.r], defs)
1014.9999999999999

But the correct value is

julia> @parameters r_new
       eq = full_equations(sys)[7]
       eq_new = Symbolics.substitute(eq.rhs, sys.act.vol₁.r => r_new)
       eq_new = ModelingToolkit.fixpoint_sub(eq_new, defs)
       Symbolics.solve_for(eq_new, r_new)
1015.0

Am I able to use the initializationsystem to resolve these floating point error issues?

[ChrisRackauckas]

It cannot resolve that because it cannot assume you are lying to it. If you state all of the equalities that you have, there is no solution to more than 2e-7. That is absolutely clear. Yes, it is floating point error. You can remove it by changing some parameters to higher precision, etc. there are things you can do. Or instead of holding some of these extra equalities, you can relax some to be guesses.

The point is that as a system, ModelingToolkit cannot tell you which equality to remove. That's fundamentally a modeler's choice. If ModelingToolkit ever decided to allow a != b when you said a ~ b, that would be a bug. If the system that you've written down is not satisfiable, it should tell you that. That's exactly what it's doing.

[ChrisRackauckas]

Okay I think I fixed all of the test cases. Some particular things to note:

1. I found an interface edge case that can show up for these simple models Trivial (length 0 state) NonlinearLeastSquaresProblem always returns success NonlinearSolve.jl#387
2. In some of the models it was clear that it was mixing guesses and u0, given the previous semantics. Now that defaults are held truly, these would error as unsatisfiable. This for example was seen in the index reduction case with the pendulum, where T => 0 was previously a guess for the algebraic variable but now that's created as T ~ 0 in the initial system, as it will satisfy what the user wrote, and thus the initialization fails. The fix of course is to properly mark that as a guess.
3. One of the test cases in the Components.jl one had the wrong sign. I had to analytically dig through to confirm, and indeed it was wrong. It only passed before because the solution ended up being a trivial zero because the initialization mixed guess semantics, and the initial 0 seems unintended. So there was some bug fixing there 😅
4. The state selection catapult test was wildly overconstrained, moving things to guesses found a reasonable solution.

All-in-all, I think that some folks may need to fix up their component-based models because now that u0 is truly "I want this initial condition", some people may get failures because what they say is just not possible. But of course, this is what we want, we want to make sure anything the user asks for is either satisfied or fails as unsatisfiable.

Though if your initialization has 0 states SciML/NonlinearSolve.jl#387 it doesn't fail 😅 but I'll get that fixed separately.

[ChrisRackauckas]

Tests all passed except master, going for it!