Errors in affect codegen

[baggepinnen]

The following system simulates correctly without the event. The event affect function of the event is empty, so I would have expected it to have no effect.

using ModelingToolkit, OrdinaryDiffEq
import ModelingToolkit: t_nounits as t, D_nounits as D

@component function ContactForce2(; name, surface=nothing)
    vars = @variables begin
        q(t) = 1
        v(t) = 0
        f(t)
        h(t)
        hd(t)
        hdd(t)
        (λ(t)=0), [irreducible=true]
    end
    pars = @parameters begin
        contact::Int = 0    # discrete.time state variable
        a0 = 100
        a1 = 20
        a2 = 1e6
    end

    equations = [
        0 ~ ifelse((contact == 1), hdd + a1*hd + a0*h + a2*h^3, λ)
        f ~ contact*λ
        D(q) ~ v

        1 * D(v) ~ -1 * 9.81 + f

        h ~ q
        hd ~ D(h)
        hdd ~ D(hd)
    ]

    function affect!(integ, u, p, _)
    end
    continuous_events = [h ~ 0] => (affect!, [h], [], [], nothing)

    ODESystem(equations, t, vars, pars; name, continuous_events)
end
@named model = ContactForce2()
model = complete(model)
ssys = structural_simplify(model)
prob = ODEProblem(ssys, [], (0.0, 5.0))
sol = solve(prob, Rodas4())

using Plots
plot(sol, layout=5)

without event an object is falling freely due to gravity

with empty event

[baggepinnen]

The problem appears related to the equation

0 ~ ifelse((contact == 1), hdd + a1*hd + a0*h + a2*h^3, λ)

Since contact === 0 in this simulation, this equation always reads 0 ~ λ. If I use 0 ~ λ as equation instead, it works as expected:

using ModelingToolkit, OrdinaryDiffEq
import ModelingToolkit: t_nounits as t, D_nounits as D

@component function ContactForce2(; name, surface=nothing)
    vars = @variables begin
        q(t) = 1
        v(t) = 0
        f(t)
        (h(t) = 1), [irreducible=true]
        hd(t)
        hdd(t)
        (λ(t)=0), [irreducible=true]
    end
    pars = @parameters begin
        contact::Int = 0    # discrete.time state variable
        a0 = 100
        a1 = 20
        a2 = 1e6
    end

    equations = [
        0 ~ λ# ifelse((contact == 1), hdd + a1*hd + a0*h + a2*h^3, λ)
        f ~ contact*λ
        D(q) ~ v

        1 * D(v) ~ -1 * 9.81 + f

        h ~ q
        hd ~ D(h)
        hdd ~ D(hd)
    ]

    function affect!(integ, u, p, _)
    end
    continuous_events = [h ~ 0] => (affect!, [h], [], [], nothing)

    ODESystem(equations, t, vars, pars; name, continuous_events)
end
@named model = ContactForce2()
model = complete(model)
ssys = structural_simplify(model)
prob = ODEProblem(ssys, [], (0.0, 5.0))
sol = solve(prob, Rodas4())

using Plots
plot(sol, layout=5)

[baggepinnen]

Even more weird stuff going on. If I set dtmax to something small, the h coordinate is reset to it's initial value each time the empty event triggers

using ModelingToolkit, OrdinaryDiffEq
import ModelingToolkit: t_nounits as t, D_nounits as D

@component function ContactForce2(; name, surface=nothing)
    vars = @variables begin
        (q(t) = 1), [irreducible=true]
        v(t) = 0
        (f(t)=0), [irreducible=true]
        (h(t) = 1), [irreducible=true]
        (hd(t) = 0)#, [irreducible=true]
        hdd(t)#, [irreducible=true]
        (λ(t)=0), [irreducible=true]
    end
    pars = @parameters begin
        contact::Int = 0    # discrete.time state variable
        a0 = 100
        a1 = 20
        a2 = 0*1e6
    end

    equations = [
        0 ~ ifelse((contact == 1), hdd + a1*hd + a0*h + a2*h^3, λ)
        f ~ contact*λ
        D(q) ~ v

        1 * D(v) ~ -1 * 9.81 + f

        h ~ q
        hd ~ D(h)
        hdd ~ D(hd)
    ]

    function affect!(integ, u, p, _)
    end
    continuous_events = [h ~ 0] => (affect!, [h], [], [], nothing)

    ODESystem(equations, t, vars, pars; name, continuous_events)
end


@named model = ContactForce2()
model = complete(model)
ssys = structural_simplify(model)
prob = ODEProblem(ssys, [], (0.0, 2.0))
sol = solve(prob, Rodas4(), dtmax=0.001)

using Plots
plot(sol, layout=6, size=(1000, 1000))

[ChrisRackauckas]

@BenChung can you take this one?

[BenChung]

This is the bug we talked about on Slack where we reinitialize the system after the affect fires using the initialization system for the initial condition, rather than the "running" condition. I probably can take it - eventually - but I don't know enough about the initialization system to really be able to fix it properly.

[ChrisRackauckas]

