ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

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

Open SLiemann opened this issue 2 years ago • 2 comments

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

SLiemann avatar Aug 16 '22 10:08 SLiemann

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.

YingboMa avatar Aug 16 '22 20:08 YingboMa

Was this handled?

ChrisRackauckas avatar Feb 22 '24 16:02 ChrisRackauckas