Help with a variation of the electrical components' example

Specific Domains Modelling & Simulations
,
1

Hello,

I have been doing my first tests with ModelingToolkit and when I tried to modify the acausal components example by adding other types of voltage source I had some problems that I don’t quite understand.

In the first example (a modification of the example in the link where I have added another voltage source), the problem has 4 states and 3 equations, so I can’ t solve it.

using ModelingToolkit, DifferentialEquations
using IfElse


@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)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; 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 ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function SourceVoltage(;name, V=1.0, freq = 20.0, fun)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters Vr,f
    eqs = [
           v ~  Vr * fun(t,f)
          ]
    extend(ODESystem(eqs, t, [], ps; defaults = Dict( Vr => V, f => freq), name=name), oneport)
end

sine(t,f) = sin(2pi * f * t)
@register sine(t,f::Number)
@named source = SourceVoltage(V = 1.0, freq = 2.5, fun = sine)
@named Ro = Resistor(R = 1.0)
@named C1 = Capacitor(C = 0.1)
@named G = Ground()

circ_eqs = [
	connect_pins(source.p,Ro.p,C1.p)
    connect_pins(source.n,Ro.n,G.g,C1.n)
         ]

@named modelCirc =  ODESystem(circ_eqs, t)
@named modelo = compose(modelCirc, [source,G,Ro,C1])

sys = structural_simplify(modelo)

Here you can see the output:

Model modelo with 3 equations
States (4):
  C1₊v(t) [defaults to 1.0]
  sine(t, source₊f)
  C1₊i(t) [defaults to 1.0]
  Ro₊i(t) [defaults to 1.0]
Parameters (4):
  source₊Vr [defaults to 1.0]
  source₊f [defaults to 2.5]
  Ro₊R [defaults to 1.0]
  C1₊C [defaults to 0.1]
Incidence matrix:
 ⋅  ⋅  ×  ⋅  ⋅  ×
 ×  ⋅  ⋅  ⋅  ×  ⋅
 ⋅  ×  ⋅  ⋅  ×  ⋅


equations(sys)
3-element Vector{Equation}:
 0 ~ Ro₊R*Ro₊i(t) - C1₊v(t)
 Differential(t)(C1₊v(t)) ~ C1₊i(t) / C1₊C
 0 ~ source₊Vr*Differential(t)(sine(t, source₊f)) + (-C1₊i(t)) / C1₊C

If instead of placing the three elements in parallel I try to solve a circuit with either an input resistance (CASE 1) or an RC circuit (CASE 2), the number of equations and states match. You can uncomment the CASE 1 or CASE 2 block to see the result.

function SourceVoltage(;name, V=1.0, freq = 20.0, fun)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters Vr,f
    eqs = [
           v ~  Vr * fun(t,f)
          ]
    extend(ODESystem(eqs, t, [], ps; defaults = Dict( Vr => V, f => freq), name=name), oneport)
end

sine(t,f) = sin(2pi * f * t)
@register sine(t,f::Number)
@named source = SourceVoltage(V = 1.0, freq = 2.5, fun = sine)
@named Ro = Resistor(R = 1.0)
@named C1 = Capacitor(C = 0.1)
@named G = Ground()
@named Ri = Resistor(R = 1.0)

#= CASE 1
circ_eqs = [
	connect_pins(source.p,Ri.p)
	connect_pins(Ri.n,Ro.p,C1.p)
    connect_pins(source.n,Ro.n,G.g,C1.n)
         ]

@named modelCirc =  ODESystem(circ_eqs, t)
@named modelo = compose(modelCirc, [source,G,Ro,C1,Ri])
=#


#= CASE 2
circ_eqs = [
	connect_pins(source.p,Ro.p)
    connect_pins(Ro.n,C1.p)
    connect_pins(C1.n, source.n, G.g)
         ]

@named modelCirc =  ODESystem(circ_eqs, t)
@named modelo = compose(modelCirc, [source,G,Ro,C1])
=#

sys = structural_simplify(modelo)
u0 = [source.v => 0.0
	  C1.v => 0.0]
prob = ODAEProblem(sys,u0, (0.0,5.0))
sol = solve(prob, Tsit5())
plot(sol, vars = [source.v, C1.i])

Returning to my initial circuit, if I define the voltage source as follows, I’m getting the expected results.

#Option 1
function SourceVoltage(;name, V=1.0, freq = 2.5)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters Vr,f
    D = Differential(t)
    eqs = [
           D(v) ~ 2*pi*f*Vr*cos(2*pi*f*t)
          ]
    extend(ODESystem(eqs, t, [], ps; defaults = Dict( Vr => V, f => freq), name=name), oneport)
end

@named source = SourceVoltage(V = 1.0, freq = 2.5)

#Option 2 (Works)
#=
function SourceVoltage(;name, V=1.0, freq = 2.5)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters Vr,f
    eqs = [
           v ~ Vr * sin(2*pi*f*t)
          ]
    extend(ODESystem(eqs, t, [], ps; defaults = Dict( Vr => V, f => freq), name=name), oneport)
end

@named source = SourceVoltage(V = 1.0, freq = 2.5)
=#

circ_eqs = [
	connect_pins(source.p,Ro.p,C1.p)
    connect_pins(source.n,Ro.n,G.g,C1.n)
         ]

I know that I am doing something wrong but I am unable to understand why an extra state appears in some cases.

Second example

The solved circuit is similar to the previous case. It consists of a voltage source, a resistor and a capacitor in parallel. As you can see, in the first case the time vector has only 7 elements and the result is always zero. If I add a resistor next to the voltage source, the problem is solved correctly. What am I missing?

using ModelingToolkit, DifferentialEquations
using IfElse

@variables t

function PWM_Voltage(;name, DutyCycle = 0.5, freq = 2.5)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters DC, f
    eqs = [
            v ~  IfElse.ifelse(mod(t,1/f) < DC*1/f,1.0,0.0)
          ]
    extend(ODESystem(eqs, t, [], ps; defaults = Dict( DC => DutyCycle, f => freq), name=name), oneport)
end

@named source = PWM_Voltage(freq = 2.5)
@named Ro = Resistor(R = 1.0)
@named C1 = Capacitor(C = 0.1)
@named G = Ground()

circ_eqs = [
	connect_pins(source.p,Ro.p,C1.p)
    connect_pins(source.n,Ro.n,G.g,C1.n)
         ]

@named modelCirc =  ODESystem(circ_eqs, t)
@named modelo = compose(modelCirc, [source,G,Ro,C1])
sys = structural_simplify(modelo)
u0 = [
	  C1.v => 0.0
	  source.v => 0.0
	 ]
prob = ODAEProblem(sys,u0,(0.0, 5.0))
sol = solve(prob, Tsit5())

#=
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 7-element Vector{Float64}:
 0.0
 9.999999999999999e-5
 0.0010999999999999998
 0.011099999999999997
 0.11109999999999996
 1.1110999999999995
 5.0
u: 7-element Vector{Vector{Float64}}:
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
=#

plot(sol)


@named Ri = Resistor(R = 1.0)
eq_circ2 = [
	connect_pins(source.p,Ri.p)
	connect_pins(Ri.n,Ro.p,C1.p)
    connect_pins(source.n,G.g,C1.n,Ro.n)
         ]

@named modelCirc_c2 =  ODESystem(eq_circ2, t)
@named modelo2 = compose(modelCirc_c2, [source,G,Ro,C1,Ri])
sys2 = structural_simplify(modelo2)
u0 = [
	C1.v => 0.0
	source.v => 0.0
]
prob2 = ODAEProblem(sys2,u0, (0.0,5.0))
sol2 = solve(prob2, Tsit5())

#=
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 232-element Vector{Float64}:
u: 232-element Vector{Vector{Float64}}:
=#

plot(sol2, vars = [source.v, C1.i,C1.v])

If possible, I would appreciate if someone has an explanation or can redirect me to a resource where I can learn how to set up the problems correctly.

Thanks.

2

In your first example, there is no need to register the driving function, remove this line:

@register sine(t,f::Number)

Registering a function indicates that the function should not be traced, so it is treated as opaque by all symbolic operations. When the derivative of this function is needed, it cannot be computed.

I’m not sure why adding the resistor fixes the second example, but I think the problem is that the never hits a change in the PWM voltage and steps right past all the edges. Try setting dtmax or tsteps Common Solver Options · DifferentialEquations.jl

3

Thanks for the reply. I tried not registering the function and it works! I was not sure if I should register it or not because in other simple circuits it worked in any case.

Looking more in detail at the second example I have seen that when simplifying the equations in the case that did not work for some reason the result of the ifelse is always 0.

#Without resistor
equations(sys)
3-element Vector{Equation}:
 0 ~ Ro₊R*Ro₊i(t) - C1₊v(t)
 Differential(t)(C1₊v(t)) ~ C1₊i(t) / C1₊C
 0 ~ (-C1₊i(t)) / C1₊C + IfElse.ifelse(mod(t, 1 / source₊f) < (source₊DC / source₊f), 0, 0)

#With resistor
equations(sys2)
3-element Vector{Equation}:
 0 ~ Ro₊R*Ro₊i(t) - C1₊v(t)
 Differential(t)(C1₊v(t)) ~ (Ri₊i(t) - Ro₊i(t)) / C1₊C
 0 ~ Ri₊R*Ri₊i(t) + C1₊v(t) - IfElse.ifelse(mod(t, 1 / source₊f) < (source₊DC / source₊f), 1.0, 0.0)
4

I think you are using an old version of MTK, I don’t see a connect_pin function in the current version. Here are the versions I used:

(jl_858WJR) pkg> st
      Status `/tmp/jl_858WJR/Project.toml`
  [0c46a032] DifferentialEquations v7.1.0
  [615f187c] IfElse v0.1.1
  [961ee093] ModelingToolkit v8.5.1
  [91a5bcdd] Plots v1.25.11

I see different equations, but similar behavior

@named source = PWM_Voltage(freq = 2.5)
@named Ro = Resistor(R = 1.0)
@named C1 = Capacitor(C = 0.1)
@named G = Ground()
eqs = [
       connect(source.p, Ro.p, C1.p),
       connect(source.n, Ro.n, C1.n, G.g),
      ]
@named circuit = ODESystem(eqs, t)
@named model = compose(circuit, [source, C1, G])
sys = structural_simplify(model)

julia> equations(sys)
2-element Vector{Equation}:
 Differential(t)(C1₊v(t)) ~ C1₊i(t) / C1₊C
 0 ~ (-C1₊i(t)) / C1₊C + IfElse.ifelse(mod(t, 1 / source₊f) < (source₊DC / source₊f), 0, 0)

I think it can be interpreted as the correct answer because the system is non-physical. You have a capacitor in parallel with a zero-impedance source that is switching state. This implies infinite current. Instead, I think the derivative of ifelse is defined to be zero everywhere, so the current is also zero. Adding a resistor in series with the source sets a finite impedance and the system can be solved.

5

I forgot to include the connect_pins function. Here is the code for it.

function connect_pins(ps...)
    eqs = [
            0 ~ sum(p->p.i, ps) #Current
          ]        
    for i = 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v) #Voltage
    end
    return eqs
end

In any case, I understand what you are saying about the input resistance. I didn’t think about it because when I tested it with SPICE i was getting results. Thanks for your help.