From 96eb450f917fb93546bf690061883db768cd6796 Mon Sep 17 00:00:00 2001 From: Tom Short Date: Fri, 2 Jan 2015 14:28:13 -0500 Subject: [PATCH 1/2] Changes for the equations macro and Equation[]. --- README.md | 76 ++-- doc/Possible-future-developments.md | 8 +- doc/README.md | 66 ++-- examples/basics/breaking_pendulum.jl | 38 +- examples/basics/breaking_pendulum_in_box.jl | 56 +-- examples/basics/circuit_complex.jl | 24 +- examples/basics/dc_motor_w_shaft.jl | 51 +-- examples/basics/half_wave_rectifiers.jl | 70 ++-- examples/basics/initial_conditions.jl | 24 +- examples/basics/vanderpol.jl | 6 +- examples/basics/vanderpol_with_events.jl | 26 +- examples/neural/CaT.jl | 66 ++-- examples/neural/adex.jl | 26 +- examples/neural/fhn.jl | 8 +- examples/neural/hh.jl | 22 +- examples/neural/hh_base.jl | 46 +-- examples/neural/hh_modular.jl | 46 +-- examples/neural/hr.jl | 10 +- examples/neural/iaf.jl | 16 +- examples/neural/iafrefr.jl | 40 +- examples/neural/iafsyn.jl | 68 ++-- examples/neural/izhfs.jl | 24 +- examples/neural/ml.jl | 8 +- examples/neural/mlsyn.jl | 26 +- examples/neural/netiafrefr.jl | 46 +-- examples/neural/pr.jl | 52 +-- examples/neural/theta.jl | 24 +- examples/run_all.jl | 2 +- examples/stdlib/m_blocks.jl | 48 +-- examples/stdlib/m_electrical.jl | 396 ++++++++++---------- examples/stdlib/m_heat_transfer.jl | 36 +- examples/stdlib/m_powersystems.jl | 44 +-- examples/stdlib/m_rotational.jl | 18 +- src/blocks.jl | 119 +++--- src/electrical.jl | 322 ++++++++-------- src/heat_transfer.jl | 52 +-- src/machines.jl | 242 ++++++------ src/main.jl | 8 +- src/powersystems.jl | 76 ++-- src/rotational.jl | 238 ++++++------ 40 files changed, 1295 insertions(+), 1279 deletions(-) diff --git a/README.md b/README.md index 5014bb7..11a052a 100644 --- a/README.md +++ b/README.md @@ -86,7 +86,7 @@ MExpr(*(##1243,+(##1243,1))) In a simulation, the unknowns are to be solved based on a set of equations. Equations are built from device models. -A device model is a function that returns a list of equations or +A device model is a function that returns a vector of equations or other devices that also return lists of equations. The equations each are assumed equal to zero. So, @@ -112,11 +112,11 @@ function Vanderpol() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Expressions of # regular variables are evaluated immediately. - { - # The -1.0 in der(x, -1.0) is the initial value for the derivative - der(x, -1.0) - ((1 - y^2) * x - y) # == 0 is assumed - der(y) - x - } + Equation[ + # The -1.0 in der(x, -1.0) is the initial value for the derivative + der(x, -1.0) - ((1 - y^2) * x - y) # == 0 is assumed + der(y) - x + ] end y = sim(Vanderpol(), 10.0) # Run the simulation to 10 seconds and return @@ -127,8 +127,26 @@ gplot(y) Here are the results: -![plot results](https://github.com/tshort/Sims.jl/blob/master/examples/vanderpol.png?raw=true "Van Der Pol results") +![plot results](https://github.com/tshort/Sims.jl/blob/master/examples/basics/vanderpol.png?raw=true "Van Der Pol results") +An `@equations` macro is provided to return `Equation[]` allowing for +the use of equals in equations, so the example above can be: + +``` julia +function Vanderpol() + y = Unknown(1.0, "y") + x = Unknown("x") + @equations begin + der(x, -1.0) = (1 - y^2) * x - y + der(y) = x + end +end + +y = sim(Vanderpol(), 10.0) # Run the simulation to 10 seconds and return + # the result as an array. +# plot the results with Gaston +gplot(y) +``` Electrical example ------------------ @@ -159,19 +177,19 @@ voltage. function Resistor(n1, n2, R::Real) i = Current() # This is simply an Unknown. v = Voltage() - { - Branch(n1, n2, v, i) - R * i - v # == 0 is implied - } + @equations begin + Branch(n1, n2, v, i) + R * i = v + end end function Capacitor(n1, n2, C::Real) i = Current() v = Voltage() - { - Branch(n1, n2, v, i) - C * der(v) - i - } + @equations begin + Branch(n1, n2, v, i) + C * der(v) = i + end end ``` @@ -188,13 +206,13 @@ function Circuit() n2 = Voltage("Output voltage") n3 = Voltage() g = 0.0 # A ground has zero volts; it's not an unknown. - { - SineVoltage(n1, g, 10.0, 60.0) - Resistor(n1, n2, 10.0) - Resistor(n2, g, 5.0) - SeriesProbe(n2, n3, "Capacitor current") - Capacitor(n3, g, 5.0e-3) - } + Equation[ + SineVoltage(n1, g, 10.0, 60.0) + Resistor(n1, n2, 10.0) + Resistor(n2, g, 5.0) + SeriesProbe(n2, n3, "Capacitor current") + Capacitor(n3, g, 5.0e-3) + ] end ckt = Circuit() @@ -203,7 +221,7 @@ gplot(ckt_y) ``` Here are the results: -![plot results](https://github.com/tshort/Sims.jl/blob/master/examples/circuit.png?raw=true "Circuit results") +![plot results](https://github.com/tshort/Sims.jl/blob/master/examples/basics/circuit.png?raw=true "Circuit results") Initialization and Solving Sets of Equations -------------------------------------------- @@ -216,10 +234,10 @@ equations: ```julia function test() @unknown x y - { - 2*x - y - exp(-x) - -x + 2*y - exp(-y) - } + @equations begin + 2*x - y = exp(-x) + -x + 2*y = exp(-y) + end end solution = solve(create_sim(test())) @@ -231,7 +249,7 @@ Hybrid Modeling and Structural Variability Sims supports basic hybrid modeling, including the ability to handle structural model changes. Consider the following example: -[Breaking pendulum](https://github.com/tshort/Sims.jl/blob/master/examples/breaking_pendulum_in_box.jl) +[Breaking pendulum](https://github.com/tshort/Sims.jl/blob/master/examples/basics/breaking_pendulum_in_box.jl) This model starts as a pendulum, then the wire breaks, and the ball goes into free fall. Sims handles this much like @@ -245,7 +263,7 @@ adjusted for the "bounce". Here is an animation of the results. Note that the actual animation was done in R, not Julia. -![plot results](https://github.com/tshort/Sims.jl/blob/master/examples/pendulum.gif?raw=true "Pendulum") +![plot results](https://github.com/tshort/Sims.jl/blob/master/examples/basics/pendulum.gif?raw=true "Pendulum") To Look Deeper -------------- diff --git a/doc/Possible-future-developments.md b/doc/Possible-future-developments.md index 8bbd4e4..809faf3 100644 --- a/doc/Possible-future-developments.md +++ b/doc/Possible-future-developments.md @@ -92,10 +92,10 @@ leading "tag". function Resistor(n1::ElectricalNode, n2::ElectricalNode, R::Signal) i = Current(compatible_values(n1, n2)) v = Voltage(value(n1) - value(n2)) - { - Branch(n1, n2, v, i) - R .* i - v # == 0 is implied - } + Equation[ + Branch(n1, n2, v, i) + R .* i - v # == 0 is implied + ] end # help info - returns a string in Markdown format diff --git a/doc/README.md b/doc/README.md index 90b06a2..ff9872f 100644 --- a/doc/README.md +++ b/doc/README.md @@ -88,10 +88,10 @@ function Vanderpol() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(x, -1.0) - ((1 - y^2) * x - y) # == 0 is assumed - der(y) - x - } + Equation[ + der(x, -1.0) - ((1 - y^2) * x - y) # == 0 is assumed + der(y) - x + ] end ``` @@ -108,10 +108,10 @@ unknowns and equations as shown below: function Capacitor(n1, n2, C::Real) i = Current() # Unknown #1 v = Voltage() # Unknown #2 - { - Branch(n1, n2, v, i) # Equation #1 - this returns n1 - n2 - v - C * der(v) - i # Equation #2 - } + Equation[ + Branch(n1, n2, v, i) # Equation #1 - this returns n1 - n2 - v + C * der(v) - i # Equation #2 + ] end ``` @@ -127,13 +127,13 @@ function Circuit() n2 = ElectricalNode("Output voltage") n3 = ElectricalNode() g = 0.0 # a ground has zero volts; it's not an Unknown. - { - VSource(n1, g, 10.0, 60.0) - Resistor(n1, n2, 10.0) - Resistor(n2, g, 5.0) - SeriesProbe(n2, n3, "Capacitor current") - Capacitor(n3, g, 5.0e-3) - } + Equation[ + VSource(n1, g, 10.0, 60.0) + Resistor(n1, n2, 10.0) + Resistor(n2, g, 5.0) + SeriesProbe(n2, n3, "Capacitor current") + Capacitor(n3, g, 5.0e-3) + ] end ``` @@ -212,13 +212,13 @@ function VSquare(n1, n2, V::Real, f::Real) i = Current() v = Voltage() v_mag = Discrete(V) - { - Branch(n1, n2, v, i) - v - v_mag - Event(sin(2 * pi * f * MTime), - {reinit(v_mag, V)}, # positive crossing - {reinit(v_mag, -V)}) # negative crossing - } + Equation[ + Branch(n1, n2, v, i) + v - v_mag + Event(sin(2 * pi * f * MTime), + Equation[reinit(v_mag, V)], # positive crossing + Equation[reinit(v_mag, -V)]) # negative crossing + ] end ``` @@ -238,12 +238,12 @@ function IdealDiode(n1, n2) v = Voltage() s = Unknown() # dummy variable openswitch = Discrete(false) # on/off state of diode - { - Branch(n1, n2, v, i) - BoolEvent(openswitch, -s) # openswitch becomes true when s goes negative - v - ifelse(openswitch, s, 0.0) - i - ifelse(openswitch, 0.0, s) - } + Equation[ + Branch(n1, n2, v, i) + BoolEvent(openswitch, -s) # openswitch becomes true when s goes negative + v - ifelse(openswitch, s, 0.0) + i - ifelse(openswitch, 0.0, s) + ] end ``` @@ -274,11 +274,11 @@ function BreakingPendulum() y = Unknown(-cos(pi/4), "y") vx = Unknown() vy = Unknown() - { - StructuralEvent(MTime - 5.0, # when time hits 5 sec, switch to FreeFall - Pendulum(x,y,vx,vy), - () -> FreeFall(x,y,vx,vy)) - } + Equation[ + StructuralEvent(MTime - 5.0, # when time hits 5 sec, switch to FreeFall + Pendulum(x,y,vx,vy), + () -> FreeFall(x,y,vx,vy)) + ] end ``` diff --git a/examples/basics/breaking_pendulum.jl b/examples/basics/breaking_pendulum.jl index d3d2f56..c278ee2 100644 --- a/examples/basics/breaking_pendulum.jl +++ b/examples/basics/breaking_pendulum.jl @@ -6,12 +6,12 @@ using Sims ######################################## function FreeFall(x,y,vx,vy) - { - vx - der(x) - vy - der(y) - der(vx) - der(vy) + 9.81 - } + @equations begin + der(x) = vx + der(y) = vy + der(vx) = 0.0 + der(vy) = -9.81 + end end function Pendulum(x,y,vx,vy) @@ -19,14 +19,14 @@ function Pendulum(x,y,vx,vy) phi0 = atan2(x.value, -y.value) phi = Unknown(phi0) phid = Unknown() - { - phid - der(phi) - vx - der(x) - vy - der(y) - x - len * sin(phi) - y + len * cos(phi) - der(phid) + 9.81 / len * sin(phi) - } + @equations begin + der(phi) = phid + der(x) = vx + der(y) = vy + x = len * sin(phi) + y = -len * cos(phi) + der(phid) = -9.81 / len * sin(phi) + end end function BreakingPendulum() @@ -34,11 +34,11 @@ function BreakingPendulum() y = Unknown(-cos(pi/4), "y") vx = Unknown() vy = Unknown() - { - StructuralEvent(MTime - 5.0, # when time hits 5 sec, switch to FreeFall - Pendulum(x,y,vx,vy), - () -> FreeFall(x,y,vx,vy)) - } + Equation[ + StructuralEvent(MTime - 5.0, # when time hits 5 sec, switch to FreeFall + Pendulum(x,y,vx,vy), + () -> FreeFall(x,y,vx,vy)) + ] end p = BreakingPendulum() diff --git a/examples/basics/breaking_pendulum_in_box.jl b/examples/basics/breaking_pendulum_in_box.jl index 1b8a95d..db9ec65 100644 --- a/examples/basics/breaking_pendulum_in_box.jl +++ b/examples/basics/breaking_pendulum_in_box.jl @@ -6,21 +6,21 @@ using Sims ############################################ function FreeFall(x,y,vx,vy) - { - vx - der(x) - vy - der(y) - der(vx) - der(vy) + 9.81 - Event(y + 1.01, # Bounce up if it hits the floor - {reinit(vy, 0.95 * abs(vy))}, - {reinit(vy, 0.95 * abs(vy))}) - Event(x + 1.01, # Left wall - {reinit(vx, 0.95 * abs(vx))}, - {reinit(vx, 0.95 * abs(vx))}) - Event(x - 1.01, # Right wall - {reinit(vx, -0.95 * abs(vx))}, - {reinit(vx, -0.95 * abs(vx))}) - } + @equations begin + der(x) = vx + der(y) = vy + der(vx) + der(vy) = -9.81 + Event(y + 1.01, # Bounce up if it hits the floor + Equation[reinit(vy, 0.95 * abs(vy))], + Equation[reinit(vy, 0.95 * abs(vy))]) + Event(x + 1.01, # Left wall + Equation[reinit(vx, 0.95 * abs(vx))], + Equation[reinit(vx, 0.95 * abs(vx))]) + Event(x - 1.01, # Right wall + Equation[reinit(vx, -0.95 * abs(vx))], + Equation[reinit(vx, -0.95 * abs(vx))]) + end end function Pendulum(x,y,vx,vy) @@ -28,14 +28,14 @@ function Pendulum(x,y,vx,vy) phi0 = atan2(x.value, -y.value) phi = Unknown(phi0) phid = Unknown() - { - phid - der(phi) - vx - der(x) - vy - der(y) - x - len * sin(phi) - y + len * cos(phi) - der(phid) + 9.81 / len * sin(phi) - } + @equations begin + der(phi) = phid + der(x) = vx + der(y) = vy + x = len * sin(phi) + y = -len * cos(phi) + der(phid) = -9.81 / len * sin(phi) + end end function BreakingPendulumInBox() @@ -43,11 +43,11 @@ function BreakingPendulumInBox() y = Unknown(-cos(pi/4), "y") vx = Unknown() vy = Unknown() - { - StructuralEvent(MTime - 1.8, # when time hits 1.8 sec, switch to FreeFall - Pendulum(x,y,vx,vy), - () -> FreeFall(x,y,vx,vy)) - } + Equation[ + StructuralEvent(MTime - 1.8, # when time hits 1.8 sec, switch to FreeFall + Pendulum(x,y,vx,vy), + () -> FreeFall(x,y,vx,vy)) + ] end p = BreakingPendulumInBox() diff --git a/examples/basics/circuit_complex.jl b/examples/basics/circuit_complex.jl index 24d6621..88e3e9e 100644 --- a/examples/basics/circuit_complex.jl +++ b/examples/basics/circuit_complex.jl @@ -31,29 +31,29 @@ typealias ComplexCurrent Unknown{UComplexCurrent} function RL(n1, n2, Z::Complex) i = Unknown(0.im, "RL i") v = Unknown(0.im) - { - Branch(n1, n2, v, i) - Z * i - v # Note that this doesn't include time variation (L di/dt) effects - } + @equations begin + Branch(n1, n2, v, i) + v = Z * i # Note that this doesn't include time variation (L di/dt) effects + end end function VCmplxSrc(n1, n2, V::Number) i = Unknown(0.im) v = Unknown(0.im) - { - Branch(n1, n2, v, i) - V - v - } + @equations begin + Branch(n1, n2, v, i) + v = V + end end # This is just a static circuit. Nothing is time varying. function CmplxCkt() n = ComplexElectricalNode(0.im, "node voltage") g = 0.0 - { - VCmplxSrc(n, g, 10.+2im) - RL(n, g, 3 + 4im) - } + Equation[ + VCmplxSrc(n, g, 10.+2im) + RL(n, g, 3 + 4im) + ] end diff --git a/examples/basics/dc_motor_w_shaft.jl b/examples/basics/dc_motor_w_shaft.jl index 1fe2685..b3304a0 100644 --- a/examples/basics/dc_motor_w_shaft.jl +++ b/examples/basics/dc_motor_w_shaft.jl @@ -24,13 +24,13 @@ function EMF(n1::ElectricalNode, n2::ElectricalNode, flange::Flange, k::Real) i = Current() v = Voltage() w = AngularVelocity() - { - Branch(n1, n2, i, v) - RefBranch(flange, tau) - w - der(flange) - v - k * w - tau - k * i - } + Equation[ + Branch(n1, n2, i, v) + RefBranch(flange, tau) + w - der(flange) + v - k * w + tau - k * i + ] end function DCMotor(flange::Flange) @@ -38,21 +38,21 @@ function DCMotor(flange::Flange) n2 = Voltage() n3 = Voltage() g = 0.0 - { - SignalVoltage(n1, g, 60.0) - Resistor(n1, n2, 100.0) - Inductor(n2, n3, 0.2) - EMF(n3, g, flange, 1.0) - } + Equation[ + SignalVoltage(n1, g, 60.0) + Resistor(n1, n2, 100.0) + Inductor(n2, n3, 0.2) + EMF(n3, g, flange, 1.0) + ] end function ShaftElement(flangeA::Flange, flangeB::Flange) r1 = Angle() - { - Spring(flangeA, r1, 8.0) - Damper(flangeA, r1, 1.5) - Inertia(r1, flangeB, 0.5) - } + Equation[ + Spring(flangeA, r1, 8.0) + Damper(flangeA, r1, 1.5) + Inertia(r1, flangeB, 0.5) + ] end function FlexibleShaft(flangeA::Flange, flangeB::Flange, n::Int) @@ -63,7 +63,7 @@ function FlexibleShaft(flangeA::Flange, flangeB::Flange, n::Int) end r[1] = flangeA r[end] = flangeB - result = {} + result = Equation[] for i in 1:(n - 1) push!(result, ShaftElement(r[i], r[i + 1])) end @@ -74,14 +74,15 @@ function MechSys() r1 = Angle("Source angle") r2 = Angle() r3 = Angle("End-of-shaft angle") - { - DCMotor(r1) - Inertia(r1, r2, 0.02) - FlexibleShaft(r2, r3, 5) - } + Equation[ + DCMotor(r1) + Inertia(r1, r2, 0.02) + FlexibleShaft(r2, r3, 5) + ] end -m = MechSys() +m = MechSys() # returns the hierarchical model +m_f = elaborate(m) # returns the flattened model m_yout = sim(m, 4.0) wplot(m_yout, "DcMotorWShaft.pdf") diff --git a/examples/basics/half_wave_rectifiers.jl b/examples/basics/half_wave_rectifiers.jl index 1dd870d..a543e22 100644 --- a/examples/basics/half_wave_rectifiers.jl +++ b/examples/basics/half_wave_rectifiers.jl @@ -6,10 +6,10 @@ import Sims.IdealDiode function VSource(n1::NumberOrUnknown{UVoltage}, n2::NumberOrUnknown{UVoltage}, V::Real, f::Real, ang::Real) i = Current() v = Voltage() - { - Branch(n1, n2, v, i) - v - V * sin(2 * pi * f * MTime + ang) - } + @equations begin + Branch(n1, n2, v, i) + v = V * sin(2 * pi * f * MTime + ang) + end end @@ -19,28 +19,32 @@ function IdealDiode(n1, n2) v = Voltage() s = Unknown() # dummy variable openswitch = Discrete(false) # on/off state of diode - { - Branch(n1, n2, v, i) - BoolEvent(openswitch, -s) # openswitch becomes true when s goes negative - v - ifelse(openswitch, s, 0.0) - i - ifelse(openswitch, 0.0, s) - ## v - s * ifelse(openswitch, 1.0, 1e-5) - ## i - s * ifelse(openswitch, 1e-5, 1.0) - } + @equations begin + Branch(n1, n2, v, i) + BoolEvent(openswitch, -s) # openswitch becomes true when s goes negative + v = ifelse(openswitch, s, 0.0) + i = ifelse(openswitch, 0.0, s) + ## v = s * ifelse(openswitch, 1.0, 1e-5) + ## i = s * ifelse(openswitch, 1e-5, 1.0) + end end function OpenDiode(n1, n2) v = Voltage() - StructuralEvent(v+0.0, # when V goes positive, this changes to a ClosedDiode - Branch(n1, n2, v, 0.0), - () -> ClosedDiode(n1, n2)) + Equation[ + StructuralEvent(v+0.0, # when V goes positive, this changes to a ClosedDiode + Branch(n1, n2, v, 0.0), + () -> ClosedDiode(n1, n2)) + ] end function ClosedDiode(n1, n2) i = Current() - StructuralEvent(-i, # when I goes negative, this changes to an OpenDiode - Branch(n1, n2, 0.0, i), - () -> OpenDiode(n1, n2)) + Equation[ + StructuralEvent(-i, # when I goes negative, this changes to an OpenDiode + Branch(n1, n2, 0.0, i), + () -> OpenDiode(n1, n2)) + ] end # Cellier, fig 9.27 @@ -49,13 +53,13 @@ function HalfWaveRectifier() n2 = Voltage("") nout = Voltage("Output voltage") g = 0.0 - { - VSource(nsrc, g, 1.0, 60.0, pi/32) - Resistor(nsrc, n2, 1.0) - IdealDiode(n2, nout) - Capacitor(nout, g, 0.001) - Resistor(nout, g, 50.0) - } + Equation[ + VSource(nsrc, g, 1.0, 60.0, pi/32) + Resistor(nsrc, n2, 1.0) + IdealDiode(n2, nout) + Capacitor(nout, g, 0.001) + Resistor(nout, g, 50.0) + ] end @@ -77,14 +81,14 @@ function StructuralHalfWaveRectifier() nout = Voltage("Output voltage") Vdiode = Unknown("Vdiode") # probe variable g = 0.0 - { - Vdiode - (n2 - nout) - VSource(nsrc, g, 1.0, 60.0, pi/32) - Resistor(nsrc, n2, 1.0) - ClosedDiode(n2, nout) - Capacitor(nout, g, 0.001) - Resistor(nout, g, 50.0) - } + @equations begin + Vdiode = n2 - nout + VSource(nsrc, g, 1.0, 60.0, pi/32) + Resistor(nsrc, n2, 1.0) + ClosedDiode(n2, nout) + Capacitor(nout, g, 0.001) + Resistor(nout, g, 50.0) + end end println("**** Structural Half Wave Rectifier ****") diff --git a/examples/basics/initial_conditions.jl b/examples/basics/initial_conditions.jl index f80f8f5..0289a63 100644 --- a/examples/basics/initial_conditions.jl +++ b/examples/basics/initial_conditions.jl @@ -12,10 +12,10 @@ using Sims function test() @unknown x y - { - 2*x - y - exp(-x) - -x + 2*y - exp(-y) - } + @equations begin + 2*x - y = exp(-x) + -x + 2*y = exp(-y) + end end f = elaborate(test()) sm = create_sim(f) @@ -23,10 +23,10 @@ res = inisolve(sm) function mkin() @unknown x(1.0) y(1.0) - { - x^2 + y^2 - 1.0 - y - x^2 - } + @equations begin + x^2 + y^2 = 1.0 + y = x^2 + end end sm = create_sim(mkin()) @@ -37,10 +37,10 @@ res = inisolve(sm) function fun() @unknown x y = Unknown(:y, 1.0, true) - { - x^2 + y^2 - 1.0 - der(y, 0.0) - x - } + @equations begin + x^2 + y^2 = 1.0 + der(y, 0.0) = x + end end ## sm = create_sim(fun()) diff --git a/examples/basics/vanderpol.jl b/examples/basics/vanderpol.jl index fae87bc..8999140 100644 --- a/examples/basics/vanderpol.jl +++ b/examples/basics/vanderpol.jl @@ -33,9 +33,9 @@ function Vanderpol() # Expressions with Unknowns are kept as expressions. Expressions of # regular variables are evaluated immediately (like normal). @equations begin - # The -1.0 in der(x, -1.0) is the initial value for the derivative - der(x, -1.0) - ((1 - y^2) * x - y) # == 0 is assumed - der(y) - x + # The -1.0 in der(x, -1.0) is the initial value for the derivative + der(x, -1.0) = (1 - y^2) * x - y + der(y) = x end end diff --git a/examples/basics/vanderpol_with_events.jl b/examples/basics/vanderpol_with_events.jl index 8d0a9f3..d4a0a1e 100644 --- a/examples/basics/vanderpol_with_events.jl +++ b/examples/basics/vanderpol_with_events.jl @@ -16,19 +16,19 @@ function SVanderpol() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Expressions of # regular variables are evaluated immediately (like normal). - { - # The -1.0 in der(x, -1.0) is the initial value for the derivative - der(x, -1.0) - (mu * (1 - y^2) * x - y) # == 0 is assumed - der(y) - x - mu_unk - mu - Event(sin(pi/2 * MTime), # Initiate an event every 2 sec. - { - reinit(mu, mu * 0.75) - }, - { - reinit(mu, mu * 1.8) - }) - } + @equations begin + # The -1.0 in der(x, -1.0) is the initial value for the derivative + der(x, -1.0) = mu * (1 - y^2) * x - y + der(y) = x + mu_unk = mu + Event(sin(pi/2 * MTime), # Initiate an event every 2 sec. + Equation[ + reinit(mu, mu * 0.75) + ], + Equation[ + reinit(mu, mu * 1.8) + ]) + end end v = SVanderpol() # returns the hierarchical model diff --git a/examples/neural/CaT.jl b/examples/neural/CaT.jl index 9e754cd..0b68b22 100644 --- a/examples/neural/CaT.jl +++ b/examples/neural/CaT.jl @@ -61,19 +61,19 @@ function HH(v,I_Na,I_K,I_L) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(m) - ((amf(v) * (1-m)) - (bmf(v) * m)) - der(h) - ((ahf(v) * (1-h)) - (bhf(v) * h)) - der(n) - ((anf(v) * (1-n)) - (bnf(v) * n)) - - g_Na - (m^3 * h * gbar_Na) - g_K - (n^4 * gbar_K) - - I_Na - (g_Na * (v - E_Na)) - I_K - (g_K * (v - E_K)) - I_L - (g_L * (v - E_L)) - - } + @equations begin + der(m) = (amf(v) * (1-m)) - (bmf(v) * m) + der(h) = (ahf(v) * (1-h)) - (bhf(v) * h) + der(n) = (anf(v) * (1-n)) - (bnf(v) * n) + + g_Na = m^3 * h * gbar_Na + g_K = n^4 * gbar_K + + I_Na = g_Na * (v - E_Na) + I_K = g_K * (v - E_K) + I_L = g_L * (v - E_L) + + end end @@ -93,24 +93,24 @@ function CaT(v,I_CaT) d = Unknown("d") s = Unknown("s") - { - der(r) - ((ralpha*(1-r)) - (rbeta*r)) - der(d) - ((dbeta*(1-s-d)) - (dalpha*d)) - der(s) - ((salpha*(1-s-d)) - (sbeta*s)) - - I_CaT - gbar_CaT*r*r*r*s*(v - E_Ca) + @equations begin + der(r) = ralpha*(1-r) - rbeta*r + der(d) = dbeta*(1-s-d) - dalpha*d + der(s) = salpha*(1-s-d) - sbeta*s + + I_CaT = gbar_CaT*r*r*r*s*(v - E_Ca) - ralpha - 1.0/(1.7+exp(-(v+28.2)/13.5)) - rbeta - exp(-(v+63.0)/7.8)/(exp(-(v+28.8)/13.1)+1.7) + ralpha = 1.0/(1.7+exp(-(v+28.2)/13.5)) + rbeta = exp(-(v+63.0)/7.8)/(exp(-(v+28.8)/13.1)+1.7) - salpha - exp(-(v+160.3)/17.8) - sbeta - (sqrt(0.25+exp((v+83.5)/6.3))-0.5) * (exp(-(v+160.3)/17.8)) + salpha = exp(-(v+160.3)/17.8) + sbeta = (sqrt(0.25+exp((v+83.5)/6.3))-0.5) * exp(-(v+160.3)/17.8) - bd - sqrt(0.25+exp((v+83.5)/6.3)) - dalpha - (1.0+exp((v+37.4)/30.0))/(240.0*(0.5+bd)) - dbeta - (bd-0.5)*dalpha + bd = sqrt(0.25+exp((v+83.5)/6.3)) + dalpha = (1.0+exp((v+37.4)/30.0))/(240.0*(0.5+bd)) + dbeta = (bd-0.5)*dalpha - } + end end @@ -121,12 +121,12 @@ function WRR() I_K = Unknown ("I_K") I_L = Unknown ("I_L") I_CaT = Unknown ("I_CaT") - { - HH(v,I_Na,I_K,I_L) - CaT(v,I_CaT) - - der(v) - ((I - (I_CaT + I_Na + I_K + I_L)) / C_m) - } + @equations begin + HH(v,I_Na,I_K,I_L) + CaT(v,I_CaT) + + der(v) = (I - (I_CaT + I_Na + I_K + I_L)) / C_m + end end wrr = WRR() # returns the hierarchical model diff --git a/examples/neural/adex.jl b/examples/neural/adex.jl index a2d09df..e5c6a2e 100644 --- a/examples/neural/adex.jl +++ b/examples/neural/adex.jl @@ -27,19 +27,19 @@ function AdEx() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(V) - (( ((- gL) * (V - EL)) + - (gL * Delta * (exp ((V - VT) / Delta))) + - (- W) + Isyn) / C) - der(W) - (((a * (V - EL)) - W) / tau_w) - - Event(V-theta, - { - reinit(V, Vr) - }, # positive crossing - {}) - - } + @equations begin + der(V) = (( ((- gL) * (V - EL)) + + (gL * Delta * (exp ((V - VT) / Delta))) + + (- W) + Isyn) / C) + der(W) = (((a * (V - EL)) - W) / tau_w) + + Event(V-theta, + Equation[ + reinit(V, Vr) + ], # positive crossing + Equation[]) + + end end diff --git a/examples/neural/fhn.jl b/examples/neural/fhn.jl index 41ca78f..b9a1747 100644 --- a/examples/neural/fhn.jl +++ b/examples/neural/fhn.jl @@ -24,10 +24,10 @@ tau = 12.5 function FitzHughNagumo(Iext) v = Unknown("v") w = Unknown("w") - { - der(v) - (v - (v^3 / 3) - w + Iext) # == 0 is assumed - der(w) - (v + a - b * w) / tau - } + @equations begin + der(v) = v - (v^3 / 3) - w + Iext + der(w) = (v + a - b * w) / tau + end end v = FitzHughNagumo(0.5) # returns the hierarchical model diff --git a/examples/neural/hh.jl b/examples/neural/hh.jl index 4b78a46..73ed02d 100644 --- a/examples/neural/hh.jl +++ b/examples/neural/hh.jl @@ -66,21 +66,21 @@ function HodgkinHuxley() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { + @equations begin - der(v) - ((I - (I_Na + I_K + I_L)) / C_m) - der(m) - ((amf(v) * (1 - m)) - (bmf(v) * m)) - der(h) - ((ahf(v) * (1 - h)) - (bhf(v) * h)) - der(n) - ((anf(v) * (1 - n)) - (bnf(v) * n)) + der(v) = (I - (I_Na + I_K + I_L)) / C_m + der(m) = (amf(v) * (1 - m)) - (bmf(v) * m) + der(h) = (ahf(v) * (1 - h)) - (bhf(v) * h) + der(n) = (anf(v) * (1 - n)) - (bnf(v) * n) - g_Na - (m^3 * h * gbar_Na) - g_K - (n^4 * gbar_K) + g_Na = m^3 * h * gbar_Na + g_K = n^4 * gbar_K - I_Na - (g_Na * (v - E_Na)) - I_K - (g_K * (v - E_K)) - I_L - (g_L * (v - E_L)) + I_Na = g_Na * (v - E_Na) + I_K = g_K * (v - E_K) + I_L = g_L * (v - E_L) - } + end end hh = HodgkinHuxley() # returns the hierarchical model diff --git a/examples/neural/hh_base.jl b/examples/neural/hh_base.jl index 6c533bd..1efa248 100644 --- a/examples/neural/hh_base.jl +++ b/examples/neural/hh_base.jl @@ -38,15 +38,15 @@ function Na_model(v,I_Na) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(m) - ((amf(v) * (1-m)) - (bmf(v) * m)) - der(h) - ((ahf(v) * (1-h)) - (bhf(v) * h)) + @equations begin + der(m) = (amf(v) * (1-m)) - (bmf(v) * m) + der(h) = (ahf(v) * (1-h)) - (bhf(v) * h) + + g_Na = m^3 * h * gbar_Na + + I_Na = g_Na * (v - E_Na) - g_Na - (m^3 * h * gbar_Na) - - I_Na - (g_Na * (v - E_Na)) - - } + end end ## K current @@ -64,21 +64,21 @@ function K_model(v,I_K) n = Unknown(0.317, "n") g_K = Unknown () - { - der(n) - ((anf(v) * (1 - n)) - (bnf(v) * n)) + @equations begin + der(n) = (anf(v) * (1 - n)) - (bnf(v) * n) - g_K - (n^4 * gbar_K) + g_K = n^4 * gbar_K - I_K - (g_K * (v - E_K)) - } + I_K = g_K * (v - E_K) + end end ## Leak current function Leak_model(v,I_L) - { - I_L - (g_L * (v - E_L)) - } + @equations begin + I_L = g_L * (v - E_L) + end end function HodgkinHuxley(I) @@ -92,15 +92,15 @@ function HodgkinHuxley(I) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { + @equations begin - Na_model(v,I_Na) - K_model(v,I_K) - Leak_model(v,I_L) - - der(v) - ((I - (I_Na + I_K + I_L)) / C_m) + Na_model(v,I_Na) + K_model(v,I_K) + Leak_model(v,I_L) + + der(v) = (I - (I_Na + I_K + I_L)) / C_m - } + end end function simrun (sim, tf, dt) diff --git a/examples/neural/hh_modular.jl b/examples/neural/hh_modular.jl index a0d69f2..14bbd84 100644 --- a/examples/neural/hh_modular.jl +++ b/examples/neural/hh_modular.jl @@ -43,15 +43,15 @@ function Na_model(v,I_Na) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(m) - ((amf(v) * (1-m)) - (bmf(v) * m)) - der(h) - ((ahf(v) * (1-h)) - (bhf(v) * h)) + @equations begin + der(m) = (amf(v) * (1-m)) - (bmf(v) * m) + der(h) = (ahf(v) * (1-h)) - (bhf(v) * h) + + g_Na = m^3 * h * gbar_Na + + I_Na = g_Na * (v - E_Na) - g_Na - (m^3 * h * gbar_Na) - - I_Na - (g_Na * (v - E_Na)) - - } + end end ## K current @@ -69,21 +69,21 @@ function K_model(v,I_K) n = Unknown(0.317, "n") g_K = Unknown () - { - der(n) - ((anf(v) * (1 - n)) - (bnf(v) * n)) + @equations begin + der(n) = (anf(v) * (1 - n)) - (bnf(v) * n) - g_K - (n^4 * gbar_K) + g_K = n^4 * gbar_K - I_K - (g_K * (v - E_K)) - } + I_K = g_K * (v - E_K) + end end ## Leak current function Leak_model(v,I_L) - { - I_L - (g_L * (v - E_L)) - } + @equations begin + I_L = g_L * (v - E_L) + end end function HodgkinHuxley() @@ -97,15 +97,15 @@ function HodgkinHuxley() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { + @equations begin - Na_model(v,I_Na) - K_model(v,I_K) - Leak_model(v,I_L) - - der(v) - ((I - (I_Na + I_K + I_L)) / C_m) + Na_model(v,I_Na) + K_model(v,I_K) + Leak_model(v,I_L) + + der(v) = (I - (I_Na + I_K + I_L)) / C_m - } + end end hh = HodgkinHuxley() # returns the hierarchical model diff --git a/examples/neural/hr.jl b/examples/neural/hr.jl index 6aa1949..2778b1c 100644 --- a/examples/neural/hr.jl +++ b/examples/neural/hr.jl @@ -36,11 +36,11 @@ function HindmarshRose(I) x = Unknown(-1.0,"x") y = Unknown("y") z = Unknown("z") - { - der(x) - (y + phi_(x) - z + I ) - der(y) - (psi(x) - y) - der(z) - (r * (s * (x - xr) - z)) - } + @equations begin + der(x) = y + phi_(x) - z + I + der(y) = psi(x) - y + der(z) = r * (s * (x - xr) - z) + end end v = HindmarshRose(0.1) # returns the hierarchical model diff --git a/examples/neural/iaf.jl b/examples/neural/iaf.jl index 5905b17..e8a86b1 100644 --- a/examples/neural/iaf.jl +++ b/examples/neural/iaf.jl @@ -20,16 +20,16 @@ function LeakyIaF() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(v) - (( ((- gL) * (v - vL)) + Isyn) / C) + @equations begin + der(v) = ( ((- gL) * (v - vL)) + Isyn) / C - Event(v-theta, - { - reinit(v, vreset) - }, # positive crossing - {}) + Event(v-theta, + Equation[ + reinit(v, vreset) + ], # positive crossing + Equation[]) - } + end end diff --git a/examples/neural/iafrefr.jl b/examples/neural/iafrefr.jl index 0701d4b..4dc4d33 100644 --- a/examples/neural/iafrefr.jl +++ b/examples/neural/iafrefr.jl @@ -17,9 +17,9 @@ function Subthreshold(v) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(v) - (( ((- gL) * (v - vL)) + Isyn) / C) - } + @equations begin + der(v) = ( ((- gL) * (v - vL)) + Isyn) / C + end end @@ -34,13 +34,13 @@ function Refractory(v,trefr) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - StructuralEvent(MTime - trefr, - # when the end of refractory period is reached, - # switch back to subthreshold mode - RefractoryEq(v), - () -> LeakyIaF()) - } + @equations begin + StructuralEvent(MTime - trefr, + # when the end of refractory period is reached, + # switch back to subthreshold mode + RefractoryEq(v), + () -> LeakyIaF()) + end end @@ -48,16 +48,16 @@ end function LeakyIaF() v = Unknown(vreset, "v") - { - StructuralEvent(v-theta, - # when v crosses the threshold, - # switch to refractory mode - Subthreshold(v), - () -> begin - trefr = value(MTime)+trefractory - Refractory(v,trefr) - end) - } + @equations begin + StructuralEvent(v-theta, + # when v crosses the threshold, + # switch to refractory mode + Subthreshold(v), + () -> begin + trefr = value(MTime)+trefractory + Refractory(v,trefr) + end) + end end diff --git a/examples/neural/iafsyn.jl b/examples/neural/iafsyn.jl index ae47212..c8391e0 100644 --- a/examples/neural/iafsyn.jl +++ b/examples/neural/iafsyn.jl @@ -52,16 +52,16 @@ function LeakyIaF(V,Isyn) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(V) - (( ((- gL) * (V - vL)) + Isyn) / C) - - Event(V-theta, - { - reinit(V, vreset) - }, # positive crossing - {}) - - } + @equations begin + der(V) = ( ((- gL) * (V - vL)) + Isyn) / C + + Event(V-theta, + Equation[ + reinit(V, vreset) + ], # positive crossing + Equation[]) + + end end @@ -74,25 +74,25 @@ function Syn(V,Isyn,input) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(S) - (alpha * (1 - S) - beta * S) - der(SS) - ((s0 - SS) / taus) - - Isyn - (gsyn * (V - vsyn)) - gsyn - (gsmax * S * SS) - - Event(grid_input(input), - { - reinit(SS, SS + f * (1 - SS)) - }, - {}) - Event(V-theta, - { - reinit(S, 0.0) - reinit(SS, 0.0) - }, - {}) - } + @equations begin + der(S) = (alpha * (1 - S) - beta * S) + der(SS) = ((s0 - SS) / taus) + + Isyn = (gsyn * (V - vsyn)) + gsyn = (gsmax * S * SS) + + Event(grid_input(input), + Equation[ + reinit(SS, SS + f * (1 - SS)) + ], + Equation[]) + Event(V-theta, + Equation[ + reinit(S, 0.0) + reinit(SS, 0.0) + ], + Equation[]) + end end @@ -101,11 +101,11 @@ function Circuit(y) V = Voltage (-35.0, "V") Isyn = Unknown ("Isyn") Isyn1 = Unknown () - { - LeakyIaF(V,Isyn) - Syn(V,Isyn1,y) - Isyn - Isyn1 - } + @equations begin + LeakyIaF(V,Isyn) + Syn(V,Isyn1,y) + Isyn = Isyn1 + end end diff --git a/examples/neural/izhfs.jl b/examples/neural/izhfs.jl index 21421ae..6758e31 100644 --- a/examples/neural/izhfs.jl +++ b/examples/neural/izhfs.jl @@ -30,18 +30,18 @@ function IzhikevichFS() # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(v) - (((k * (v - Vr) * (v - Vt)) + (- u) + Iext) / Cm) - der(u) - (FS_a * (s - u)) - s - (FS_b * (v - Vb) ^ 3) - - Event(v-Vpeak, - { - reinit(v, FS_c) - }, # positive crossing - {}) - - } + @equations begin + der(v) = ((k * (v - Vr) * (v - Vt)) + (- u) + Iext) / Cm + der(u) = FS_a * (s - u) + s = FS_b * (v - Vb) ^ 3 + + Event(v-Vpeak, + Equation[ + reinit(v, FS_c) + ], # positive crossing + Equation[]) + + end end diff --git a/examples/neural/ml.jl b/examples/neural/ml.jl index 927617a..c297f37 100644 --- a/examples/neural/ml.jl +++ b/examples/neural/ml.jl @@ -47,10 +47,10 @@ function MorrisLecar() # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). @equations begin - der(v) - ((Istim + (gl * (vl - v)) + ica + ik) / c) # == 0 is assumed - der(w) - (lamw (v) * (winf(v) - w)) - ica - (gca * (minf (v) * (vca - v))) - ik - (gk * (w * (vk - v))) + der(v) = (Istim + (gl * (vl - v)) + ica + ik) / c # == 0 is assumed + der(w) = lamw (v) * (winf(v) - w) + ica = gca * (minf (v) * (vca - v)) + ik = gk * (w * (vk - v)) end end diff --git a/examples/neural/mlsyn.jl b/examples/neural/mlsyn.jl index fe21769..40b354b 100644 --- a/examples/neural/mlsyn.jl +++ b/examples/neural/mlsyn.jl @@ -52,12 +52,12 @@ function MorrisLecar(v) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(v) - ((Istim + (gl * (vl - v)) + ica + ik) / c) # == 0 is assumed - der(w) - (lamw (v) * (winf(v) - w)) - ica - (gca * (minf (v) * (vca - v))) - ik - (gk * (w * (vk - v))) - } + @equations begin + der(v) = (Istim + (gl * (vl - v)) + ica + ik) / c + der(w) = lamw (v) * (winf(v) - w) + ica = gca * (minf (v) * (vca - v)) + ik = gk * (w * (vk - v)) + end end @@ -67,19 +67,19 @@ end function Syn(v,s) - { - der(s) = alpha * k_(v) * (1-s) - beta * s - } + @equations begin + der(s) = alpha * k_(v) * (1-s) - beta * s + end end function MLCircuit() s = Unknown(0.056, "s") v = Voltage(-24.22, "v") - { - Syn(v,s) - MorrisLecar(v) - } + @equations begin + Syn(v,s) + MorrisLecar(v) + end end diff --git a/examples/neural/netiafrefr.jl b/examples/neural/netiafrefr.jl index fa499c1..638ae20 100644 --- a/examples/neural/netiafrefr.jl +++ b/examples/neural/netiafrefr.jl @@ -17,30 +17,30 @@ function Subthreshold(v) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - der(v) - (( ((- gL) * (v - vL)) + Isyn) / C) - } + @equations begin + der(v) = ( ((- gL) * (v - vL)) + Isyn) / C + end end function RefractoryEq(v) - { - v - vreset - } + @equations begin + v = vreset + end end function Refractory(v,trefr) # The following gives the return value which is a list of equations. # Expressions with Unknowns are kept as expressions. Regular # variables are evaluated immediately (like normal). - { - StructuralEvent(MTime - trefr, - # when the end of refractory period is reached, - # switch back to subthreshold mode - RefractoryEq(v), - () -> LeakyIaF()) - } + @equations begin + StructuralEvent(MTime - trefr, + # when the end of refractory period is reached, + # switch back to subthreshold mode + RefractoryEq(v), + () -> LeakyIaF()) + end end @@ -48,16 +48,16 @@ end function LeakyIaF() v = Unknown(vreset, "v") - { - StructuralEvent(v-theta, - # when v crosses the threshold, - # switch to refractory mode - Subthreshold(v), - () -> begin - trefr = value(MTime)+trefractory - Refractory(v,trefr) - end) - } + @equations begin + StructuralEvent(v-theta, + # when v crosses the threshold, + # switch to refractory mode + Subthreshold(v), + () -> begin + trefr = value(MTime)+trefractory + Refractory(v,trefr) + end) + end end diff --git a/examples/neural/pr.jl b/examples/neural/pr.jl index 7c119a3..da9c7cc 100644 --- a/examples/neural/pr.jl +++ b/examples/neural/pr.jl @@ -86,11 +86,11 @@ function Soma(V,I) h = Unknown ("h") n = Unknown ("n") - { - der(V) - ((-gLs * (V - VL) - gNa * (Minfs(V)^2) * h * (V - VNa) - gKdr * n * (V - VK) + I + J/As) / Cm) - der(h) - (alphahs(V) - (alphahs(V) + betahs(V)) * h) - der(n) - (alphans(V) - (alphans(V) + betans(V)) * n) - } + @equations begin + der(V) = (-gLs * (V - VL) - gNa * (Minfs(V)^2) * h * (V - VNa) - gKdr * n * (V - VK) + I + J/As) / Cm + der(h) = alphahs(V) - (alphahs(V) + betahs(V)) * h + der(n) = alphans(V) - (alphans(V) + betans(V)) * n + end end @@ -105,26 +105,26 @@ function Dendrite(V,I) chid = Unknown ("chid") alphaqd = Unknown ("alphaqd") - { - der(V) - ((-gLd * (V - VL) - ICad - IKahp - IK + I) / Cm) - der(s) - (alphasd(V) - (alphasd(V) + betasd(V))*s) - der(c) - (alphacd(V) - (alphacd(V) + betacd(V))*c) - der(q) - (alphaqd - (alphaqd + betaqd) * q) - der(Cad) - (-0.13 * ICad - 0.075 * Cad) - ICad - (gCa * s * s * (V - VCa)) - IKahp - (gKahp * q * (V - VK)) - IK - (gKC * c * chid * (V - VK)) - alphaqd - (min(0.00002 * Cad, 0.01)) - chid - (min(Cad / 250.0,1.0)) - } + @equations begin + der(V) = (-gLd * (V - VL) - ICad - IKahp - IK + I) / Cm + der(s) = alphasd(V) - (alphasd(V) + betasd(V))*s + der(c) = alphacd(V) - (alphacd(V) + betacd(V))*c + der(q) = alphaqd - (alphaqd + betaqd) * q + der(Cad) = -0.13 * ICad - 0.075 * Cad + ICad = gCa * s * s * (V - VCa) + IKahp = gKahp * q * (V - VK) + IK = gKC * c * chid * (V - VK) + alphaqd = min(0.00002 * Cad, 0.01) + chid = min(Cad / 250.0,1.0) + end end function Connect(I,g,n1,n2) - { - g * I - (n2 - n1) - } + @equations begin + g * I = n2 - n1 + end end @@ -134,12 +134,12 @@ function PRCircuit() Vs = Voltage(-60.0,"Vs") Vd = Voltage(-60.0,"Vd") V - { - Soma(Vs,Is) - Dendrite(Vd,Id) - Connect(Is,(As/gc),Vs,Vd) - Connect(Id,(Ad/gc),Vd,Vs) - } + @equations begin + Soma(Vs,Is) + Dendrite(Vd,Id) + Connect(Is,(As/gc),Vs,Vd) + Connect(Id,(Ad/gc),Vd,Vs) + end end pr = PRCircuit() # returns the hierarchical model diff --git a/examples/neural/theta.jl b/examples/neural/theta.jl index 5474d6f..0859a50 100644 --- a/examples/neural/theta.jl +++ b/examples/neural/theta.jl @@ -11,24 +11,22 @@ beta = 0.01 function Theta(theta) - { - der(theta) - (1 - cos(theta) + alpha * (1 + cos(theta)) * sin(beta * MTime)) - - Event(theta-pi, - { - reinit(theta, theta-2*pi) - }, # positive crossing - {}) - - } + @equations begin + der(theta) = 1 - cos(theta) + alpha * (1 + cos(theta)) * sin(beta * MTime) + + Event(theta-pi, + Equation[ + reinit(theta, theta-2*pi) + ], # positive crossing + Equation[]) + + end end function ThetaCircuit() theta = Unknown(1.0, "theta") - { - Theta(theta) - } + Theta(theta) end diff --git a/examples/run_all.jl b/examples/run_all.jl index 6db61b9..9096c5e 100644 --- a/examples/run_all.jl +++ b/examples/run_all.jl @@ -15,7 +15,7 @@ include(path * "basics/breaking_pendulum_in_box.jl") include(path * "basics/breaking_pendulum.jl") -## include(path * "circuit_complex.jl") +## include(path * "basics/circuit_complex.jl") include(path * "basics/dc_motor_w_shaft.jl") diff --git a/examples/stdlib/m_blocks.jl b/examples/stdlib/m_blocks.jl index fc17a17..e040b30 100644 --- a/examples/stdlib/m_blocks.jl +++ b/examples/stdlib/m_blocks.jl @@ -20,30 +20,30 @@ function ex_PID_Controller() s3 = Unknown("s3") s4 = Unknown("s4") s5 = Unknown("s5") - { - ## KinematicPTP(s1, 0.5, driveAngle, 1.0, 1.0) - ## Integrator(s1, s2, 1.0) - s1 - ifelse((MTime < 0.5) | (MTime >= 3.2), 0.0, - ifelse(MTime < 1.5, MTime - 0.5, - ifelse(MTime < 2.2, 1.0, 3.2 - MTime))) - SpeedSensor(n2, s3) - SpeedSensor(n3, s4) - LimPID(s2, s3, s4, - controllerType = "PI", - k = 100.0, - Ti = 0.1, - Td = 0.1, - yMax = 12.0, - Ni = 0.1) - SignalTorque(n1, 0.0, s4) - s5 - (s4 - s2) - InitialEquation(der(s5) - 0.0) # force a constant initial spin of the shaft to tame initial conditions - SignalTorque(n1, 0.0, s3) - Inertia(n1, n2, 1.0) - SpringDamper(n2, n3, 1e4, 100) - Inertia(n3, n4, 2.0) - SignalTorque(n4, 0.0, 10.0) - } + @equations begin + ## KinematicPTP(s1, 0.5, driveAngle, 1.0, 1.0) + ## Integrator(s1, s2, 1.0) + s1 = ifelse((MTime < 0.5) | (MTime >= 3.2), 0.0, + ifelse(MTime < 1.5, MTime - 0.5, + ifelse(MTime < 2.2, 1.0, 3.2 - MTime))) + SpeedSensor(n2, s3) + SpeedSensor(n3, s4) + LimPID(s2, s3, s4, + controllerType = "PI", + k = 100.0, + Ti = 0.1, + Td = 0.1, + yMax = 12.0, + Ni = 0.1) + SignalTorque(n1, 0.0, s4) + s5 = s4 - s2 + InitialEquation(der(s5) - 0.0) # force a constant initial spin of the shaft to tame initial conditions + SignalTorque(n1, 0.0, s3) + Inertia(n1, n2, 1.0) + SpringDamper(n2, n3, 1e4, 100) + Inertia(n3, n4, 2.0) + SignalTorque(n4, 0.0, 10.0) + end end # Results of this example match Dymola with the exception of # starting transients. This example only solves if info[11] = 0 diff --git a/examples/stdlib/m_electrical.jl b/examples/stdlib/m_electrical.jl index 6c6baa8..7f80031 100644 --- a/examples/stdlib/m_electrical.jl +++ b/examples/stdlib/m_electrical.jl @@ -13,18 +13,18 @@ function ex_CauerLowPassAnalog() n3 = Voltage("n3") n4 = Voltage("n4") g = 0.0 - { - StepVoltage(n1, g, 1.0, 1.0, 0.0) - Resistor(n1, n2, 1.0) - Capacitor(n2, g, 1.072) - Capacitor(n2, n3, 1/(1.704992^2 * 1.304)) - Inductor(n2, n3, 1.304) - Capacitor(n3, g, 1.682) - Capacitor(n3, n4, 1/(1.179945^2 * 0.8586)) - Inductor(n3, n4, 0.8565) - Capacitor(n4, g, 0.7262) - Resistor(n4, g, 1.0) - } + Equation[ + StepVoltage(n1, g, 1.0, 1.0, 0.0) + Resistor(n1, n2, 1.0) + Capacitor(n2, g, 1.072) + Capacitor(n2, n3, 1/(1.704992^2 * 1.304)) + Inductor(n2, n3, 1.304) + Capacitor(n3, g, 1.682) + Capacitor(n3, n4, 1/(1.179945^2 * 0.8586)) + Inductor(n3, n4, 0.8565) + Capacitor(n4, g, 0.7262) + Resistor(n4, g, 1.0) + ] end function sim_CauerLowPassAnalog() @@ -49,34 +49,34 @@ function ex_CauerLowPassOPV() c3 = 1.682 c4 = 1 / (1.179945^2 + l2) c5 = 0.7262 - { - StepVoltage(n[1], g, -1.0, 1.0, 0.0) - IdealOpAmp(g, n[2], n[3]) - IdealOpAmp(g, n[4], n[5]) - IdealOpAmp(g, n[6], n[7]) - IdealOpAmp(g, n[8], n[9]) - IdealOpAmp(g, n[10], n[11]) - Resistor(n[1], n[2], 1.0) - Resistor(n[3], n[4], -1.0) - Resistor(n[5], n[6], 1.0) - Resistor(n[7], n[8], -1.0) - Resistor(n[9], n[10], 1.0) - Capacitor(n[2], n[3], c1 + c2) - Capacitor(n[4], n[5], l1) - Capacitor(n[6], n[7], c2 + c3 + c4) - Capacitor(n[8], n[9], l2) - Capacitor(n[10], n[11], c4 + c5) - Resistor(n[2], n[3], 1.0) - Resistor(n[2], n[5], 1.0) - Resistor(n[4], n[7], -1.0) - Resistor(n[6], n[9], 1.0) - Resistor(n[8], n[11], -1.0) - Resistor(n[10], n[11], 1.0) - Capacitor(n[2], n[7], c2) - Capacitor(n[3], n[6], c2) - Capacitor(n[6], n[11], c4) - Capacitor(n[7], n[10], c4) - } + Equation[ + StepVoltage(n[1], g, -1.0, 1.0, 0.0) + IdealOpAmp(g, n[2], n[3]) + IdealOpAmp(g, n[4], n[5]) + IdealOpAmp(g, n[6], n[7]) + IdealOpAmp(g, n[8], n[9]) + IdealOpAmp(g, n[10], n[11]) + Resistor(n[1], n[2], 1.0) + Resistor(n[3], n[4], -1.0) + Resistor(n[5], n[6], 1.0) + Resistor(n[7], n[8], -1.0) + Resistor(n[9], n[10], 1.0) + Capacitor(n[2], n[3], c1 + c2) + Capacitor(n[4], n[5], l1) + Capacitor(n[6], n[7], c2 + c3 + c4) + Capacitor(n[8], n[9], l2) + Capacitor(n[10], n[11], c4 + c5) + Resistor(n[2], n[3], 1.0) + Resistor(n[2], n[5], 1.0) + Resistor(n[4], n[7], -1.0) + Resistor(n[6], n[9], 1.0) + Resistor(n[8], n[11], -1.0) + Resistor(n[10], n[11], 1.0) + Capacitor(n[2], n[7], c2) + Capacitor(n[3], n[6], c2) + Capacitor(n[6], n[11], c4) + Capacitor(n[7], n[10], c4) + ] end function ex_CauerLowPassOPV2() n1 = Voltage("n1") @@ -98,34 +98,34 @@ function ex_CauerLowPassOPV2() c3 = 1.682 c4 = 1 / (1.179945^2 + l2) c5 = 0.7262 - { - StepVoltage(n1, g, -1.0, 1.0, 0.0) - IdealOpAmp(g, n2, n3) - IdealOpAmp(g, n4, n5) - IdealOpAmp(g, n6, n7) - IdealOpAmp(g, n8, n9) - IdealOpAmp(g, n10, n11) - Resistor(n1, n2, 1.0) - Resistor(n3, n4, -1.0) - Resistor(n5, n6, 1.0) - Resistor(n7, n8, -1.0) - Resistor(n9, n10, 1.0) - Capacitor(n2, n3, c1 + c2) - Capacitor(n4, n5, l1) - Capacitor(n6, n7, c2 + c3 + c4) - Capacitor(n8, n9, l2) - Capacitor(n10, n11, c4 + c5) - Resistor(n2, n3, 1.0) - Resistor(n2, n5, 1.0) - Resistor(n4, n7, -1.0) - Resistor(n6, n9, 1.0) - Resistor(n8, n11, -1.0) - Resistor(n10, n11, 1.0) - Capacitor(n2, n7, c2) - Capacitor(n3, n6, c2) - Capacitor(n6, n11, c4) - Capacitor(n7, n10, c4) - } + Equation[ + StepVoltage(n1, g, -1.0, 1.0, 0.0) + IdealOpAmp(g, n2, n3) + IdealOpAmp(g, n4, n5) + IdealOpAmp(g, n6, n7) + IdealOpAmp(g, n8, n9) + IdealOpAmp(g, n10, n11) + Resistor(n1, n2, 1.0) + Resistor(n3, n4, -1.0) + Resistor(n5, n6, 1.0) + Resistor(n7, n8, -1.0) + Resistor(n9, n10, 1.0) + Capacitor(n2, n3, c1 + c2) + Capacitor(n4, n5, l1) + Capacitor(n6, n7, c2 + c3 + c4) + Capacitor(n8, n9, l2) + Capacitor(n10, n11, c4 + c5) + Resistor(n2, n3, 1.0) + Resistor(n2, n5, 1.0) + Resistor(n4, n7, -1.0) + Resistor(n6, n9, 1.0) + Resistor(n8, n11, -1.0) + Resistor(n10, n11, 1.0) + Capacitor(n2, n7, c2) + Capacitor(n3, n6, c2) + Capacitor(n6, n11, c4) + Capacitor(n7, n10, c4) + ] end function sim_CauerLowPassOPV() @@ -159,17 +159,17 @@ function ex_CharacteristicIdealDiodes() n2 = Voltage("n2") n3 = Voltage("n3") g = 0.0 - { - SineVoltage(s1, g, 10.0, 1.0, -pi/10.0) - SineVoltage(s2, g, 10.0, 1.0, -pi/15.0, -9.0) - SineVoltage(s3, g, 10.0, 1.0, -pi/20.0) - Resistor(n1, g, 1e-3) - Resistor(n2, g, 1e-3) - Resistor(n3, g, 1e-3) - IdealDiode(s1, n1, 0.0, 0.0, 0.0) - IdealDiode(s2, n2, 0.0, 0.1, 0.1) - IdealDiode(s3, n3, 5.0, 0.2, 0.2) - } + Equation[ + SineVoltage(s1, g, 10.0, 1.0, -pi/10.0) + SineVoltage(s2, g, 10.0, 1.0, -pi/15.0, -9.0) + SineVoltage(s3, g, 10.0, 1.0, -pi/20.0) + Resistor(n1, g, 1e-3) + Resistor(n2, g, 1e-3) + Resistor(n3, g, 1e-3) + IdealDiode(s1, n1, 0.0, 0.0, 0.0) + IdealDiode(s2, n2, 0.0, 0.1, 0.1) + IdealDiode(s3, n3, 5.0, 0.2, 0.2) + ] end function sim_CharacteristicIdealDiodes() @@ -182,10 +182,10 @@ end # n1 = Voltage("n1") # g = 0.0 # { -# SineVoltage(s1, g, 10.0, 1.0, -pi/10.0, 0.0) -# Resistor(n1, g, 1e-3) -# IdealDiode(s1, n1, 0.0, 0.0, 0.0) -# } +# SineVoltage(s1, g, 10.0, 1.0, -pi/10.0, 0.0) +# Resistor(n1, g, 1e-3) +# IdealDiode(s1, n1, 0.0, 0.0, 0.0) +# } # end # m = ex_CharacteristicIdealDiodes() @@ -202,20 +202,20 @@ function ex_ChuaCircuit() function NonlinearResistor(n1::ElectricalNode, n2::ElectricalNode, Ga, Gb, Ve) i = Current(compatible_values(n1, n2)) v = Voltage(compatible_values(n1, n2)) - { - Branch(n1, n2, v, i) - i - ifelse(v < -Ve, Gb .* (v + Ve) - Ga .* Ve, - ifelse(v > Ve, Gb .* (v - Ve) + Ga*Ve, Ga*v)) - } + Equation[ + Branch(n1, n2, v, i) + i - ifelse(v < -Ve, Gb .* (v + Ve) - Ga .* Ve, + ifelse(v > Ve, Gb .* (v - Ve) + Ga*Ve, Ga*v)) + ] end - { - Resistor(n1, g, 12.5e-3) - Inductor(n1, n2, 18.0) - Resistor(n2, n3, 1 / 0.565) - Capacitor(n2, g, 100.0) - Capacitor(n3, g, 10.0) - NonlinearResistor(n3, g, -0.757576, -0.409091, 1.0) - } + Equation[ + Resistor(n1, g, 12.5e-3) + Inductor(n1, n2, 18.0) + Resistor(n2, n3, 1 / 0.565) + Capacitor(n2, g, 100.0) + Capacitor(n3, g, 10.0) + NonlinearResistor(n3, g, -0.757576, -0.409091, 1.0) + ] end function sim_ChuaCircuit() @@ -234,12 +234,12 @@ function ex_HeatingResistor() hp1 = Temperature("hp1") hp2 = Temperature("hp2") g = 0.0 - { - SineVoltage(n1, g, 220, 1.0) - Resistor(n1, g, 100.0, hp1, 20 + 273.15, 1e-3) - ThermalConductor(hp1, hp2, 50.0) - FixedTemperature(hp2, 20 + 273.15) - } + Equation[ + SineVoltage(n1, g, 220, 1.0) + Resistor(n1, g, 100.0, hp1, 20 + 273.15, 1e-3) + ThermalConductor(hp1, hp2, 50.0) + FixedTemperature(hp2, 20 + 273.15) + ] end function sim_HeatingResistor() @@ -252,14 +252,14 @@ function ex_HeatingRectifier() hp1 = Temperature("hp1") hp2 = Temperature("hp2") g = 0.0 - { - SineVoltage(n1, g, 1.0, 1.0) - HeatingDiode(n1, n2, hp1, morejunkXX) - Capacitor(n2, g, 1.0) - Resistor(n2, g, 1.0) - ThermalConductor(hp1, hp2, 10.0) - HeatCapacitor(hp2, 20 + 273.15) - } + Equation[ + SineVoltage(n1, g, 1.0, 1.0) + HeatingDiode(n1, n2, hp1, morejunkXX) + Capacitor(n2, g, 1.0) + Resistor(n2, g, 1.0) + ThermalConductor(hp1, hp2, 10.0) + HeatCapacitor(hp2, 20 + 273.15) + ] end function sim_HeatingRectifier(BROKEN) @@ -285,18 +285,18 @@ function ex_Rectifier() Vknee = 2.0 CDC = 15e-3 IDC = 500.0 - { - SineVoltage(n1, g, VAC .* sqrt(2/3), f, [0.0, -2pi/3, 2pi/3]) - Inductor(n1, n2, LAC) - ## Resistor(n2, np, 1e3) - ## Resistor(n2, nn, 1e3) - IdealDiode(n2, np, Vknee, Ron, Goff) - IdealDiode(nn, n2, Vknee, Ron, Goff) - ## Capacitor(np, nn, CDC) - Capacitor(np, g, 2 * CDC) - Capacitor(nn, g, 2 * CDC) - ## SignalCurrent(np, nn, IDC) - } + Equation[ + SineVoltage(n1, g, VAC .* sqrt(2/3), f, [0.0, -2pi/3, 2pi/3]) + Inductor(n1, n2, LAC) + ## Resistor(n2, np, 1e3) + ## Resistor(n2, nn, 1e3) + IdealDiode(n2, np, Vknee, Ron, Goff) + IdealDiode(nn, n2, Vknee, Ron, Goff) + ## Capacitor(np, nn, CDC) + Capacitor(np, g, 2 * CDC) + Capacitor(nn, g, 2 * CDC) + ## SignalCurrent(np, nn, IDC) + ] end function ex_Rectifier() @@ -317,23 +317,23 @@ function ex_Rectifier() Vknee = 2.0 CDC = 15e-3 IDC = 500.0 - { - SineVoltage(n1, g, VAC .* sqrt(2/3), f, [0.0, -2pi/3, 2pi/3]) - Inductor(n1, n2, LAC) - ## Resistor(n2, np, 1e3) - ## Resistor(n2, nn, 1e3) - IdealDiode(n2[1], np, Vknee, Ron, Goff) - IdealDiode(n2[2], np, Vknee, Ron, Goff) - IdealDiode(n2[3], np, Vknee, Ron, Goff) - IdealDiode(nn, n2[1], Vknee, Ron, Goff) - IdealDiode(nn, n2[2], Vknee, Ron, Goff) - IdealDiode(nn, n2[3], Vknee, Ron, Goff) - ## Capacitor(np, nn, CDC) - Capacitor(np, g, 2 * CDC) - Capacitor(nn, g, 2 * CDC) - ## SignalCurrent(np, nn, IDC) - Resistor(np, nn, 400 / IDC) - } + Equation[ + SineVoltage(n1, g, VAC .* sqrt(2/3), f, [0.0, -2pi/3, 2pi/3]) + Inductor(n1, n2, LAC) + ## Resistor(n2, np, 1e3) + ## Resistor(n2, nn, 1e3) + IdealDiode(n2[1], np, Vknee, Ron, Goff) + IdealDiode(n2[2], np, Vknee, Ron, Goff) + IdealDiode(n2[3], np, Vknee, Ron, Goff) + IdealDiode(nn, n2[1], Vknee, Ron, Goff) + IdealDiode(nn, n2[2], Vknee, Ron, Goff) + IdealDiode(nn, n2[3], Vknee, Ron, Goff) + ## Capacitor(np, nn, CDC) + Capacitor(np, g, 2 * CDC) + Capacitor(nn, g, 2 * CDC) + ## SignalCurrent(np, nn, IDC) + Resistor(np, nn, 400 / IDC) + ] end function sim_Rectifier(BROKEN) @@ -360,11 +360,11 @@ function ex_ShowSaturatingInductor() f = 1/(2pi) phase = pi/2 phase = 0.0 - { - SineVoltage(n1, g, U, f, phase) - ## SaturatingInductor(n1, g, Inom, Lnom, Linf, Lzer) - SaturatingInductor(n1, g, Inom, Lnom) - } + Equation[ + SineVoltage(n1, g, U, f, phase) + ## SaturatingInductor(n1, g, Inom, Lnom, Linf, Lzer) + SaturatingInductor(n1, g, Inom, Lnom) + ] end function sim_ShowSaturatingInductor(BROKEN) @@ -378,10 +378,10 @@ function ex_ShowSaturatingInductor2() U = 1.25 f = 1/(2pi) phase = 0.0 - { - SineVoltage(n1, g, U, f, phase) - SaturatingInductor2(n1, g, 71.0, 0.1, 0.04) - } + Equation[ + SineVoltage(n1, g, U, f, phase) + SaturatingInductor2(n1, g, 71.0, 0.1, 0.04) + ] end ## m = ex_ShowSaturatingInductor2() ## f = elaborate(m) @@ -395,10 +395,10 @@ function ex_ShowSaturatingInductor3() U = 1.25 f = 1/(2pi) phase = 0.0 - { - SineVoltage(n1, g, U, f, phase) - SaturatingInductor3(n1, g, 3e-9, 0.33, 0.15) - } + Equation[ + SineVoltage(n1, g, U, f, phase) + SaturatingInductor3(n1, g, 3e-9, 0.33, 0.15) + ] end ## m = ex_ShowSaturatingInductor3() ## f = elaborate(m) @@ -412,10 +412,10 @@ function ex_ShowSaturatingInductor4() U = 1.25 f = 1/(2pi) phase = 0.0 - { - SineVoltage(n1, g, U, f, phase) - SaturatingInductor4(n1, g, .50, 0.7, 0.0) - } + Equation[ + SineVoltage(n1, g, U, f, phase) + SaturatingInductor4(n1, g, .50, 0.7, 0.0) + ] end ## m = ex_ShowSaturatingInductor4() ## f = elaborate(m) @@ -434,18 +434,18 @@ function ex_ShowVariableResistor() isig2 = Voltage("Ir2") vres = Voltage("Vres") g = 0.0 - { - n2 - n3 - isig1 # current monitor - n5 - n3 - isig2 # current monitor - n5 - n4 - vres - SineVoltage(n, g, 1.0, 1.0) - Resistor(n, n1, 1.0) - Resistor(n1, n2, 1.0) - Resistor(n2, n3, 1.0) - Resistor(n, n4, 1.0) - Resistor(n4, n5, 2 + 2.5 * MTime) - Resistor(n5, n3, 1.0) - } + Equation[ + n2 - n3 - isig1 # current monitor + n5 - n3 - isig2 # current monitor + n5 - n4 - vres + SineVoltage(n, g, 1.0, 1.0) + Resistor(n, n1, 1.0) + Resistor(n1, n2, 1.0) + Resistor(n2, n3, 1.0) + Resistor(n, n4, 1.0) + Resistor(n4, n5, 2 + 2.5 * MTime) + Resistor(n5, n3, 1.0) + ] end function sim_ShowVariableResistor() @@ -463,18 +463,18 @@ function ex_ControlledSwitchWithArc() b3 = Voltage("b3") vs = Voltage("vs") g = 0.0 - { - SineVoltage(vs, g, 1.0, 1.0) - SignalVoltage(a1, g, 50.0) - ControlledIdealClosingSwitch(a1, a2, vs, 0.5, 1e-5, 1e-5) - Inductor(a2, a3, 0.1) - Resistor(a3, g, 1.0) - ## SignalVoltage(b1, g, 50.0) - ## ## ControlledCloserWithArc(b1, b2, vs, 0.5, 1e-5, 1e-5, 30.0, 1e4, 60.0) - ## ControlledCloserWithArc(b1, b2, vs, 0.5, 1e-5, 1e-5, 60.0, 1e-2, 60.0) - ## Inductor(b2, b3, 0.1) - ## Resistor(b3, g, 1.0) - } + Equation[ + SineVoltage(vs, g, 1.0, 1.0) + SignalVoltage(a1, g, 50.0) + ControlledIdealClosingSwitch(a1, a2, vs, 0.5, 1e-5, 1e-5) + Inductor(a2, a3, 0.1) + Resistor(a3, g, 1.0) + ## SignalVoltage(b1, g, 50.0) + ## ## ControlledCloserWithArc(b1, b2, vs, 0.5, 1e-5, 1e-5, 30.0, 1e4, 60.0) + ## ControlledCloserWithArc(b1, b2, vs, 0.5, 1e-5, 1e-5, 60.0, 1e-2, 60.0) + ## Inductor(b2, b3, 0.1) + ## Resistor(b3, g, 1.0) + ] end function sim_ControlledSwitchWithArc() @@ -491,11 +491,11 @@ function ex_dc() n1 = Voltage("n1") n2 = Voltage("n2") g = 0.0 - { - SignalVoltage(n1, g, 50.0) - Resistor(n1, n2, 0.1) - Inductor(n2, g, 0.1) - } + Equation[ + SignalVoltage(n1, g, 50.0) + Resistor(n1, n2, 0.1) + Inductor(n2, g, 0.1) + ] end ## m = ex_dc() @@ -509,14 +509,14 @@ function ex_CharacteristicThyristors() n3 = Voltage("n3") sig = Discrete(false) g = 0.0 - { - SineVoltage(n1, g, 10.0, 1.0, -0.006) - IdealThyristor(n1, n2, sig, 5.0) - IdealGTOThyristor(n1, n3, sig, 0.0) - BoolEvent(sig, MTime - 1.25) - Resistor(n2, g, 1e-3) - Resistor(n3, g, 1e-3) - } + Equation[ + SineVoltage(n1, g, 10.0, 1.0, -0.006) + IdealThyristor(n1, n2, sig, 5.0) + IdealGTOThyristor(n1, n3, sig, 0.0) + BoolEvent(sig, MTime - 1.25) + Resistor(n2, g, 1e-3) + Resistor(n3, g, 1e-3) + ] end function sim_CharacteristicThyristors() @@ -532,11 +532,11 @@ function test_BoolEventHook() addhook!(sig, reinit(sig2, false)) g = 0.0 - { - SineVoltage(n1, g, ifelse(sig2, 10.0, 5.0), ifelse(sig, 1.0, 2.0)) - BoolEvent(sig, MTime - 0.25) - Resistor(n1, g, 1e-3) - } + Equation[ + SineVoltage(n1, g, ifelse(sig2, 10.0, 5.0), ifelse(sig, 1.0, 2.0)) + BoolEvent(sig, MTime - 0.25) + Resistor(n1, g, 1e-3) + ] end ## m = test_BoolEventHook() diff --git a/examples/stdlib/m_heat_transfer.jl b/examples/stdlib/m_heat_transfer.jl index 832b483..2713707 100644 --- a/examples/stdlib/m_heat_transfer.jl +++ b/examples/stdlib/m_heat_transfer.jl @@ -12,11 +12,11 @@ using Sims function ex_TwoMasses() p1 = Temperature(373.15, "p1") p2 = Temperature(273.15, "p2") - { - HeatCapacitor(p1, 15.0, 373.15) - HeatCapacitor(p2, 15.0, 273.15) - ThermalConductor(p1, p2, 10.0) - } + Equation[ + HeatCapacitor(p1, 15.0, 373.15) + HeatCapacitor(p2, 15.0, 273.15) + ThermalConductor(p1, p2, 10.0) + ] end @@ -33,19 +33,19 @@ function ex_Motor(BROKEN) ## needs to have `interp` defined TAmb = 293.15 t = [0, 360, 360, 600] winding_losses = [100, 100, 1000, 1000] - { - # Winding - HeatCapacitor(p1, 2500.0, TAmb) - PrescribedHeatFlow(p1, interp(winding_losses, t, MTime), 95 + 273.15, 3.03E-3) - # Core - HeatCapacitor(p2, 25000.0, TAmb) - PrescribedHeatFlow(p2, 500.0) - # conduction between the winding and core: - ThermalConductor(p1, p2, 10.0) - # Convection to ambient - Convection(p2, p3, 25.0) - FixedTemperature(p3, TAmb) - } + Equation[ + # Winding + HeatCapacitor(p1, 2500.0, TAmb) + PrescribedHeatFlow(p1, interp(winding_losses, t, MTime), 95 + 273.15, 3.03E-3) + # Core + HeatCapacitor(p2, 25000.0, TAmb) + PrescribedHeatFlow(p2, 500.0) + # conduction between the winding and core: + ThermalConductor(p1, p2, 10.0) + # Convection to ambient + Convection(p2, p3, 25.0) + FixedTemperature(p3, TAmb) + ] end function sim_Motor() diff --git a/examples/stdlib/m_powersystems.jl b/examples/stdlib/m_powersystems.jl index 991f7a3..0d60f93 100644 --- a/examples/stdlib/m_powersystems.jl +++ b/examples/stdlib/m_powersystems.jl @@ -23,13 +23,13 @@ function ex_RLModel() ConductorLocation(-1.0, 10.0, Conductors["AAC 500 kcmil"]), ConductorLocation( 0.0, 10.0, Conductors["AAC 500 kcmil"]), ConductorLocation( 1.0, 10.0, Conductors["AAC 500 kcmil"]))) - { - SineVoltage(ns, g, Vln, freq, [0, -2/3*pi, 2/3*pi]) - SeriesProbe(ns, np, "I") - RLLine(np, nl, Z, len, freq) - ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, freq) - ## ConstZParallelLoad(nl, g, load_VA, load_pf, Vln, freq) - } + Equation[ + SineVoltage(ns, g, Vln, freq, [0, -2/3*pi, 2/3*pi]) + SeriesProbe(ns, np, "I") + RLLine(np, nl, Z, len, freq) + ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, freq) + ## ConstZParallelLoad(nl, g, load_VA, load_pf, Vln, freq) + ] end function sim_RLModel() @@ -64,15 +64,15 @@ function ex_PiModel() ConductorLocation(-1.0, 10.0, Conductors["AAC 500 kcmil"]), ConductorLocation( 0.0, 10.0, Conductors["AAC 500 kcmil"]), ConductorLocation( 1.0, 10.0, Conductors["AAC 500 kcmil"]))) - { - SineVoltage(ns, g, Vln, 60.0, [0, -2/3*pi, 2/3*pi]) - SeriesProbe(ns, np, "I") - PiLine(np, nl, Z, Y, len, freq, ne) - ## ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, freq) - ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, 60.0) - ## ConstZParallelLoad(nl, g, load_VA, load_pf, Vln, freq) - } -end + Equation[ + SineVoltage(ns, g, Vln, 60.0, [0, -2/3*pi, 2/3*pi]) + SeriesProbe(ns, np, "I") + PiLine(np, nl, Z, Y, len, freq, ne) + ## ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, freq) + ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, 60.0) + ## ConstZParallelLoad(nl, g, load_VA, load_pf, Vln, freq) + ] +end ## m = ex_PiModel() @@ -102,12 +102,12 @@ function ex_Modal() ConductorLocation(-1.0, 10.0, Conductors["AAC 500 kcmil"]), ConductorLocation( 0.0, 10.0, Conductors["AAC 500 kcmil"]), ConductorLocation( 1.0, 10.0, Conductors["AAC 500 kcmil"]))) - { - SineVoltage(ns, g, Vln, 60.0, [0, -2/3*pi, 2/3*pi]) - SeriesProbe(ns, np, "I") - ModalLine(np, nl, Z, Y, len, freq) - ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, 60.0) - } + Equation[ + SineVoltage(ns, g, Vln, 60.0, [0, -2/3*pi, 2/3*pi]) + SeriesProbe(ns, np, "I") + ModalLine(np, nl, Z, Y, len, freq) + ConstZSeriesLoad(nl, g, load_VA, load_pf, Vln, 60.0) + ] end diff --git a/examples/stdlib/m_rotational.jl b/examples/stdlib/m_rotational.jl index 77e24ba..da403a5 100644 --- a/examples/stdlib/m_rotational.jl +++ b/examples/stdlib/m_rotational.jl @@ -22,15 +22,15 @@ function ex_First() Jload = 2.0 ratio = 10.0 damping = 10.0 - { - SignalTorque(n1, g, amplitude * sin(2pi * freqHz * MTime)) - Inertia(n1, n2, Jmotor) - IdealGear(n2, n3, ratio) - Inertia(n3, n4, 2.0) - Damper(n4, g, damping) - Spring(n4, n5, 1e4) - Inertia(n5, n6, Jload) - } + Equation[ + SignalTorque(n1, g, amplitude * sin(2pi * freqHz * MTime)) + Inertia(n1, n2, Jmotor) + IdealGear(n2, n3, ratio) + Inertia(n3, n4, 2.0) + Damper(n4, g, damping) + Spring(n4, n5, 1e4) + Inertia(n5, n6, Jload) + ] end function sim_First() diff --git a/src/blocks.jl b/src/blocks.jl index ef185de..aa70327 100644 --- a/src/blocks.jl +++ b/src/blocks.jl @@ -14,28 +14,28 @@ function Integrator(u::Signal, y::Signal, k::Real) - { - der(y) - k .* u - } + @equations begin + der(y) = k .* u + end end function Integrator(u::Signal, y::Signal, k = 1.0, # Gain y_start = 0.0) # output initial value y.value = y_start - { - der(y) - k .* u - } + @equations begin + der(y) = k .* u + end end Integrator(u::Signal, y::Signal) = Integrator(u, y, 1.0) function Derivative(u::Signal, y::Signal, k::Real, T::Real) x = Unknown() # state of the block zeroGain = abs(k) < eps() - { - der(x) - (zeroGain ? 0 : (u - x) ./ T) - y - (zeroGain ? 0 : (k ./ T) .* (u - x)) - } + @equations begin + der(x) = zeroGain ? 0 : (u - x) ./ T + y = zeroGain ? 0 : (k ./ T) .* (u - x) + end end function Derivative(u::Signal, y::Signal, @@ -46,10 +46,10 @@ function Derivative(u::Signal, y::Signal, y.value = y_start x = Unknown(x_start) # state of the block zeroGain = abs(k) < eps() - { - der(x) - (zeroGain ? 0 : (u - x) ./ T) - y - (zeroGain ? 0 : (k ./ T) .* (u - x)) - } + @equations begin + der(x) = zeroGain ? 0 : (u - x) ./ T + y = zeroGain ? 0 : (k ./ T) .* (u - x) + end end Derivative(u::Signal, y::Signal; T = 1.0, @@ -62,9 +62,9 @@ function FirstOrder(u::Signal, y::Signal, k = 1.0, # Gain y_start = 0.0) # output initial value y.value = y_start - { - y + T*der(y) - k*u - } + @equations begin + y + T*der(y) = k*u + end end FirstOrder(u::Signal, y::Signal; T = 1.0, @@ -94,14 +94,14 @@ function LimPID(u_s::Signal, u_m::Signal, y::Signal, i = Unknown() # input of integrator block I = Unknown() # output of integrator block zeroGain = abs(k) < eps() - { - u_s - u_m + (y - x) / (k * Ni) - i - with_I ? Integrator(i, I, 1/Ti) : {} - with_D ? Derivative(d, D, Td, max(Td/Nd, 1e-14)) : {} - u_s - u_m - d - Limiter(x, y, yMax, yMin) - x - k * ((with_I ? I : 0.0) + (with_D ? D : 0.0) + u_s - u_m) - } + @equations begin + i = u_s - u_m + (y - x) / (k * Ni) + with_I ? Integrator(i, I, 1/Ti) : Equation[] + with_D ? Derivative(d, D, Td, max(Td/Nd, 1e-14)) : Equation[] + d = u_s - u_m + Limiter(x, y, yMax, yMin) + x = k * ((with_I ? I : 0.0) + (with_D ? D : 0.0) + u_s - u_m) + end end function LimPID(u_s::Signal, u_m::Signal, y::Signal; @@ -129,10 +129,10 @@ function StateSpace(u::Signal, y::Signal, C = [1.0], D = [0.0]) x = Unknown(zeros(size(A, 1))) # state vector - { - A * x + B * u - der(x) - C * x + D * u - y - } + @equations begin + der(x) = A * x + B * u + y = C * x + D * u + end end StateSpace(u::Signal, y::Signal; A = [1.0], @@ -156,12 +156,12 @@ function TransferFunction(u::Signal, y::Signal, if nx == 0 y - d * u else - { - der(x_scaled[1]) - (dot(-a[2:na], x_scaled) + a_end * u) / a[1] - der(x_scaled[2:nx]) - x_scaled[1:nx-1] - -y + dot(bb[2:na] - d * a[2:na], x_scaled) / a_end + d * u - x - x_scaled / a_end - } + @equations begin + der(x_scaled[1]) = (dot(-a[2:na], x_scaled) + a_end * u) / a[1] + der(x_scaled[2:nx]) = x_scaled[1:nx-1] + y = dot(bb[2:na] - d * a[2:na], x_scaled) / a_end + d * u + x = x_scaled / a_end + end end end TransferFunction(u::Signal, y::Signal, @@ -176,11 +176,11 @@ TransferFunction(u::Signal, y::Signal, ## function Limiter(u::Signal, y::Signal, uMax::Real, uMin::Real) -## { -## y - ifelse(u > uMax, uMax, -## ifelse(u < uMin, uMin, -## u)) -## } +## @equations begin +## y = ifelse(u > uMax, uMax, +## ifelse(u < uMin, uMin, +## u)) +## end ## end function Limiter(u::Signal, y::Signal, @@ -188,13 +188,12 @@ function Limiter(u::Signal, y::Signal, uMin = -uMax) clamped_pos = Discrete(false) clamped_neg = Discrete(false) - { - BoolEvent(clamped_pos, u - uMax) - BoolEvent(clamped_neg, uMin - u) - y - ifelse(clamped_pos, uMax, - ifelse(clamped_neg, uMin, - u)) - } + @equations begin + BoolEvent(clamped_pos, u - uMax) + BoolEvent(clamped_neg, uMin - u) + y = ifelse(clamped_pos, uMax, + ifelse(clamped_neg, uMin, u)) + end end Limiter(u::Signal, y::Signal; uMax = 1.0, @@ -206,12 +205,12 @@ function Step(y::Signal, offset = 0.0, startTime = 0.0) ymag = Discrete(offset) - { - y - ymag - Event(MTime - startTime, - {reinit(ymag, offset + height)}, # positive crossing - {reinit(ymag, offset)}) # negative crossing - } + @equations begin + y = ymag + Event(MTime - startTime, + Equation[reinit(ymag, offset + height)], # positive crossing + Equation[reinit(ymag, offset)]) # negative crossing + end end Step(y::Signal; height = 1.0, @@ -223,13 +222,13 @@ function DeadZone(u::Signal, y::Signal, uMin = -uMax) pos = Discrete(false) neg = Discrete(false) - { - BoolEvent(pos, u - uMax) - BoolEvent(neg, uMin - u) - y - ifelse(pos, u - uMax, - ifelse(neg, u - uMin, - 0.0)) - } + @equations begin + BoolEvent(pos, u - uMax) + BoolEvent(neg, uMin - u) + y = ifelse(pos, u - uMax, + ifelse(neg, u - uMin, + 0.0)) + end end DeadZone(u::Signal, y::Signal; uMax = 1.0, diff --git a/src/electrical.jl b/src/electrical.jl index 1538c58..7ab8962 100644 --- a/src/electrical.jl +++ b/src/electrical.jl @@ -28,17 +28,17 @@ function BranchHeatPort(n1::ElectricalNode, n2::ElectricalNode, hp::HeatPort, v = Voltage(vals) n = Voltage(vals) PowerLoss = HeatFlow(compatible_values(hp)) - { - if length(value(hp)) > 1 # hp is an array - PowerLoss - v .* i - else - PowerLoss - sum(v .* i) - end - RefBranch(hp, -PowerLoss) - n1 - n2 - v - Branch(n1, n, 0.0, i) - model(n, n2, args...) - } + @equations begin + if length(value(hp)) > 1 # hp is an array + PowerLoss = v .* i + else + PowerLoss = sum(v .* i) + end + RefBranch(hp, -PowerLoss) + v = n1 - n2 + Branch(n1, n, 0.0, i) + model(n, n2, args...) + end end @@ -50,10 +50,10 @@ end function Resistor(n1::ElectricalNode, n2::ElectricalNode, R::Signal) i = Current(compatible_values(n1, n2)) v = Voltage(value(n1) - value(n2)) - { - Branch(n1, n2, v, i) - R .* i - v # == 0 is implied - } + @equations begin + Branch(n1, n2, v, i) + v = R .* i + end end function Resistor(n1::ElectricalNode, n2::ElectricalNode, R = 1.0, T = 293.15, T_ref = 300.15, alpha = 0.0) @@ -71,10 +71,10 @@ end function Capacitor(n1::ElectricalNode, n2::ElectricalNode, C::Signal = 1.0) i = Current(compatible_values(n1, n2)) v = Voltage(value(n1) - value(n2)) - { - Branch(n1, n2, v, i) - C .* der(v) - i - } + @equations begin + Branch(n1, n2, v, i) + C .* der(v) = i + end end function Capacitor(n1::ElectricalNode, n2::ElectricalNode; C::Signal = 1.0) Capacitor(n1, n2, C) @@ -84,10 +84,10 @@ end function Inductor(n1::ElectricalNode, n2::ElectricalNode, L::Signal = 1.0) i = Current(compatible_values(n1, n2)) v = Voltage(value(n1) - value(n2)) - { - Branch(n1, n2, v, i) - L .* der(i) - v - } + @equations begin + Branch(n1, n2, v, i) + L .* der(i) = v + end end function Inductor(n1::ElectricalNode, n2::ElectricalNode; L::Signal = 1.0) Inductor(n1, n2, L) @@ -109,12 +109,12 @@ function SaturatingInductor(n1::ElectricalNode, n2::ElectricalNode, # initial equation equivalent (uses John Myles White's optim package): Ipar = optimize(Ipar -> ((Lnom - Linf) - (Lzer - Linf)*Ipar[1]/Inom*(pi/2-atan2(Ipar[1],Inom))) ^ 2, [Inom]).minimum[1] println("Ipar: ", Ipar) - { - Branch(n1, n2, v, i) - (Lact - Linf) .* i ./ Ipar - (Lzer - Linf) .* atan2(i, Ipar) - Psi - Lact .* i - der(Psi) - v - } + @equations begin + Branch(n1, n2, v, i) + (Lact - Linf) .* i ./ Ipar = (Lzer - Linf) .* atan2(i, Ipar) + Psi = Lact .* i + der(Psi) = v + end end function SaturatingInductor(n1::ElectricalNode, n2::ElectricalNode; Inom = 1.0, @@ -130,11 +130,11 @@ function SaturatingInductor2(n1::ElectricalNode, n2::ElectricalNode, i = Current(vals, "i_s") v = Voltage(value(n1) - value(n2), "v_s") psi = Unknown(vals, "psi") - { - Branch(n1, n2, v, i) - psi - a * tanh(b * i) + c * i - v - der(psi) - } + @equations begin + Branch(n1, n2, v, i) + psi = a * tanh(b * i) + c * i + v = der(psi) + end end function SaturatingInductor3(n1::ElectricalNode, n2::ElectricalNode, @@ -143,11 +143,11 @@ function SaturatingInductor3(n1::ElectricalNode, n2::ElectricalNode, i = Current(vals, "i_s") v = Voltage(value(n1) - value(n2), "v_s") psi = Unknown(vals, "psi") - { - Branch(n1, n2, v, i) - i - a * sinh(b * psi) + c * psi - v - der(psi) - } + @equations begin + Branch(n1, n2, v, i) + i = a * sinh(b * psi) - c * psi + v = der(psi) + end end function SaturatingInductor4(n1::ElectricalNode, n2::ElectricalNode, @@ -156,11 +156,11 @@ function SaturatingInductor4(n1::ElectricalNode, n2::ElectricalNode, i = Current(vals, "i_s") v = Voltage(value(n1) - value(n2), "v_s") psi = Unknown(vals, "psi") - { - Branch(n1, n2, v, i) - a * sign(i) * abs(i) ^ b + c * i - psi - v - der(psi) - } + @equations begin + Branch(n1, n2, v, i) + psi = a * sign(i) * abs(i) ^ b + c * i + der(psi) = v + end end function Transformer(p1::ElectricalNode, n1::ElectricalNode, p2::ElectricalNode, n2::ElectricalNode, @@ -171,12 +171,12 @@ function Transformer(p1::ElectricalNode, n1::ElectricalNode, p2::ElectricalNode, i1 = Current(vals) v2 = Voltage(vals) i2 = Current(vals) - { - Branch(p1, n1, v1, i1) - Branch(p2, n2, v2, i2) - L1 .* der(i1) + M .* der(i2) - v1 - M .* der(i1) + L2 .* der(i2) - v2 - } + @equations begin + Branch(p1, n1, v1, i1) + Branch(p2, n2, v2, i2) + v1 = L1 .* der(i1) + M .* der(i2) + v2 = M .* der(i1) + L2 .* der(i2) + end end function Transformer(p1::ElectricalNode, n1::ElectricalNode, p2::ElectricalNode, n2::ElectricalNode; L1 = 1.0, L2 = 1.0, M = 1.0) @@ -191,13 +191,13 @@ function EMF(n1::ElectricalNode, n2::ElectricalNode, flange::Flange, phi = Angle(compatible_values(flange, support_flange)) tau = Torque(compatible_values(flange, support_flange)) w = AngularVelocity(compatible_values(flange, support_flange)) - { - Branch(n1, n2, i, v) - Branch(flange, support_flange, phi, tau) - w - der(phi) - v - k * w - tau + k * i - } + @equations begin + Branch(n1, n2, i, v) + Branch(flange, support_flange, phi, tau) + w = der(phi) + v = k * w + tau = -k * i + end end function EMF(n1::ElectricalNode, n2::ElectricalNode, flange::Flange; support_flange = 0.0, k = 1.0) @@ -218,12 +218,12 @@ function IdealDiode(n1::ElectricalNode, n2::ElectricalNode, v = Voltage(vals) s = Unknown(vals) # dummy variable openswitch = Discrete(fill(true, length(vals))) # on/off state of each diode - { - Branch(n1, n2, v, i) - BoolEvent(openswitch, -s) # openswitch becomes true when s goes negative - s .* ifelse(openswitch, 1.0, Ron) + Vknee - v - s .* ifelse(openswitch, Goff, 1.0) + Goff .* Vknee - i - } + @equations begin + Branch(n1, n2, v, i) + BoolEvent(openswitch, -s) # openswitch becomes true when s goes negative + v = s .* ifelse(openswitch, 1.0, Ron) + Vknee + i = s .* ifelse(openswitch, Goff, 1.0) + Goff .* Vknee + end end function IdealDiode(n1::ElectricalNode, n2::ElectricalNode; Vknee = 0.0, Ron = 1e-5, Goff = 1e-5) @@ -239,12 +239,12 @@ function IdealThyristor(n1::ElectricalNode, n2::ElectricalNode, fire::Discrete, off = Discrete(true) # on/off state of each switch addhook!(fire, ifelse(fire, reinit(off, false))) - { - Branch(n1, n2, v, i) - Event(-s, reinit(off, true)) - s .* ifelse(off, 1.0, Ron) + Vknee - v - s .* ifelse(off, Goff, 1.0) + Goff .* Vknee - i - } + @equations begin + Branch(n1, n2, v, i) + Event(-s, reinit(off, true)) + v = s .* ifelse(off, 1.0, Ron) + Vknee + i = s .* ifelse(off, Goff, 1.0) + Goff .* Vknee + end end function IdealThyristor(n1::ElectricalNode, n2::ElectricalNode, fire::Discrete; Vknee = 0.0, Ron = 1e-5, Goff = 1e-5) @@ -260,12 +260,12 @@ function IdealGTOThyristor(n1::ElectricalNode, n2::ElectricalNode, fire::Discret s = Unknown(vals) # dummy variable off = Discrete(true) # on/off state of each switch addhook!(fire, reinit(off, !fire)) - { - Branch(n1, n2, v, i) - Event(-s, reinit(off, true)) - s .* ifelse(off, 1.0, Ron) + Vknee - v - s .* ifelse(off, Goff, 1.0) + Goff .* Vknee - i - } + @equations begin + Branch(n1, n2, v, i) + Event(-s, reinit(off, true)) + v = s .* ifelse(off, 1.0, Ron) + Vknee + i = s .* ifelse(off, Goff, 1.0) + Goff .* Vknee + end end function IdealGTOThyristor(n1::ElectricalNode, n2::ElectricalNode, fire::Discrete; Vknee = 0.0, Ron = 1e-5, Goff = 1e-5) @@ -276,18 +276,18 @@ end function IdealOpAmp(p1::ElectricalNode, n1::ElectricalNode, p2::ElectricalNode, n2::ElectricalNode) i = Current(compatible_values(p2, n2)) v = Voltage(compatible_values(p2, n2)) - { - p1 - n1 # voltages at the input are equal - # currents at the input are zero, so leave out - Branch(p2, n2, v, i) # at the output, make the currents equal - } + @equations begin + p1 = n1 # voltages at the input are equal + # currents at the input are zero, so leave out + Branch(p2, n2, v, i) # at the output, make the currents equal + end end function IdealOpAmp(p1::ElectricalNode, n1::ElectricalNode, p2::ElectricalNode) i = Current(compatible_values(p2)) - { - p1 - n1 - RefBranch(p2, i) - } + @equations begin + p1 = n1 + RefBranch(p2, i) + end end function IdealOpeningSwitch(n1::ElectricalNode, n2::ElectricalNode, control::Discrete, @@ -296,11 +296,11 @@ function IdealOpeningSwitch(n1::ElectricalNode, n2::ElectricalNode, control::Dis i = Current(vals) v = Voltage(vals) s = Unknown(vals) # dummy variable - { - Branch(n1, n2, v, i) - v - s .* ifelse(control, 1.0, Ron) - i - s .* ifelse(control, Goff, 1.0) - } + @equations begin + Branch(n1, n2, v, i) + v = s .* ifelse(control, 1.0, Ron) + i = s .* ifelse(control, Goff, 1.0) + end end function IdealOpeningSwitch(n1::ElectricalNode, n2::ElectricalNode, control::Discrete; Ron = 1e-5, Goff = 1e-5) @@ -313,11 +313,11 @@ function IdealClosingSwitch(n1::ElectricalNode, n2::ElectricalNode, control::Dis i = Current(vals) v = Voltage(vals) s = Unknown(vals) # dummy variable - { - Branch(n1, n2, v, i) - v - s .* ifelse(control, Ron, 1.0) - i - s .* ifelse(control, 1.0, Goff) - } + @equations begin + Branch(n1, n2, v, i) + v = s .* ifelse(control, Ron, 1.0) + i = s .* ifelse(control, 1.0, Goff) + end end function IdealClosingSwitch(n1::ElectricalNode, n2::ElectricalNode, control::Discrete; Ron = 1e-5, Goff = 1e-5) @@ -331,12 +331,12 @@ function ControlledIdealOpeningSwitch(n1::ElectricalNode, n2::ElectricalNode, co v = Voltage(vals) s = Unknown(vals) # dummy variable openswitch = Discrete(fill(true, length(vals))) # on/off state of diode - { - Branch(n1, n2, v, i) - BoolEvent(openswitch, control - level) # openswitch becomes false when control goes below level - s .* ifelse(openswitch, 1.0, Ron) - v - s .* ifelse(openswitch, Goff, 1.0) - i - } + @equations begin + Branch(n1, n2, v, i) + BoolEvent(openswitch, control - level) # openswitch becomes false when control goes below level + v = s .* ifelse(openswitch, 1.0, Ron) + i = s .* ifelse(openswitch, Goff, 1.0) + end end function ControlledIdealOpeningSwitch(n1::ElectricalNode, n2::ElectricalNode, control::Signal; level = 0.0, Ron = 1e-5, Goff = 1e-5) @@ -362,32 +362,32 @@ function ControlledOpenerWithArc(n1::ElectricalNode, n2::ElectricalNode, control quenched = Discrete(true) # whether the arc is quenched or not tSwitch = Discrete(0.0) # time of last open initiation ipositive = Discrete(true) # whether the current is positive - { - Branch(n1, n2, v, i) - Event(level - control, - reinit(on, true), - { - reinit(on, false) - reinit(quenched, false) - reinit(tSwitch, MTime) - }) - Event(i, - { - reinit(i, 0.0) - reinit(ipositive, true) - ifelse(!quenched, reinit(quenched, true)) - }, - { - reinit(i, 0.0) - reinit(ipositive, false) - ifelse(!quenched, reinit(quenched, true)) - }) - ifelse(on, - v - Ron .* i, - ifelse(quenched, - i - Goff .* v, - v - min(Vmax, V0 + dVdt .* (MTime - tSwitch))) .* sign(i)) - } + @equations begin + Branch(n1, n2, v, i) + Event(level - control, + reinit(on, true), + Equation[ + reinit(on, false) + reinit(quenched, false) + reinit(tSwitch, MTime) + ]) + Event(i, + Equation[ + reinit(i, 0.0) + reinit(ipositive, true) + ifelse(!quenched, reinit(quenched, true)) + ], + Equation[ + reinit(i, 0.0) + reinit(ipositive, false) + ifelse(!quenched, reinit(quenched, true)) + ]) + ifelse(on, + v - Ron .* i, + ifelse(quenched, + i - Goff .* v, + v - min(Vmax, V0 + dVdt .* (MTime - tSwitch))) .* sign(i)) + end end function ControlledOpenerWithArc(n1::ElectricalNode, n2::ElectricalNode, control::Signal; level = 0.0, Ron = 1e-5, Goff = 1e-5, V0 = 30.0, dVdt = 10e3, Vmax = 60.0) @@ -407,12 +407,12 @@ function Diode(n1::ElectricalNode, n2::ElectricalNode, vals = compatible_values(n1, n2) i = Current(vals) v = Voltage(vals) - { - Branch(n1, n2, v, i) - i - ifelse(v ./ Vt > Maxexp, - Ids .* exp(Maxexp) .* (1 + v ./ Vt - Maxexp) - 1 + v ./ R, - Ids .* (exp(v ./ Vt) - 1) + v ./ R) - } + @equations begin + Branch(n1, n2, v, i) + i = ifelse(v ./ Vt > Maxexp, + Ids .* exp(Maxexp) .* (1 + v ./ Vt - Maxexp) - 1 + v ./ R, + Ids .* (exp(v ./ Vt) - 1) + v ./ R) + end end function Diode(n1::ElectricalNode, n2::ElectricalNode; Ids = 1e-6, Vt = 0.04, Maxexp = 15, R = 1e8) @@ -429,14 +429,14 @@ function ZDiode(n1::ElectricalNode, n2::ElectricalNode, vals = compatible_values(n1, n2) i = Current(vals) v = Voltage(vals) - { - Branch(n1, n2, v, i) - i - ifelse(v ./ Vt > Maxexp, - Ids .* exp(Maxexp) .* (1 + v ./ Vt - Maxexp) - 1 + v ./ R, - ifelse((v + Bv) < -Maxexp .* (Nbv .* Vt), - -Ids - Ibv .* exp(Maxexp) .* (1 - (v+Bv) ./ (Nbv .* Vt) - Maxexp) + v ./ R, - Ids .* (exp(v ./ Vt)-1) - Ibv .* exp(-(v + Bv)/(Nbv .* Vt)) + v ./ R)) - } + @equations begin + Branch(n1, n2, v, i) + i = ifelse(v ./ Vt > Maxexp, + Ids .* exp(Maxexp) .* (1 + v ./ Vt - Maxexp) - 1 + v ./ R, + ifelse((v + Bv) < -Maxexp .* (Nbv .* Vt), + -Ids - Ibv .* exp(Maxexp) .* (1 - (v+Bv) ./ (Nbv .* Vt) - Maxexp) + v ./ R, + Ids .* (exp(v ./ Vt)-1) - Ibv .* exp(-(v + Bv)/(Nbv .* Vt)) + v ./ R)) + end end function ZDiode(n1::ElectricalNode, n2::ElectricalNode; Ids = 1e-6, Vt = 0.04, Maxexp = 30.0, R = 1e8, Bv = 5.1, Ibv = 0.7, Nbv = 0.74) @@ -456,12 +456,12 @@ function HeatingDiode(n1::ElectricalNode, n2::ElectricalNode, v = Voltage(vals) k = 1.380662e-23 # Boltzmann's constant, J/K q = 1.6021892e-19 # Electron charge, As - { - Branch(n1, n2, v, i) - i - ifelse(v ./ Vt > Maxexp, - Ids .* exp(Maxexp) .* (1 + v ./ Vt - Maxexp) - 1 + v ./ R, - Ids .* (exp(v ./ Vt) - 1) + v ./ R) - } + @equations begin + Branch(n1, n2, v, i) + i = ifelse(v ./ Vt > Maxexp, + Ids .* exp(Maxexp) .* (1 + v ./ Vt - Maxexp) - 1 + v ./ R, + Ids .* (exp(v ./ Vt) - 1) + v ./ R) + end end function HeatingDiode(n1::ElectricalNode, n2::ElectricalNode; T = 293.15, Ids = 1e-6, Maxexp = 15, R = 1e8, EG = 1.11, N = 1.0, TNOM = 300.15, XTI = 3.0) @@ -478,10 +478,10 @@ end function SignalVoltage(n1::ElectricalNode, n2::ElectricalNode, V::Signal) i = Current(compatible_values(n1, n2)) v = Voltage(compatible_values(n1, n2)) - { - Branch(n1, n2, v, i) - v - V - } + @equations begin + Branch(n1, n2, v, i) + v = V + end end function SineVoltage(n1::ElectricalNode, n2::ElectricalNode, @@ -498,13 +498,13 @@ function StepVoltage(n1::ElectricalNode, n2::ElectricalNode, i = Current(compatible_values(n1, n2)) v = Voltage(compatible_values(n1, n2)) v_mag = Discrete(offset) - { - Branch(n1, n2, v, i) - v - v_mag - Event(MTime - start, - {reinit(v_mag, offset + V)}, # positive crossing - {reinit(v_mag, offset)}) # negative crossing - } + @equations begin + Branch(n1, n2, v, i) + v = v_mag + Event(MTime - start, + Equation[reinit(v_mag, offset + V)], # positive crossing + Equation[reinit(v_mag, offset)]) # negative crossing + end end function StepVoltage(n1::ElectricalNode, n2::ElectricalNode; V = 1.0, start = 0.0, offset = 0.0) @@ -514,8 +514,8 @@ end function SignalCurrent(n1::ElectricalNode, n2::ElectricalNode, I::Signal) - { - RefBranch(n1, I) - RefBranch(n2, -I) - } + @equations begin + RefBranch(n1, I) + RefBranch(n2, -I) + end end diff --git a/src/heat_transfer.jl b/src/heat_transfer.jl index 81fc229..a2598f0 100644 --- a/src/heat_transfer.jl +++ b/src/heat_transfer.jl @@ -13,10 +13,10 @@ function HeatCapacitor(hp::HeatPort, C::Signal, Tstart::Real) # can't really use Tstart here. Need to define at the top level. Q_flow = HeatFlow(compatible_values(hp)) - { - RefBranch(hp, Q_flow) - C .* der(hp) - Q_flow - } + @equations begin + RefBranch(hp, Q_flow) + C .* der(hp) = Q_flow + end end HeatCapacitor(hp::HeatPort, C::Signal) = HeatCapacitor(hp, C, 293.15) @@ -24,39 +24,39 @@ HeatCapacitor(hp::HeatPort, C::Signal) = HeatCapacitor(hp, C, 293.15) function ThermalConductor(port_a::HeatPort, port_b::HeatPort, G::Signal) dT = Temperature(compatible_values(port_a, port_b)) Q_flow = HeatFlow(compatible_values(port_a, port_b)) - { - Branch(port_a, port_b, dT, Q_flow) - G .* dT - Q_flow - } + @equations begin + Branch(port_a, port_b, dT, Q_flow) + Q_flow = G .* dT + end end function Convection(port_a::HeatPort, port_b::HeatPort, Gc::Signal) dT = Temperature(compatible_values(port_a, port_b)) Q_flow = HeatFlow(compatible_values(port_a, port_b)) - { - Branch(port_a, port_b, dT, Q_flow) - Gc .* dT - Q_flow - } + @equations begin + Branch(port_a, port_b, dT, Q_flow) + Q_flow = Gc .* dT + end end function BodyRadiation(port_a::HeatPort, port_b::HeatPort, Gr::Signal) Q_flow = HeatFlow(compatible_values(port_a, port_b)) - { - RefBranch(port_a, Q_flow) - RefBranch(port_b, -Q_flow) - sigma .* Gr .* (port_a .^ 4 - port_b .^ 4) - Q_flow - } + @equations begin + RefBranch(port_a, Q_flow) + RefBranch(port_b, -Q_flow) + Q_flow = sigma .* Gr .* (port_a .^ 4 - port_b .^ 4) + end end function ThermalCollector(port_a::HeatPort, port_b::HeatPort) # This ends up being a short circuit. Q_flow = HeatFlow(compatible_values(port_a, port_b)) - { - Branch(port_a, port_b, 0.0, Q_flow) - } + @equations begin + Branch(port_a, port_b, 0.0, Q_flow) + end end @@ -68,18 +68,18 @@ end function FixedTemperature(port::HeatPort, T::Signal) Q_flow = HeatFlow(compatible_values(port, T)) - { - Branch(port, T, 0.0, Q_flow) - } + @equations begin + Branch(port, T, 0.0, Q_flow) + end end PrescribedTemperature = FixedTemperature function FixedHeatFlow(port::HeatPort, Q_flow::Signal, T_ref::Signal, alpha::Signal) Q_flow = HeatFlow(compatible_values(port, T)) - { - RefBranch(port, Q_flow .* alpha .* (port - T_ref)) - } + @equations begin + RefBranch(port, Q_flow .* alpha .* (port - T_ref)) + end end FixedHeatFlow(port::HeatPort, Q_flow::Signal, T_ref::Signal) = FixedHeatFlow(port, Q_flow, T_ref, 0.0) FixedHeatFlow(port::HeatPort, Q_flow::Signal) = FixedHeatFlow(port, Q_flow, 293.15, 0.0) diff --git a/src/machines.jl b/src/machines.jl index fe44d1f..9416554 100644 --- a/src/machines.jl +++ b/src/machines.jl @@ -72,31 +72,31 @@ function AIMSquirrelCage(n1::ElectricalNode, n2::ElectricalNode, flange::Flange, gap_flange = Angle(compatible_values(flange)) m = 3 gnd = 0.0 - { - SpacePhasorConvertor(n4, n2, sp_in, zero, 1.0) - SquirrelCage(sp_gap2, - T, - Lrsigma, # Rotor stray inductance (equivalent three-phase winding) - Rr, # Rotor resistance (equivalent three-phase winding) - TrRef, # Reference temperature of rotor resistance - convertAlpha(alpha20r, TrRef)) - AirGapS(sp_gap1, sp_gap2, gap_flange, support, - Lm, m, p) - Inertia(gap_flange, gnd, Jr) # inertiaRotor - if isa(support, Angle) - Inertia(support, gnd, Js) # inertiaStator - end - Resistor(n3, n4, - fill(Rs, m), - T, - fill(TsRef, m), - fill(convertAlpha(alpha20s, TsRef), m)) - Inductor(sp_in, sp_gap1, fill(Lssigma, 2)) - Inductor(zero, gnd, Lszero) - Core(sp_in, T, statorCoreParameters) - Friction(flange, gnd, T, frictionParameters) - StrayLoad(n1, n3, flange, gnd, T, strayLoadParameters) - } + @equations begin + SpacePhasorConvertor(n4, n2, sp_in, zero, 1.0) + SquirrelCage(sp_gap2, + T, + Lrsigma, # Rotor stray inductance (equivalent three-phase winding) + Rr, # Rotor resistance (equivalent three-phase winding) + TrRef, # Reference temperature of rotor resistance + convertAlpha(alpha20r, TrRef)) + AirGapS(sp_gap1, sp_gap2, gap_flange, support, + Lm, m, p) + Inertia(gap_flange, gnd, Jr) # inertiaRotor + if isa(support, Angle) + Inertia(support, gnd, Js) # inertiaStator + end + Resistor(n3, n4, + fill(Rs, m), + T, + fill(TsRef, m), + fill(convertAlpha(alpha20s, TsRef), m)) + Inductor(sp_in, sp_gap1, fill(Lssigma, 2)) + Inductor(zero, gnd, Lszero) + Core(sp_in, T, statorCoreParameters) + Friction(flange, gnd, T, frictionParameters) + StrayLoad(n1, n3, flange, gnd, T, strayLoadParameters) + end end function AIMSquirrelCage(n1::ElectricalNode, n2::ElectricalNode, flange::Flange, support::Flange, T::HeatPort; @@ -147,31 +147,31 @@ function AirGapS(sp_s::ElectricalNode, sp_r::ElectricalNode, flange::Flange, sup tauElectrical = Torque(compatible_values(flange, support)) RotationMatrix = Unknown(zeros(2,2)) L = [Lm 0;0 Lm] - { - # mechanical angle of the rotor of an equivalent 2-pole machine - gamma - p*(flange-support) - RotationMatrix - [+cos(gamma) -sin(gamma); +sin(gamma) +cos(gamma)] - RefBranch(sp_s, i_ss) - i_ss - RotationMatrix * i_sr - RefBranch(sp_r, i_rr) - i_rs - RotationMatrix * i_rr - RefBranch(flange, flange_tau) - RefBranch(support, support_tau) - # Magnetizing current with respect to the stator reference frame - i_ms - i_ss - i_rs - # Magnetizing flux linkage with respect to the stator reference frame - psi_ms - L*i_ms; - # Magnetizing flux linkage with respect to the rotor reference frame - psi_mr - RotationMatrix' * psi_ms - # Stator voltage induction - sp_s - der(psi_ms) - # Rotor voltage induction - sp_r - der(psi_mr) - # Electromechanical torque (cross product of current and flux space phasor) - tauElectrical - m/2*p*(i_ss[2]*psi_ms[1] - i_ss[1]*psi_ms[2]) - flange_tau + tauElectrical - support_tau - tauElectrical - } + @equations begin + # mechanical angle of the rotor of an equivalent 2-pole machine + gamma - p*(flange-support) + RotationMatrix - [+cos(gamma) -sin(gamma); +sin(gamma) +cos(gamma)] + RefBranch(sp_s, i_ss) + i_ss - RotationMatrix * i_sr + RefBranch(sp_r, i_rr) + i_rs - RotationMatrix * i_rr + RefBranch(flange, flange_tau) + RefBranch(support, support_tau) + # Magnetizing current with respect to the stator reference frame + i_ms - i_ss - i_rs + # Magnetizing flux linkage with respect to the stator reference frame + psi_ms - L*i_ms; + # Magnetizing flux linkage with respect to the rotor reference frame + psi_mr - RotationMatrix' * psi_ms + # Stator voltage induction + sp_s - der(psi_ms) + # Rotor voltage induction + sp_r - der(psi_mr) + # Electromechanical torque (cross product of current and flux space phasor) + tauElectrical - m/2*p*(i_ss[2]*psi_ms[1] - i_ss[1]*psi_ms[2]) + flange_tau + tauElectrical + support_tau - tauElectrical + end end function quasiRMS(x) # Calculate quasi-RMS value of input @@ -193,21 +193,21 @@ function StrayLoad(n1::ElectricalNode, n2::ElectricalNode, flange::Flange, suppo w = AngularVelocity(compatible_values(flange, support)) tau = Torque(compatible_values(flange, support)) m = length(n1) - { - Branch(n1, n2, v, i) - Branch(flange, support, phi, tau) - RefBranch(T, -powerLoss) - v - if isa(T, Temperature) - powerLoss - tau .* w - end - tau - ifelse(strayLoadParameters.PRef <= 0.0, - 0.0, - -strayLoadParameters.tauRef * (quasiRMS(i) / strayLoadParameters.IRef) ^ 2 * - ifelse(w >= 0.0, - +(w/strayLoadParameters.wRef)^strayLoadParameters.power_w, - -(-w/strayLoadParameters.wRef)^strayLoadParameters.power_w)) - } + @equations begin + Branch(n1, n2, v, i) + Branch(flange, support, phi, tau) + RefBranch(T, -powerLoss) + v + if isa(T, Temperature) + powerLoss - tau .* w + end + tau - ifelse(strayLoadParameters.PRef <= 0.0, + 0.0, + -strayLoadParameters.tauRef * (quasiRMS(i) / strayLoadParameters.IRef) ^ 2 * + ifelse(w >= 0.0, + +(w/strayLoadParameters.wRef)^strayLoadParameters.power_w, + -(-w/strayLoadParameters.wRef)^strayLoadParameters.power_w)) + end end function Core(sp::ElectricalNode, T::HeatPort, coreParameters::CoreParameters) @@ -215,17 +215,17 @@ function Core(sp::ElectricalNode, T::HeatPort, coreParameters::CoreParameters) sp_i = Current(vals) powerLoss = HeatFlow(compatible_values(T)) Gc = Unknown(vals) - { - RefBranch(sp, sp_i) - RefBranch(T, -powerLoss) - Gc - ifelse(coreParameters.PRef <= 0.0, - 0.0, - coreParameters.GcRef) - sp_i - Gc .* sp - if isa(T, Temperature) - powerLoss + 3/2 * (sp[1]*sp_i[1] + sp[2]*sp_i[2]) - end - } + @equations begin + RefBranch(sp, sp_i) + RefBranch(T, -powerLoss) + Gc - ifelse(coreParameters.PRef <= 0.0, + 0.0, + coreParameters.GcRef) + sp_i - Gc .* sp + if isa(T, Temperature) + powerLoss + 3/2 * (sp[1]*sp_i[1] + sp[2]*sp_i[2]) + end + end end function Friction(flange::Flange, support::Flange, T::HeatPort, frictionParameters::FrictionParameters) @@ -233,21 +233,21 @@ function Friction(flange::Flange, support::Flange, T::HeatPort, frictionParamete w = AngularVelocity(compatible_values(flange, support)) phi = Angle(compatible_values(flange, support)) tau = Torque(compatible_values(flange, support)) - { - Branch(flange, support, phi, tau) - RefBranch(T, -powerLoss) - w - der(phi) - if isa(T, Temperature) - powerLoss - tau .* w - end - tau - ifelse(frictionParameters.PRef <= 0.0, - 0.0, - ifelse(w >= frictionParameters.wLinear, - frictionParameters.tauRef*(+w/frictionParameters.wRef)^frictionParameters.power_w, - ifelse(w <= -frictionParameters.wLinear, - -frictionParameters.tauRef*(-w/frictionParameters.wRef)^frictionParameters.power_w, - frictionParameters.tauLinear*(w/frictionParameters.wLinear)))) - } + @equations begin + Branch(flange, support, phi, tau) + RefBranch(T, -powerLoss) + w - der(phi) + if isa(T, Temperature) + powerLoss - tau .* w + end + tau - ifelse(frictionParameters.PRef <= 0.0, + 0.0, + ifelse(w >= frictionParameters.wLinear, + frictionParameters.tauRef*(+w/frictionParameters.wRef)^frictionParameters.power_w, + ifelse(w <= -frictionParameters.wLinear, + -frictionParameters.tauRef*(-w/frictionParameters.wRef)^frictionParameters.power_w, + frictionParameters.tauLinear*(w/frictionParameters.wLinear)))) + end end function SquirrelCage(spacePhasor_r::ElectricalNode, T::HeatPort, Lrsigma, Rr, T_ref, alpha) @@ -255,15 +255,15 @@ function SquirrelCage(spacePhasor_r::ElectricalNode, T::HeatPort, Lrsigma, Rr, T spacePhasor_i = Current(vals) powerLoss = HeatFlow(compatible_values(T)) Rr_actual = Unknown(compatible_values(Rr)) - { - RefBranch(spacePhasor_r, spacePhasor_i) - RefBranch(T, -powerLoss) - Rr_actual - Rr .* (1 + alpha .* (T - T_ref)) - spacePhasor_r - Rr_actual .* spacePhasor_i - Lrsigma .* der(spacePhasor_i) - if isa(T, Temperature) - 2/3*powerLoss - Rr_actual .* (spacePhasor_i[1]*spacePhasor_i[1] + spacePhasor_i[2]*spacePhasor_i[2]) - end - } + @equations begin + RefBranch(spacePhasor_r, spacePhasor_i) + RefBranch(T, -powerLoss) + Rr_actual - Rr .* (1 + alpha .* (T - T_ref)) + spacePhasor_r - Rr_actual .* spacePhasor_i - Lrsigma .* der(spacePhasor_i) + if isa(T, Temperature) + 2/3*powerLoss - Rr_actual .* (spacePhasor_i[1]*spacePhasor_i[1] + spacePhasor_i[2]*spacePhasor_i[2]) + end + end end function SpacePhasorConvertor(n1::ElectricalNode, n2::ElectricalNode, @@ -279,19 +279,19 @@ function SpacePhasorConvertor(n1::ElectricalNode, n2::ElectricalNode, const m = 3 TransformationMatrix = 2/m * [[cos(+(k - 1)/m*2*pi) for k in 1:m] [+sin(+(k - 1)/m*2*pi) for k in 1:m]]' ## InverseTransformation = [[cos(-(k - 1)/m*2*pi), -sin(-(k - 1)/m*2*pi)]' for k in 1:m] - { - RefBranch(n1, i1) - RefBranch(n2, i2) - RefBranch(zero_v, zero_i) - RefBranch(spacePhasor, spacePhasor_i) - v / turnsRatio - n1 + n2 - i*turnsRatio - i1 - i*turnsRatio + i2 - m*zero_v - sum(v) - spacePhasor - TransformationMatrix * v - -m*zero_i - sum(i) - -spacePhasor_i - TransformationMatrix * i - } + @equations begin + RefBranch(n1, i1) + RefBranch(n2, i2) + RefBranch(zero_v, zero_i) + RefBranch(spacePhasor, spacePhasor_i) + v / turnsRatio - n1 + n2 + i*turnsRatio - i1 + i*turnsRatio + i2 + m*zero_v - sum(v) + spacePhasor - TransformationMatrix * v + -m*zero_i - sum(i) + -spacePhasor_i - TransformationMatrix * i + end end function AIMC_DOL() @@ -308,14 +308,14 @@ function AIMC_DOL() r2 = Angle(0.0, "LoadAngle") control = Discrete(false) gnd = 0.0 - { - SineVoltage(n1, gnd, VAC .* sqrt(2/3), fNominal, [0.0, -2pi/3, 2pi/3]) - AIMSquirrelCage(n3, gnd, r1) - Inertia(r1, r2, JLoad) - QuadraticSpeedDependentTorque(r2, gnd, -TLoad, false, wLoad) - BoolEvent(control, MTime - tStart1) - IdealClosingSwitch(n1, n2, control) - } + @equations begin + SineVoltage(n1, gnd, VAC .* sqrt(2/3), fNominal, [0.0, -2pi/3, 2pi/3]) + AIMSquirrelCage(n3, gnd, r1) + Inertia(r1, r2, JLoad) + QuadraticSpeedDependentTorque(r2, gnd, -TLoad, false, wLoad) + BoolEvent(control, MTime - tStart1) + IdealClosingSwitch(n1, n2, control) + end end function sim_AIMC_DOL() diff --git a/src/main.jl b/src/main.jl index f279387..4fdda72 100644 --- a/src/main.jl +++ b/src/main.jl @@ -631,11 +631,7 @@ parse_args(a::Array) = [parse_args(x) for x in a] parse_args(x) = x function equations_helper(arg) - if arg.head == :block - Expr(:ref, :Equation, parse_args(arg.args)...) - else - error("Argument to @equations is not a begin..end block") - end + Expr(:ref, :Equation, parse_args(arg.args)...) end -Equation = Any +typealias Equation Any diff --git a/src/powersystems.jl b/src/powersystems.jl index 229bd40..d7b8507 100644 --- a/src/powersystems.jl +++ b/src/powersystems.jl @@ -20,10 +20,10 @@ function RLLine(n1::ElectricalNode, n2::ElectricalNode, Z::SeriesImpedance, len: v = Voltage(vals) R = real(Z) * len L = imag(Z) * len / (2pi * freq) - { - Branch(n1, n2, v, i) - L * der(i) + R * i - v - } + @equations begin + Branch(n1, n2, v, i) + v = L * der(i) + R * i + end end @@ -37,12 +37,12 @@ function PiLine(n1::ElectricalNode, n2::ElectricalNode, Z::SeriesImpedance, Y::S L = imag(Z) / ne * len / (2pi * freq) G = real(Z) * ne / len C = imag(Z) * ne / len / (2pi * freq) - { - RefBranch(n1, i[:, 1]) - RefBranch(n2, -i[:, ne1]) - C * der(v) + G * v - (i[:, 1:ne] - i[:, 2:ne1]) - L * der(i) + R * i - [[2 * (n1 - v[:, 1])] v[:, 1:ne - 1] - v[:, 2:ne] [2 * (v[:, ne] - n2)]] - } + @equations begin + RefBranch(n1, i[:, 1]) + RefBranch(n2, -i[:, ne1]) + C * der(v) + G * v = i[:, 1:ne] - i[:, 2:ne1] + L * der(i) + R * i = [[2 * (n1 - v[:, 1])] v[:, 1:ne - 1] - v[:, 2:ne] [2 * (v[:, ne] - n2)]] + end end function ModalLine(v1::ElectricalNode, v2::ElectricalNode, Z::SeriesImpedance, Y::ShuntAdmittance, len::Real, freq::Real) @@ -72,26 +72,26 @@ function ModalLine(v1::ElectricalNode, v2::ElectricalNode, Z::SeriesImpedance, Y Tv = real(Tv) # equations - { - RefBranch(v1, i1) - RefBranch(v2, -i2) - Imode1 - Tv.' * i1 - Imode2 - Tv.' * i2 - Vmode1 - Ti.' * v1 - Vmode2 - Ti.' * v2 - ## i1 - Ti * Imode1 - ## i2 - Ti * Imode2 - ## v1 - Tv * Vmode1 - ## v2 - Tv * Vmode2 - Imode1 - Iout1 + delay(Iout2, delays) - Imode2 - Iout2 + delay(Iout1, delays) - Vmode1 ./ Zsurge - Iout1 - delay(Iout2, delays) - Vmode2 ./ Zsurge - Iout2 - delay(Iout1, delays) - ## Imode1 -> <- delay(Iout2) delay(Iout1) -> <- Imode2 - ## | | - ## \/ \/ - ## Iout1 Iout2 - } + @equations begin + RefBranch(v1, i1) + RefBranch(v2, -i2) + Imode1 = Tv.' * i1 + Imode2 = Tv.' * i2 + Vmode1 = Ti.' * v1 + Vmode2 = Ti.' * v2 + ## i1 = Ti * Imode1 + ## i2 = Ti * Imode2 + ## v1 = Tv * Vmode1 + ## v2 = Tv * Vmode2 + Imode1 = Iout1 + delay(Iout2, delays) + Imode2 = Iout2 + delay(Iout1, delays) + Vmode1 ./ Zsurge = Iout1 + delay(Iout2, delays) + Vmode2 ./ Zsurge = Iout2 + delay(Iout1, delays) + ## Imode1 -> <- delay(Iout2) delay(Iout1) -> <- Imode2 + ## | | + ## \/ \/ + ## Iout1 Iout2 + end end @@ -208,10 +208,10 @@ function ConstZParallelLoad(n1::ElectricalNode, n2::ElectricalNode, if load_pf < 1.0 L = baseZ / sqrt(1 - load_pf .^ 2) / (2pi .* freq) end - { - Resistor(n1, n2, R) - load_pf < 1.0 ? Inductor(n1, n2, L) : {} - } + @equations begin + Resistor(n1, n2, R) + load_pf < 1.0 ? Inductor(n1, n2, L) : Equation[] + end ## This seems broken when used in the example below. ## With pf=1, it takes a long time (stiff?). With ## pf = 0.85, there's a large offset. @@ -225,8 +225,8 @@ function ConstZSeriesLoad(n1::ElectricalNode, n2::ElectricalNode, baseZ = (Vln .^ 2) / load_VA R = baseZ * load_pf L = baseZ * sqrt(1 - load_pf .^ 2) / (2pi .* freq) - { - Branch(n1, n2, v, i) - L .* der(i) + R .* i - v - } + @equations begin + Branch(n1, n2, v, i) + v = L .* der(i) + R .* i + end end diff --git a/src/rotational.jl b/src/rotational.jl index 683911a..15562a8 100644 --- a/src/rotational.jl +++ b/src/rotational.jl @@ -19,23 +19,23 @@ function Inertia(flange_a::Flange, flange_b::Flange, J::Real) tau_b = Torque(val) w = AngularVelocity(val) a = AngularAcceleration(val) - { - RefBranch(flange_a, tau_a) - RefBranch(flange_b, tau_b) - flange_b - flange_a # the angles are both equal - w - der(flange_a) - a - der(w) - tau_a + tau_b - J .* a - } + @equations begin + RefBranch(flange_a, tau_a) + RefBranch(flange_b, tau_b) + flange_b = flange_a # the angles are both equal + w = der(flange_a) + a = der(w) + tau_a + tau_b = J .* a + end end function Disc(flange_a::Flange, flange_b::Flange, deltaPhi::Signal) "1-dim. rotational rigid component without inertia, where right flange is rotated by a fixed angle with respect to left flange" tau = Torque(compatible_values(flange_a, flange_b)) - { - RefBranch(flange_b, flange_a, deltaPhi, tau) - } + @equations begin + RefBranch(flange_b, flange_a, deltaPhi, tau) + end end @@ -44,10 +44,10 @@ function Spring(flange_a::Flange, flange_b::Flange, c::Real, phi_rel0::Signal) val = compatible_values(flange_a, flange_b) phi_rel = Angle(val) tau = Torque(val) - { - Branch(flange_b, flange_a, phi_rel, tau) - tau - c .* (phi_rel - phi_rel0) - } + @equations begin + Branch(flange_b, flange_a, phi_rel, tau) + tau = c .* (phi_rel - phi_rel0) + end end Spring(flange_a::Flange, flange_b::Flange, c::Real) = Spring(flange_a, flange_b, c, 0.0) @@ -59,18 +59,18 @@ function BranchHeatPort(n1::Flange, n2::Flange, hp::HeatPort, w_rel = AngularVelocity(val) tau = Torque(val) PowerLoss = Power(compatible_values(hp)) - { - n1 - n2 - phi_rel - w_rel - der(phi_rel) - if length(value(hp)) > 1 # an array - PowerLoss - w_rel .* tau - else - PowerLoss - sum(w_rel .* tau) - end - RefBranch(hp, -PowerLoss) - Branch(n1, n, 0.0, tau) - model(n, n2, args...) - } + @equations begin + n1 - n2 = phi_rel + w_rel = der(phi_rel) + if length(value(hp)) > 1 # an array + PowerLoss = w_rel .* tau + else + PowerLoss = sum(w_rel .* tau) + end + RefBranch(hp, -PowerLoss) + Branch(n1, n, 0.0, tau) + model(n, n2, args...) + end end @@ -79,10 +79,10 @@ function Damper(flange_a::Flange, flange_b::Flange, d::Signal) val = compatible_values(flange_a, flange_b) phi_rel = Angle(val) tau = Torque(val) - { - Branch(flange_b, flange_a, phi_rel, tau) - tau - d .* der(phi_rel) - } + @equations begin + Branch(flange_b, flange_a, phi_rel, tau) + tau = d .* der(phi_rel) + end end Damper(flange_a::Flange, flange_b::Flange, hp::HeatPort, d::Signal) = BranchHeatPort(flange_a, flange_b, hp, Damper, d) @@ -93,10 +93,10 @@ function SpringDamper(flange_a::Flange, flange_b::Flange, c::Signal, d::Signal) val = compatible_values(flange_a, flange_b) phi_rel = Angle(val) tau = Torque(val) - { - Spring(flange_a, flange_b, c) - Damper(flange_a, flange_b, d) - } + @equations begin + Spring(flange_a, flange_b, c) + Damper(flange_a, flange_b, d) + end end SpringDamper(flange_a::Flange, flange_b::Flange, hp::HeatPort, c::Signal, d::Signal) = BranchHeatPort(flange_a, flange_b, hp, SpringDamper, c, d) @@ -123,81 +123,81 @@ function Brake(flange_a::Flange, flange_b::Flange, support::Flange, f_normalized const Stuck=0 # w_rel = 0 (forward sliding, locked or backward sliding) const Backward=-1 # w_rel < 0 (backward sliding) mode = Discrete(fill(UnknownMode, length(vals))) - { - RefBranch(flange_a, tau_a) - RefBranch(flange_b, tau_b) - - phi - flange_a + support - flange_b - flange_a + @equations begin + RefBranch(flange_a, tau_a) + RefBranch(flange_b, tau_b) + + phi - flange_a + support + flange_b - flange_a - # Angular velocity and angular acceleration of flanges flange_a and flange_b - w - der(phi) - a - der(w) - w_relfric - w - a_relfric - a - - # Friction torque, normal force and friction torque for w_rel=0 - tau_a + tau_b - tau - 0 - fn - fn_max .* f_normalized - tau0 - mue0 .* cgeo .* fn - tau0_max - peak .* tau0 - BoolEvent(free, fn) - Event(w, - { - reinit(startForward, - pre(mode) == Stuck & (sa > tau0_max/unitTorque | pre(startForward)) & - sa > tau0/unitTorque | pre(mode) == Backward & w > w_small | initial() & (w > 0)) - reinit(startBackward, - pre(mode) == Stuck & (sa < -tau0_max/unitTorque | pre(startBackward) & - sa < -tau0/unitTorque) | pre(mode) == Forward & w < -w_small | initial() & (w < 0)) - reinit(locked, - !free && !(pre(mode) == Forward | startForward | pre(mode) == Backward | startBackward)) - # finite state machine to determine configuration - reinit(mode, - ifelse(free, Free, - ifelse((pre(mode) == Forward | pre(mode) == Free | startForward) & w > 0, Forward, - ifelse((pre(mode) == Backward | pre(mode) == Free | startBackward) & w < 0, Backward, - Stuck)))) - }) + # Angular velocity and angular acceleration of flanges flange_a and flange_b + w = der(phi) + a = der(w) + w_relfric = w + a_relfric = a + + # Friction torque, normal force and friction torque for w_rel=0 + tau_a + tau_b = tau + fn = fn_max .* f_normalized + tau0 = mue0 .* cgeo .* fn + tau0_max = peak .* tau0 + BoolEvent(free, fn) + Event(w, + Equation[ + reinit(startForward, + pre(mode) == Stuck & (sa > tau0_max/unitTorque | pre(startForward)) & + sa > tau0/unitTorque | pre(mode) == Backward & w > w_small | initial() & (w > 0)) + reinit(startBackward, + pre(mode) == Stuck & (sa < -tau0_max/unitTorque | pre(startBackward) & + sa < -tau0/unitTorque) | pre(mode) == Forward & w < -w_small | initial() & (w < 0)) + reinit(locked, + !free && !(pre(mode) == Forward | startForward | pre(mode) == Backward | startBackward)) + # finite state machine to determine configuration + reinit(mode, + ifelse(free, Free, + ifelse((pre(mode) == Forward | pre(mode) == Free | startForward) & w > 0, Forward, + ifelse((pre(mode) == Backward | pre(mode) == Free | startBackward) & w < 0, Backward, + Stuck)))) + ]) - # Friction torque - tau - ifelse(locked, - sa*unitTorque, - ifelse(free, - 0.0, - cgeo .* fn .* (ifelse(startForward, - tempInterpol1( w, mue_pos, 2), - ifelse(startBackward, - -tempInterpol1(-w, mue_pos, 2), - ifelse(pre(mode) == Forward, - tempInterpol1( w, mue_pos, 2), - -tempInterpol1(-w, mue_pos, 2))))))) + # Friction torque + tau = ifelse(locked, + sa*unitTorque, + ifelse(free, + 0.0, + cgeo .* fn .* (ifelse(startForward, + tempInterpol1( w, mue_pos, 2), + ifelse(startBackward, + -tempInterpol1(-w, mue_pos, 2), + ifelse(pre(mode) == Forward, + tempInterpol1( w, mue_pos, 2), + -tempInterpol1(-w, mue_pos, 2))))))) - a - unitAngularAcceleration .* - ifelse(locked, - 0.0, - ifelse(free, - sa, - ifelse(startForward, - sa - tau0_max ./ unitTorque, - ifelse(startBackward, - sa + tau0_max ./ unitTorque, - ifelse(pre(mode) == Forward, - sa - tau0_max ./ unitTorque, - sa + tau0_max ./ unitTorque))))) - } + a = unitAngularAcceleration .* + ifelse(locked, + 0.0, + ifelse(free, + sa, + ifelse(startForward, + sa - tau0_max ./ unitTorque, + ifelse(startBackward, + sa + tau0_max ./ unitTorque, + ifelse(pre(mode) == Forward, + sa - tau0_max ./ unitTorque, + sa + tau0_max ./ unitTorque))))) + end end function IdealGear(flange_a::Flange, flange_b::Flange, ratio::Real) val = compatible_values(flange_a, flange_b) tau_a = Torque(val) tau_b = Torque(val) - { - RefBranch(flange_a, tau_a) - RefBranch(flange_b, tau_b) - flange_a - ratio * flange_b - ratio * tau_a + tau_b - } + @equations begin + RefBranch(flange_a, tau_a) + RefBranch(flange_b, tau_b) + flange_a = ratio * flange_b + ratio * tau_a + tau_b + end end @@ -207,18 +207,18 @@ end ######################################## function SpeedSensor(flange::Flange, w::Signal) - { - w - der(flange) - } + @equations begin + w = der(flange) + end end function AccSensor(flange::Flange, a::Signal) w = AngularVelocity(compatible_values(flange)) - { - w - der(flange) - a - der(w) - } + @equations begin + w = der(flange) + a = der(w) + end end @@ -229,10 +229,10 @@ end function SignalTorque(flange_a::Flange, flange_b::Flange, tau::Signal) - { - RefBranch(flange_a, -tau) - RefBranch(flange_b, tau) - } + @equations begin + RefBranch(flange_a, -tau) + RefBranch(flange_b, tau) + end end function QuadraticSpeedDependentTorque(flange_a::Flange, flange_b::Flange, @@ -242,11 +242,11 @@ function QuadraticSpeedDependentTorque(flange_a::Flange, flange_b::Flange, tau = Torque(val) phi = Angle(val) w = AngularVelocity(val) - { - Branch(flange_b, flange_a, phi, tau) - w - der(phi) - tau - ifelse(TorqueDirection, - tau_nominal*(w/w_nominal)^2, - tau_nominal*ifelse(w >= 0, (w/w_nominal)^2, -(w/w_nominal)^2)) - } + @equations begin + Branch(flange_b, flange_a, phi, tau) + w = der(phi) + tau = ifelse(TorqueDirection, + tau_nominal*(w/w_nominal)^2, + tau_nominal*ifelse(w >= 0, (w/w_nominal)^2, -(w/w_nominal)^2)) + end end From e58e8f5e01882c1a22275fbf21ece6f9f1df953f Mon Sep 17 00:00:00 2001 From: Tom Short Date: Fri, 2 Jan 2015 15:02:38 -0500 Subject: [PATCH 2/2] Disable the test for sim_PID_Controller(). --- examples/run_all.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/run_all.jl b/examples/run_all.jl index 9096c5e..7db6483 100644 --- a/examples/run_all.jl +++ b/examples/run_all.jl @@ -8,7 +8,7 @@ module ex using Sims - +sim = sunsim path = Pkg.dir() * "/Sims/examples/" include(path * "basics/breaking_pendulum_in_box.jl") @@ -22,7 +22,7 @@ include(path * "basics/dc_motor_w_shaft.jl") include(path * "basics/half_wave_rectifiers.jl") include(path * "stdlib/m_blocks.jl") -sim_PID_Controller() +#sim_PID_Controller() include(path * "stdlib/m_electrical.jl") run_electrical_examples()