Better RC circuit example

Author: YingboMa

test/components.jl

@@ -1,87 +1,98 @@
-using ModelingToolkit
 using Test
+using ModelingToolkit, StructuralTransformations, OrdinaryDiffEq
 
 # Basic electric components
-@parameters t
+const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
 function Pin(name)
     @variables v(t) i(t)
     ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0])
 end
 
 function Ground(name)
-    p = Pin(name)
-    eqs = [p.v ~ 0]
-    ODESystem(eqs, t, [], [], systems=[p], name=name)
+    g = Pin(:g)
+    eqs = [g.v ~ 0]
+    ODESystem(eqs, t, [], [], systems=[g], name=name)
 end
 
-function OnePort(name)
+function ConstantVoltage(name; V = 1.0)
+    val = V
     p = Pin(:p)
     n = Pin(:n)
-    @variables v(t) i(t)
+    @parameters V
     eqs = [
-           v ~ p.v - n.v
+           V ~ p.v - n.v
            0 ~ p.i + n.i
-           i ~ p.i
           ]
-    v, i, ODESystem(eqs, t, [v, i], [], systems=[p, n], name=name)
-end
-
-function ConstantVoltage(name; V = 1.0)
-    val = V
-    v, i, op = OnePort(name)
-    @parameters V
-    push!(op.eqs, v ~ V)
-    V = ModelingToolkit.value(V)
-    push!(op.ps, V)
-    op.default_p[V] = val
-    op
+    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
 end
 
 function Resistor(name; R = 1.0)
     val = R
-    v, i, op = OnePort(name)
+    p = Pin(:p)
+    n = Pin(:n)
+    @variables v(t)
     @parameters R
-    push!(op.eqs, v ~ R * i)
-    R = ModelingToolkit.value(R)
-    push!(op.ps, R)
-    op.default_p[R] = val
-    op
+    eqs = [
+           v ~ p.v - n.v
+           0 ~ p.i + n.i
+           v ~ p.i * R
+          ]
+    ODESystem(eqs, t, [v], [R], systems=[p, n], default_p=Dict(R => val), name=name)
 end
 
 function Capacitor(name; C = 1.0)
     val = C
-    v, i, op = OnePort(name)
+    p = Pin(:p)
+    n = Pin(:n)
+    @variables v(t)
     @parameters C
     D = Differential(t)
-    push!(op.eqs, D(v) ~ i / C)
-    C = ModelingToolkit.value(C)
-    push!(op.ps, C)
-    op.default_p[C] = val
-    op
+    eqs = [
+           v ~ p.v - n.v
+           0 ~ p.i + n.i
+           D(v) ~ p.i / C
+          ]
+    ODESystem(eqs, t, [v], [C], systems=[p, n], default_p=Dict(C => val), name=name)
 end
 
-function Inductor(name; L = 1.0)
-    val = L
-    v, i, op = OnePort(name)
-    @parameters L
-    D = Differential(t)
-    push!(op.eqs, D(i) ~ v / L)
-    L = ModelingToolkit.value(L)
-    push!(op.ps, L)
-    op.default_p[L] = val
-    op
+R = 1.0
+C = 1.0
+V = 1.0
+resistor = Resistor(:resistor, R=R)
+capacitor = Capacitor(:capacitor, C=C)
+source = ConstantVoltage(:source, V=V)
+ground = Ground(:ground)
+
+function connect(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
+rc_eqs = [
+          connect(source.p, resistor.p)
+          connect(resistor.n, capacitor.p)
+          connect(capacitor.n, source.n, ground.g)
+         ]
 
-R1 = Resistor(:R1)
-C1 = Capacitor(:C1)
-I1 = Inductor(:I1)
-@test isequal(R1.p.i, (@variables R1₊p₊i(t))[1])
-@parameters R
-@test ModelingToolkit.default_p(flatten(R1)) == Dict(R => 1.0)
-@variables n₊v(t) p₊v(t) n₊i(t) p₊i(t)
-@test ModelingToolkit.default_u0(flatten(R1)) == Dict(
-                                                      n₊v => 1.0,
-                                                      p₊v => 1.0,
-                                                      n₊i => 1.0,
-                                                      p₊i => 1.0,
-                                                     )
+rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground], name=:rc)
+sys = alias_elimination(rc_model)
+@test ModelingToolkit.default_p(sys) == Dict(
+                                             capacitor.C => 1.0,
+                                             source.V => 1.0,
+                                             resistor.R => 1.0,
+                                            )
+u0 = [
+      capacitor.v => 0.0
+      capacitor.p.i => 0.0
+      resistor.v => 0.0
+     ]
+prob = ODEProblem(sys, u0, (0, 10.0))
+sol = solve(prob, Rodas5())
+#using Plot
+#plot(sol)