solving problem with discrete events

[lixiao00093]

I think this relates to my previous issues: #3043, and has some relationship to SciML/OrdinaryDiffEq.jl#2441.

• If I solve it without discrete events, it can be solved with the default solver, but if I use Rodas4() it shows "Unstable" and the solution is unsuccessful.
• If I solve it with discrete events, when I use the default solver, I get the error 'ERROR: UndefRefError: access to undefined reference'. Even if you then switch to another solver and run it again, the initialisation error remains. But if you use Rodas4() to solve the problem from the beginning, it will show "Unstable", the solver is not successful either.

I am stuck with this problem for a long time, no matter I use the latest version 9.41.0 or the older version of the ModelingToolkit.

What should I do? Should I modify my codes? Or this is a bug?

• Minimal Reproducible Example 👇

using ModelingToolkit, DifferentialEquations, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D_t

# --------system models --------
@connector  AcBus begin
    @variables begin
        vm(t),  [guess = 1.0]
        va(t),  [guess = 0.0]
        pin(t), [connect = Flow, guess = 0.0]
        qin(t), [connect = Flow, guess = 0.0]
    end
end

@mtkmodel Zipload_component begin
    @parameters begin 
        u = 1,[description = "Enable Signal"]
        g
        b
    end
    @components begin
        bus = AcBus()
    end
    @equations begin
        bus.pin ~ u * (  g * bus.vm ^2)
        bus.qin ~ u * (- b * bus.vm ^2)
    end
end