Oh I see. This needs to be handled through what I had mentioned with the ImplicitDiscreteSystem thing if we're to make it robust.

[BenChung]

Yeah. My feeling - for now - is that it's better to not run initialization (or error if there is an initialization) after an affect so that the user doesn't get surprised by it and takes on the responsibility for ensuring that the requirements of a DAE are satisfied, and then we need to revisit this problem once we have better handling for system "live" initialization.

[ChrisRackauckas]

That's simply impossible.

[baggepinnen]

Why is running the initialization system required at all? If there are no algebraic equations, the state cannot be in feasible. If there are algebraic equations, the DAE solver will have to find a root to progress anyways and will solve the problem then?

The initialization system would be the incorrect system to solve anyways, parameters might have changed and equations for initial condition should no longer be included, it would be a completely different system.

[ChrisRackauckas]

You always have to run initialization after any affect! because otherwise you cannot guarantee a consistent state, and if you don't have a consistent state then you cannot take a step of the integrator.

The initialization system would be the incorrect system to solve anyways, parameters might have changed and equations for initial condition should no longer be included, it would be a completely different system.

That's the point of the implicit discrete form.

[HKruenaegel]

I have the same probleme with a discrete Callback, states are set to initial value every time the callback is fired.

What can I do to get the same behaviour as with MTK v8.73?

[HKruenaegel]

I have a model which is developed with MTK v8.72.0. With v8.75.0 it becomes slow an interrupts without in error after 3 seconds. With v9.19.0 it becomes much slower and interrupts with error at the same time. With >v9.36 it has the Problem described above.

I will try to shrink down the system for a example that I can share. But it seems that it's something with the descrete callback. Are there similar known Issues yet?

[BenChung]

@HKruenaegel We should be fixing the incorrect initialization system issue that causes the system to go to all-0s after an affect fires soon(^tm) when my development branch gets merged. I haven't had an opportunity to look into performance, though, so can't comment on that.

[HKruenaegel]

I tried again to run my model with MTKv9.30.0. If I give all states in u0 and the guesses are empty the model works, but as mentioned before it is 30 times slower than with MTKv8.73.2.

Is that a general observation with MTK9, @ChrisRackauckas?

[HKruenaegel]

Here is another example for the problem initially mentioned in this issue. This is a system which represents this circuit:

The resistor Rd is a ideal diode switch, therfore a Callback is defined, which sets the resistance to 1e4 if the voltage over Rd is <0 and to 1e-3 if it is >0.

The result of the simulation gives:

The current through the inductor is correct, also the voltage of the capacitor v3 during conducting state of the diode. But when the callback for switching off the diode is fired the capacitor voltage v3 is set to zero, which is not correct.

For this example MTK v.9.40.0 is used.

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

@parameters Rd Rsw C1 L1 V1 Rl
@variables I_L1(t) I_V1(t) v1(t) v2(t) [irreducible = true] v3(t) [irreducible = true]

eqs =         [I_L1 + I_V1~ 0
-I_L1 + v2/Rsw + v2/Rd - v3/Rd~ 0
C1*D(v3) + v3/Rl - v2/Rd + v3/Rd~ 0
v1~ 10*sin(2*pi*50*t)
-D(I_L1)*L1 + v1 - v2~ 0]

function D_on(i, u, p,ctx)
    i.ps[p.Rd]=1e-3
end
function D_off(i, u, p,ctx)
    i.ps[p.Rd]=1e3
end
c = ModelingToolkit.SymbolicContinuousCallback(
    [v2-v3~0], (D_on, [v2,v3], [Rd], [], nothing);affect_neg =(D_off, [v2,v3], [Rd], [], nothing),
    rootfind = SciMLBase.LeftRootFind)

@mtkbuild pend = ODESystem(eqs, t;continuous_events=c)

u0 = [I_L1=>0.0,
v3=>0.0,
]

g=[v1=>0.0,
v2=>0.0,
I_V1=>0.0]
p = [        Rd => 0.001,
Rsw=>1e6,
C1=>0.01,
L1=>0.001,
V1=>10.0,
Rl=>20.0]
prob = ODEProblem(pend, u0, (0.0, 1.5), p, guesses = g)

@time sol = solve(prob, Rodas5P(),dtmax=1e-4)

[HKruenaegel]

Ok, the trick is to use the initializealg NoInit(). Then there is no problem with states becoming initial and also the the simulation runs much faster.

@time sol = solve(prob, Rodas5P(),dtmax=1e-4;initializealg = NoInit())

[ChrisRackauckas]

This is now handled by the CheckInit default.

[HKruenaegel]

The code I posted earlier produces the same behaviour, I get the same wrong result.

[baggepinnen]

Yeah, the code I posted does not work either

[ChrisRackauckas]

@BenChung your change to checkinit default was already in?

[BenChung]

Should have in #3144. I'll look at this tomorrow; also working on getting the changes to initialize and finalize packaged out independently of ImperativeAffect.

[HKruenaegel]

Additionally I can not use the workaround with initializealg = NoInit():

