SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.39k stars 200 forks source link

Wrong solution for observed variables in case of parameter change (Acausal modeling) #1755

Open SLiemann opened 1 year ago

SLiemann commented 1 year ago

Hello,

below you can find a modified RC example, where I inserted an AC voltage source which decreases its voltage at t=0.04 s. You can see that the plotted solution is wrong if the voltage state of the source is an observed variable, while it is correct if it is not observed

In addition, the voltage of the resistor is also wrong. This could cause some problems if you want to use the solution for something else, e.g. calculating its power over time.

Is the solution of the observed variables calculated with the final set of parameters at the end of the simulation? And do you know how to implement this correctly without changing it to an observed variable? Or did I miss something?

Thank you very much for your help!

My MWE:

using ModelingToolkit, DifferentialEquations, Plots

function VoltageStep(;irr=false)
    @variables t
    @connector function Pin(;name)
        sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow] 
        ODESystem(Equation[], t, sts, []; name=name)
    end

    function Ground(;name)
        @named g = Pin()
        eqs = [g.v ~ 0]
        compose(ODESystem(eqs, t, [], []; name=name), g)
    end

    function OnePort(;name,irv = false, iri = false)
        @named p = Pin()
        @named n = Pin()
        @variables v(t)=1.0 [irreducible=irv] 
        @variables i(t)=1.0 [connect = Flow, irreducible=iri]
        eqs = [
            v ~ p.v - n.v
            0 ~ p.i + n.i
            i ~ p.i
            ]
        compose(ODESystem(eqs, t, [v,i], []; name=name), p, n)
    end

    function Resistor(;name, R = 1.0)
        @named oneport = OnePort()
        @unpack v, i = oneport
        ps = @parameters R=R
        eqs = [
            v ~ i * R
            ]
        extend(ODESystem(eqs, t, [], ps; name=name), oneport)
    end

    function Capacitor(;name, C = 1.0)
        @named oneport = OnePort()
        @unpack v, i = oneport
        ps = @parameters C=C
        D = Differential(t)
        eqs = [
            D(v) ~ i / C
            ]
        extend(ODESystem(eqs, t, [], ps; name=name), oneport)
    end

    function ACStepVoltage(;name, V = 1.0, freq = 1.0, phase = 0.0,irreducible = false)
        @named oneport = OnePort(irv=irreducible)
        @unpack v = oneport
        ps = @parameters V=V freq=freq phase=phase
        eqs = [
                v ~ V*sin(2*pi*freq*t+phase)
            ]
        extend(ODESystem(eqs, t, [], ps; name=name), oneport)
    end

    @named AC = ACStepVoltage(V=230*sqrt(2),freq = 50.0,irreducible=irr) #voltage state is either observed or not
    @named R = Resistor(R=10)
    @named C = Capacitor(C=1e-3)
    @named ground = Ground()

    rc_eqs = [
        connect(AC.p,R.p)
        connect(R.n,C.p)
        connect(C.n,ground.g,AC.n)
    ]

    step = [0.1] => [AC.V ~ 0.3*230*sqrt(2)]

    @named _rc_model = ODESystem(rc_eqs, t,[],[]; systems = [AC,R,C,ground],discrete_events = [step])
    sys = structural_simplify(_rc_model)

    u0 = zeros(length(states(sys)))

    prob = ODEProblem(sys, u0, (0.0, 0.2));
    sol = solve(prob,Rodas4(autodiff=false),dtmax = 1e-4,tstops=[0.1]); 
    plot!(sol,vars=[AC.v,R.v,C.v],layout=(1,3))
end

VoltageStep(irr=false) #not correct
VoltageStep(irr=true)  #correct

MTK_bug_observed

(@v1.7) pkg> st
      Status `C:\Users\liemann\.julia\environments\v1.7\Project.toml`
  [0c46a032] DifferentialEquations v7.2.0
  [615f187c] IfElse v0.1.1
  [961ee093] ModelingToolkit v8.18.9
  [91a5bcdd] Plots v1.31.7
YingboMa commented 1 year ago

Thanks for reporting the bug. The event system is still under development. This bug happens because the event system doesn't communicate the event information to the structural transformation code.

ChrisRackauckas commented 5 months ago

Was this handled?