@mtkmodel Piline_component begin
    @parameters begin
        u = 1 ,[description = "Enable Signal"]
        rft
        xft
        g0 = 0.0
        b0
        tap_ratio = 1.0
        gft =  rft/(rft^2 + xft^2)
        bft = -xft/(rft^2 + xft^2)
        m_real = tap_ratio 
        m2 = tap_ratio^2
    end

    @components begin
        bus_fr = AcBus()
        bus_to = AcBus()
    end

    @equations begin
        bus_fr.pin*m2 ~ u*( g0/2*bus_fr.vm^2 + gft*bus_fr.vm^2
                            - gft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
                            - bft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
        -bus_fr.qin*m2 ~ u*(b0/2*bus_fr.vm^2 + bft*bus_fr.vm^2
                            - bft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
                            + gft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
        
        bus_to.pin*m2 ~ u*( g0/2*m2*bus_to.vm^2 + gft*m2*bus_to.vm^2
                            - gft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
                            + bft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
        -bus_to.qin*m2 ~ u*(b0/2*m2*bus_to.vm^2 + bft*m2*bus_to.vm^2 
                            - bft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va) 
                            - gft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
    end
end

@mtkmodel Syn2_component begin
    @components begin
        bus = AcBus()
    end

    @variables begin
        delta(t), [description = "rotor angle"]
        omega(t), [description = "rotor speed", guess = 1.0]

        vf(t)
        tm(t)
        p(t)
        q(t)

        _Id(t), [guess = 0.0]
        _Iq(t), [guess = 0.0]
    end

    @parameters begin 
        u = 1, [description = "Enable Signal"]
        fn
        H
        Damp
        ra
        xd1
        Kw
        Kp
        iM = u / (2 * H)
        c1 = ra / (ra^2 + xd1^2)
        c2 = xd1 / (ra^2 + xd1^2)
        c3 = c2
        vm0
        va0
    end

    @equations begin
        _Id ~ u*(-(bus.vm * sin(delta - bus.va) * c1) - (bus.vm * cos(delta - bus.va) * c3) + c3 * vf)
        _Iq ~ u*( (bus.vm * sin(delta - bus.va) * c2) - (bus.vm * cos(delta - bus.va) * c1) + c1 * vf)

        p ~ bus.vm * sin(delta - bus.va) * _Id + bus.vm * cos(delta - bus.va) * _Iq
        q ~ bus.vm * cos(delta - bus.va) * _Id - bus.vm * sin(delta - bus.va) * _Iq
        bus.pin ~ -p                                # for AcBus
        bus.qin ~ -q                                # for AcBus

        D_t(delta) ~ u * 2 * π * fn * (omega - 1)
        D_t(omega) ~ iM * (tm - p - ra * (_Id^2 + _Iq^2) - Damp * (omega - 1))
        D_t(tm) ~ 0
        D_t(vf) ~ 0
    end
end

@mtkmodel Fault_component begin
    @parameters begin
        u = 0,[description="fault status"]
        rf
        xf
        yg = rf/(rf^2 + xf^2)
        yb = -xf/(rf^2 + xf^2)
    end

    @components begin
        bus = AcBus()
    end

    @equations begin
        bus.pin ~  u * yg * bus.vm^2
        bus.qin ~ -u * yb * bus.vm^2
    end
end

# -------colection of instances of models--------
Bus_Dict = Dict{Symbol,ODESystem}(Symbol(1)=>AcBus(;name = Symbol("Bus 1")),Symbol(2)=>AcBus(;name = Symbol("Bus 2")))
device_list::Dict{Symbol,ODESystem} = Dict(
    :line => Piline_component(;name = :line, rft = 0.01,xft = 0.01,g0 = 0,b0 = 0.01,tap_ratio = 1.0),
    :syn => Syn2_component(;name = :syn1, fn = 50.0, H = 1.0, Damp = 0.8, ra = 0.01, xd1 = 0.1, Kw=1.0, Kp=1.0, va0 = 0.0 , vm0=1.0),
    :zipload => Zipload_component(;name = :zipload1, g = 1.0, b = 0.0),
    :fault => Fault_component(;name = :fault, rf = 0.0, xf = 0.1)
)

# -------connection of models--------
eqs = [ connect(Bus_Dict[Symbol(1)], device_list[:syn].bus),
        connect(Bus_Dict[Symbol(1)], device_list[:line].bus_fr),
        connect(Bus_Dict[Symbol(1)], device_list[:fault].bus),
        
        connect(Bus_Dict[Symbol(2)], device_list[:zipload].bus),
        connect(Bus_Dict[Symbol(2)], device_list[:line].bus_to)]

# -------initialization of models--------
ini_eqs = [device_list[:syn].omega ~ 1,
    device_list[:syn].tm - device_list[:syn].p - device_list[:syn].ra * (device_list[:syn]._Id^2 + device_list[:syn]._Iq^2) ~ 0,
    device_list[:syn].bus.va ~ device_list[:syn].va0, 
    device_list[:syn].bus.vm ~ device_list[:syn].vm0]

# -------change the parameters of models--------
# -------Apply and clear fault in the models--------
u_temp_record = Dict{Int,Vector{Float64}}()
function affect1!(integ, u, p, ctx)
    println("apply fault at t = ", integ.t)
    if integ.ps[p.u] != ctx
        integ.ps[p.u] = ctx
        u_temp_record[1] = copy(integ.u[integ.differential_vars.==false])
    end
end
function affect2!(integ, u, p, ctx)
    println("clear fault at t = ", integ.t)
    if integ.ps[p.u] != ctx
        integ.ps[p.u] = ctx
        #offer the initial value of algebraic variables for the reinitialization
        integ.u[integ.differential_vars.==false] = u_temp_record[1]
    end
end

fault_function = [(t==1.0) => (affect1!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 1.0),
                 (t==1.05) => (affect2!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 0.0)]

# -------set the initial guessvalues of models--------
function set_ini_syn2(; ra, xd1)
    pg0 = 1.0
    qg0 = 0.0
    vm0 = 1.0
    va0 = 0.0

    V = vm0 .* exp.(va0 .* 1im)
    S = pg0 - qg0 .* 1im
    I = S ./ conj.(V)
    E = V + (ra + xd1 .* 1im) .* I
    delta = imag(log.(E ./ abs.(E) + 0im))
    delta0 = delta
    omega0 = 1.0

    # d- and q-axis voltages and currents
    vdq = V .* 1im .* exp.(-delta .* 1im)
    idq = I .* 1im .* exp.(-delta .* 1im)
    vd = real(vdq)
    vq = imag(vdq)
    Id = real(idq)
    Iq = imag(idq)
    # mechanical torques/powers
    tm0 = (vq + ra .* Iq) .* Iq + (vd + ra .* Id) .* Id
    vf0 =  vq + ra .* Iq + xd1 .* Id
    return delta0, omega0, tm0, vf0, pg0, qg0
end
u0_guess = Dict(zip([device_list[:syn].delta, device_list[:syn].omega, device_list[:syn].tm, device_list[:syn].vf, device_list[:syn].p, device_list[:syn].q], set_ini_syn2(ra = 0.01, xd1 = 0.1)))

# -------simplify the system--------
test_system_init = ODESystem(eqs, t,[],[]; systems = collect(values(merge(device_list, Bus_Dict))), name = :test_system_init, discrete_events = fault_function)
test_system_init_simplify = structural_simplify(test_system_init, fully_determined = true)

# -------format the system to solve--------
prob = ODEProblem(test_system_init_simplify,[], (0,20),[]; guesses = u0_guess, initialization_eqs = ini_eqs, fully_determined = true, jac=true,sparse=true)
# sol = solve(prob, Rodas4(), tstops = [1.0,1.05])
sol = solve(prob, tstops = [1.0,1.05])
plot(sol)

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [6e4b80f9] BenchmarkTools v1.5.0
  [336ed68f] CSV v0.10.14
  [a93c6f00] DataFrames v1.7.0
  [0c46a032] DifferentialEquations v7.14.0
  [5789e2e9] FileIO v1.16.3
  [961ee093] ModelingToolkit v9.41.0
  [8913a72c] NonlinearSolve v3.14.0
  [1dea7af3] OrdinaryDiffEq v6.89.0
  [91a5bcdd] Plots v1.40.8
  [0bca4576] SciMLBase v2.54.0
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [2f01184e] SparseArrays v1.10.0

What I mean is that solving the problem with discrete event by default solver:

# -------simplify the system--------
test_system_init = ODESystem(eqs, t,[],[]; systems = collect(values(merge(device_list, Bus_Dict))), name = :test_system_init, discrete_events = reflect)
test_system_init_simplify = structural_simplify(test_system_init, fully_determined = true)
# -------format the system to solve--------
prob = ODEProblem(test_system_init_simplify,[], (0,2),[]; guesses = u0_guess, initialization_eqs = ini_eqs, fully_determined = true, jac=true,sparse=true)
sol = solve(prob, tstops = [1.0,1.05])

The output:

apply fault at t = 1.0
apply fault at t = 1.0
ERROR: UndefRefError: access to undefined reference

if solving the problem with discrete event by Rodas4() solver:

# -------simplify the system--------
test_system_init = ODESystem(eqs, t,[],[]; systems = collect(values(merge(device_list, Bus_Dict))), name = :test_system_init, discrete_events = reflect)
test_system_init_simplify = structural_simplify(test_system_init, fully_determined = true)
# -------format the system to solve--------
prob = ODEProblem(test_system_init_simplify,[], (0,2),[]; guesses = u0_guess, initialization_eqs = ini_eqs, fully_determined = true, jac=true,sparse=true)
sol = solve(prob, Rodas4(), tstops = [1.0,1.05])

The output:

Warning: At t=1.0e-6, dt was forced below floating point epsilon 2.117582368135751e-22, and step error estimate = 6.659475516269642. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
retcode: Unstable
Interpolation: specialized 3rd order "free" stiffness-aware interpolation

[lixiao00093]

Is this a bug?
I am quite sure this can work in ModelingToolkit v9.38. in the following version.

using ModelingToolkit, DifferentialEquations, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D_t

# --------system models --------
@connector  AcBus begin
    @variables begin
        vm(t) = 1.0
        va(t) = 0.0
        pin(t) = 0.0, [connect = Flow]
        qin(t) = 0.0, [connect = Flow]
    end
end

@mtkmodel Zipload_component begin
    @parameters begin 
        u = 1,[description = "Enable Signal"]
        g
        b
    end
    @components begin
        bus = AcBus()
    end
    @equations begin
        bus.pin ~ u * (  g * bus.vm ^2)
        bus.qin ~ u * (- b * bus.vm ^2)
    end
end

@mtkmodel Piline_component begin
    @parameters begin
        u = 1 ,[description = "Enable Signal"]
        rft
        xft
        g0 = 0.0
        b0
        tap_ratio = 1.0
        gft =  rft/(rft^2 + xft^2)
        bft = -xft/(rft^2 + xft^2)
        m_real = tap_ratio 
        m2 = tap_ratio^2
    end

    @components begin
        bus_fr = AcBus()
        bus_to = AcBus()
    end

    @equations begin
        bus_fr.pin*m2 ~ u*( g0/2*bus_fr.vm^2 + gft*bus_fr.vm^2
                            - gft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
                            - bft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
        -bus_fr.qin*m2 ~ u*(b0/2*bus_fr.vm^2 + bft*bus_fr.vm^2
                            - bft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
                            + gft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
        
        bus_to.pin*m2 ~ u*( g0/2*m2*bus_to.vm^2 + gft*m2*bus_to.vm^2
                            - gft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
                            + bft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
        -bus_to.qin*m2 ~ u*(b0/2*m2*bus_to.vm^2 + bft*m2*bus_to.vm^2 
                            - bft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va) 
                            - gft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
    end
end

@mtkmodel Syn2_component begin
    @components begin
        bus = AcBus()
    end

    @variables begin
        delta(t), [description = "rotor angle"]
        omega(t) = 1.0, [description = "rotor speed"]

        vf(t)
        tm(t)
        p(t)
        q(t)

        _Id(t) = 0.0
        _Iq(t) = 0.0
    end

    @parameters begin
        time_flag 
        u = 1, [description = "Enable Signal"]
        fn
        H
        Damp
        ra
        xd1
        Kw
        Kp
        iM = u / (2 * H)
        c1 = ra / (ra^2 + xd1^2)
        c2 = xd1 / (ra^2 + xd1^2)
        c3 = c2
        vm0
        va0
    end

    @equations begin
        _Id ~ u*(-(bus.vm * sin(delta - bus.va) * c1) - (bus.vm * cos(delta - bus.va) * c3) + c3 * vf)
        _Iq ~ u*( (bus.vm * sin(delta - bus.va) * c2) - (bus.vm * cos(delta - bus.va) * c1) + c1 * vf)

        p ~ bus.vm * sin(delta - bus.va) * _Id + bus.vm * cos(delta - bus.va) * _Iq
        q ~ bus.vm * cos(delta - bus.va) * _Id - bus.vm * sin(delta - bus.va) * _Iq
        bus.pin ~ -p                                # for AcBus
        bus.qin ~ -q                                # for AcBus

        D_t(delta) ~ u * 2 * π * fn * (omega - 1)
        D_t(omega) ~ iM * (tm - p - ra * (_Id^2 + _Iq^2) - Damp * (omega - 1))
        D_t(tm) + (time_flag<0)* (bus.va-va0) ~ 0
        D_t(vf) + (time_flag<0)* (bus.vm-vm0) ~ 0
    end
end

@mtkmodel Fault_component begin
    @parameters begin
        u = 0,[description="fault status"]
        rf
        xf
        yg = rf/(rf^2 + xf^2)
        yb = -xf/(rf^2 + xf^2)
    end

    @components begin
        bus = AcBus()
    end

    @equations begin
        bus.pin ~  u * yg * bus.vm^2
        bus.qin ~ -u * yb * bus.vm^2
    end
end

# -------colection of instances of models--------
Bus_Dict = Dict{Symbol,ODESystem}(Symbol(1)=>AcBus(;name = Symbol("Bus 1")),Symbol(2)=>AcBus(;name = Symbol("Bus 2")))
device_list::Dict{Symbol,ODESystem} = Dict(
    :line => Piline_component(;name = :line, rft = 0.01,xft = 0.01,g0 = 0,b0 = 0.01,tap_ratio = 1.0),
    :syn => Syn2_component(;name = :syn1, time_flag = -1, fn = 50.0, H = 1.0, Damp = 0.8, ra = 0.01, xd1 = 0.1, Kw=1.0, Kp=1.0, va0 = 0.0 , vm0=1.0),
    :zipload => Zipload_component(;name = :zipload1, g = 1.0, b = 0.0),
    :fault => Fault_component(;name = :fault, rf = 0.0, xf = 0.1)
)

# -------connection of models--------
eqs = [ connect(Bus_Dict[Symbol(1)], device_list[:syn].bus),
        connect(Bus_Dict[Symbol(1)], device_list[:line].bus_fr),
        connect(Bus_Dict[Symbol(1)], device_list[:fault].bus),
        
        connect(Bus_Dict[Symbol(2)], device_list[:zipload].bus),
        connect(Bus_Dict[Symbol(2)], device_list[:line].bus_to)]

# -------Apply and clear fault in the models--------
u_temp_record = Dict{Int,Vector{Float64}}()
function affect1!(integ, u, p, ctx)
    if integ.ps[p.u] != ctx
        println("apply fault at t = ", integ.t)
        integ.ps[p.u] = ctx
        u_temp_record[1] = copy(integ.u[integ.differential_vars.==false])
    end
end
function affect2!(integ, u, p, ctx)
    if integ.ps[p.u] != ctx
        println("clear fault at t = ", integ.t)
        integ.ps[p.u] = ctx
        #offer the initial value of algebraic variables for the reinitialization
        integ.u[integ.differential_vars.==false] = u_temp_record[1]
    end
end

fault_function = [(t==1.0) => (affect1!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 1.0),
                 (t==1.05) => (affect2!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 0.0)]

# -------set the initial guessvalues of models--------
function set_ini_syn2(; ra, xd1)
    pg0 = 1.0
    qg0 = 0.0
    vm0 = 1.0
    va0 = 0.0

    V = vm0 .* exp.(va0 .* 1im)
    S = pg0 - qg0 .* 1im
    I = S ./ conj.(V)
    E = V + (ra + xd1 .* 1im) .* I
    delta = imag(log.(E ./ abs.(E) + 0im))
    delta0 = delta
    omega0 = 1.0

    # d- and q-axis voltages and currents
    vdq = V .* 1im .* exp.(-delta .* 1im)
    idq = I .* 1im .* exp.(-delta .* 1im)
    vd = real(vdq)
    vq = imag(vdq)
    Id = real(idq)
    Iq = imag(idq)
    # mechanical torques/powers
    tm0 = (vq + ra .* Iq) .* Iq + (vd + ra .* Id) .* Id
    vf0 =  vq + ra .* Iq + xd1 .* Id
    return delta0, omega0, tm0, vf0, pg0, qg0
end
u0_guess = Dict(zip([device_list[:syn].delta, device_list[:syn].omega, device_list[:syn].tm, device_list[:syn].vf, device_list[:syn].p, device_list[:syn].q], set_ini_syn2(ra = 0.01, xd1 = 0.1)))

# -------simplify the system--------
test_system_init = ODESystem(eqs, t,[],[]; systems = collect(values(merge(device_list, Bus_Dict))), name = :test_system_init, discrete_events = fault_function)
test_system_init_simplify = structural_simplify(test_system_init, fully_determined = true)


test_system_xy0 = SteadyStateProblem(test_system_init_simplify, u0_guess; jac = true, sparse = true)

sol = DifferentialEquations.solve(test_system_xy0)
xy0_sol_dict = Dict{Symbolics.Num, Any}()
for var in unknowns(test_system_xy0.f.sys)
    xy0_sol_dict[var] = sol[var]
end
for var in observed(test_system_xy0.f.sys)
    xy0_sol_dict[var.lhs] = sol[var.lhs]
end

defaults(test_system_init_simplify)[device_list[:syn].time_flag] = 1

# -------format the system to solve--------
prob = ODEProblem(test_system_init_simplify, xy0_sol_dict, (0,20); fully_determined = true, jac=true,sparse=true)
# sol = solve(prob, Rodas4(),  tstops = [1.0,1.05])
sol = solve(prob, dtmax = 0.01, tstops = [1.0,1.05])
plot(sol)

plot(sol.t,sol[device_list[:syn].delta])
plot(sol.t,sol[device_list[:syn].omega])
plot(sol.t,sol[device_list[:zipload].bus.vm])

But this way cannot work in the later version, because the way of initialisation has changed. But I do not know why I can not run my code after I modify that like the top side of this issue.

When I run the above codes, the output of Pkg.status()

  [6e4b80f9] BenchmarkTools v1.5.0
  [336ed68f] CSV v0.10.14
⌃ [a93c6f00] DataFrames v1.6.1
  [0c46a032] DifferentialEquations v7.14.0
  [5789e2e9] FileIO v1.16.3
⌃ [961ee093] ModelingToolkit v9.38.0
  [91a5bcdd] Plots v1.40.8
⌃ [0bca4576] SciMLBase v2.53.1
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging

[lixiao00093]

If nobody answers me, I would like to report that as a bug issue and close this one. This question is making me feel crazy.

[ChrisRackauckas]

There's just a lot of extra things in here so it will take some time to narrow it down before we can give a reasonable response. Is everything in that script actually needed to reproduce the issue? The plotting stuff for example seems unnecessary. If you can simplify the issue you'll get a quicker response, otherwise it will have to wait until I have the time to simplify the issue.

[lixiao00093]

There's just a lot of extra things in here so it will take some time to narrow it down before we can give a reasonable response. Is everything in that script actually needed to reproduce the issue? The plotting stuff for example seems unnecessary. If you can simplify the issue you'll get a quicker response, otherwise it will have to wait until I have the time to simplify the issue.

Thank you so much! Your response gives me the confidence to solve it.
Although I have simplified my case, I will continue to simplify this issue and update here. I know this will save a lot of time in solving this problem. Thank you very much!

[ChrisRackauckas]

Yes it will definitely help, since someone has to simplify it. If I have to do it then I'll need to wait until I have a good chunk of time open, but if it's simplified I can pretty say what's going on. My guess from a quick look is that it's the initialization things being fixed in #2922