Skip to content

Weird results trying to simulate an RC circuit with AC power source #1152

@AayushSabharwal

Description

@AayushSabharwal

@ChrisRackauckas
I tried to extend the https://mtk.sciml.ai/stable/tutorials/acausal_components/ example, and tried playing around with the above circuit. Following is my code:

using ModelingToolkit, DifferentialEquations

@parameters t

# a point in the circuit with a given voltage and current
@connector function Pin(; name)
    sts = @variables v(t)=1.0 i(t)=1.0
    ODESystem(Equation[], t, sts, []; name)
end

# ground pin, where voltage is 0
function Ground(; name)
    @named g = Pin()
    eqs = [g.v ~ 0.]
    compose(ODESystem(eqs, t, [], []; name), g)
end

# abstraction for all 2-pin components
function OnePort(; name)
    # two pins
    @named p = Pin()
    @named n = Pin()
    # component has a voltage across it, and current through
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
        v ~ p.v - n.v    # KVL
        0. ~ p.i + n.i   # KCL
        i ~ p.i          # Current through component is current through +ve pin
    ]
    compose(ODESystem(eqs, t, sts, []; name), p, n)
end

function Resistor(; name, R = 1.0)
    # 2 pin component
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
        v ~ i * R
    ]
    extend(ODESystem(eqs, t, [], ps; 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), oneport)
end

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

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

# V = Vm sin(ωt+ϕ)
function ACVoltageSource(; name, Vm = 1.0, ω = 2π, ϕ = 0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters Vm=Vm ω=ω ϕ=ϕ
    eqs = [
        v ~ Vm * sin* t + ϕ)
    ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end

function ModelingToolkit.connect(::Type{Pin}, ps...)
    eqs = [
        0. ~ sum(p->p.i, ps)    # KCL
    ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end
    return eqs
end

@named acsrc = ACVoltageSource(Vm = 10, ω = 20π)
@named res = Resistor(R = 5.)
@named cap = Capacitor(C = 5e-6)
@named gnd = Ground()
eqs = [
           connect(acsrc.p, res.p)
           connect(res.n, cap.p)
           connect(cap.n, gnd.g, acsrc.n)
       ]
@named acmodel = compose(ODESystem(eqs, t), [res, cap, acsrc, gnd])
sys = structural_simplify(acmodel)
u0 = [
           cap.v => 0.0,
           res.i => 0.0,
       ]
prob = ODEProblem(sys, u0, (0., 50.))
sol = solve(prob, Rosenbrock23())

plot(sol) initially gave:
image

Which doesn't look right. However, I restarted my REPL and got this:
image

Which makes sense. I'm not entirely sure what went wrong initially.
There were some other intermediate plots with other forms of weird output with other solvers, but since I was working in a REPL, I didn't save them

Metadata

Metadata

Assignees

Labels

No labels
No labels

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions