Hi there, I’m new to Julia and modelingtoolkit, and I am trying to use it to implement a 2 element windkessel model consisting of a capacitor and resistor in parallel
My code is as follows:
using ModelingToolkit
using Plots
using DifferentialEquations
@variables t
D = Differential(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
eqs = [ D(v) ~ i / C ]
extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end
@component function DrivenCurrent(; name, I=1.0, fun)
@named oneport = OnePort()
@unpack i = oneport
ps = @parameters I = I
eqs = [
i ~ I * fun(t)
]
extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end
@component function WK2(; name, Rp=1.0, C=1.0)
@named in = Pin()
@named out = Pin()
sts = @variables Δv(t) = 0.0 i(t) = 0.0
# Parameters are inherited from subcomponents
ps = []
# These are the components the subsystem is made of:
@named Rp = Resistor(R=Rp)
@named C = Capacitor(C=C)
eqs = [
Δv ~ out.v - in.v
0 ~ in.i + out.i
i ~ in.i
connect(in, Rp.p, C.p)
connect(Rp.n, C.n, out)
]
# and finaly compose the system
compose(ODESystem(eqs, t, sts, ps; name=name), in, out, Rp, C)
end
@mtkmodel WK2model begin
@components begin
wk2 = WK2(Rp = 0.9, C = 1.1)
source = DrivenCurrent(I = 1.0,fun=sin)
ground = Ground()
end
@equations begin
connect(source.n, wk2.in)
connect(source.p, wk2.out, ground.g)
end
end
@mtkbuild wk2model = WK2model()
u0 = [
wk2model.wk2.in.v => 0.0
]
prob = ODEProblem(wk2model, u0, (0, 10.0))
sol_2_elem = solve(prob, RK4())
plot(sol_2_elem)
Which yields the following plot:
The ODE has been solved for the voltage across the capacitor as is seen in the graph (wk2.C.v(t)), however I would like it to have been solved for the voltage across the whole system (something like wk2.in.v(t)). I know that in this case they are one and the same because the capacitor and resistor are in parallel, however I would like to figure this out as I am hoping to build more complicated circuits in which this is not the case.
Is it possible to specify which voltage you wish to solve for when solving the ODE?
Thank you!