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?