@time sol = solve(prob, Rodas5P();initializealg = NoInit())

I get the error:

ERROR: OrdinaryDiffEqCore.CheckInitFailureError(89.33829573052752, 1.0e-6)

[ChrisRackauckas]

That means it's working. It's stopping you from getting a wrong result. In particular, after the callback the algebraic conditions are checked and you have a residual of 89.33. You previously had NoInit() which disabled this check and let the integration proceed, which can introduce unbounded error in the result which is why we do not recommend using that.

You can set NoInit in the callback to force it to skip this check, but I highly do not recommend that (and will probably regret mentioning that this is even possible) because that will allow incorrect / non error checked solutions through. At this time, the recommendation is to ensure that the affect results in a consistent system.

[BenChung]

I tried NoInit with Fredrik's example and it exhibits exactly the same behavior that it does with CheckInit so there's something else happening here. I think that either the underlying solver is still doing the init process beyond what's called for with CheckInit/NoInit or there's some other influence (potentially something to do with root finding?). I'm currently trying to replicate this without MTK.

[ChrisRackauckas]

To clarify, there are two inits. There is an init at the start of the DAE integration, that's the normal reinitialization. If you don't override it, sol = solve(prob, Rodas4()), then you're good. There's also a reinitialization after each callback. That one should now default to CheckInit(), which will check whether the algebraic constraints are satisfied before continuing the integration, and error if they are not.

It's that part that is and should be erroring in @HKruenaegel 's example because the change to Rd means -I_L1 + v2/Rsw + v2/Rd - v3/Rd~ 0 is no longer true, so you have to change some other values in order for that equation to be satisfied. A BrownBasicInit() is likely a good idea for this case if changing v2 to compensate for the parameter change is the expected behavior. Since it's impossible to guess the expected behavior, we simply error there.

[HKruenaegel]

Ok, so you mean I should use:

@time sol = solve(prob, Rodas5P();initializealg = BrownBasicInit())

But that gives me the error:

ERROR: UndefVarError: `BrownBasicInit` not defined

Also

@time sol = solve(prob, Rodas5P();initializealg = SciMLBase.BrownBasicInit())

fails.

[ChrisRackauckas]

@BenChung is this solved now?

[BenChung]

Should be fixed now; #3197 just adds tests that the examples given all work now.

[HKruenaegel]

The example given above is now running, but if I add a discrete event to change the Resistance of Rsw, I get the error again.

ERROR: CheckInit specified but initialization not satisfied. normresid = 56692.09744593064 > abstol = 1.0e-6

@ChrisRackauckas You mentioned that BrownBasicInit() could be the better alg for reinitialization, but how can I change the reinitialization alg?

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

@parameters Rd Rsw C1 L1 V1 Rl
@variables I_L1(t) I_V1(t) v1(t) v2(t) [irreducible = true] v3(t) [irreducible = true]

eqs =         [I_L1 + I_V1~ 0
-I_L1 + v2/Rsw + v2/Rd - v3/Rd~ 0
C1*D(v3) + v3/Rl - v2/Rd + v3/Rd~ 0
v1~ 10*sin(2*pi*50*t)
-D(I_L1)*L1 + v1 - v2~ 0]

function D_on(i, u, p,ctx)
    i.ps[p.Rd]=1e-3
end
function D_off(i, u, p,ctx)
    i.ps[p.Rd]=1e3
end
c = ModelingToolkit.SymbolicContinuousCallback(
    [v2-v3~0], (D_on, [v2,v3], [Rd], [], nothing);affect_neg =(D_off, [v2,v3], [Rd], [], nothing),
    rootfind = SciMLBase.LeftRootFind)

function bb_affect!(integ, u, p, ctx)
    if integ.ps[p.Rsw]==1e6
        integ.ps[p.Rsw]=1e-3
    else
        integ.ps[p.Rsw]=1e6
    end
end
d= 0.00016666666666666666 => (bb_affect!, [],[Rsw], [], nothing)

@mtkbuild pend = ODESystem(eqs, t;continuous_events=c,discrete_events=d)

u0 = [I_L1=>0.0,
v3=>100.0,
]

g=[v1=>0.0,
v2=>0.0,
I_V1=>0.0]
p = [        Rd => 0.001,
Rsw=>1e6,
C1=>0.01,
L1=>0.001,
V1=>100.0,
Rl=>1.0]
prob = ODEProblem(pend, u0, (0.0, 1), p, guesses = g)

@time sol = solve(prob, Rodas5P(),dtmax=1e-4)

[HKruenaegel]

Ok, I think it is answered in 3254.

[ChrisRackauckas]

Yes indeed that's what's going on here with the full details. I hope we're done with the v10 by the end of 2024. Basically, we have a lot of callback issues right now, but at least they boil down to the same fundamental problem, which is how to do the correct reinitialization after a callback, and so we just need to do the breaking change that makes this clearer / more explicit and then I think most callback issues will close at the same time. It's frustrating because that means there hasn't been the image of much progress this month, but because it's breaking it wil effectively all land at once.