Merge branch 'bgc/multi2d' into bgc/revolute

Author: Brad Carman

docs/pages.jl

@@ -4,6 +4,7 @@ pages = [
         "RC Circuit" => "tutorials/rc_circuit.md",
         "Custom Components" => "tutorials/custom_component.md",
         "Thermal Conduction Model" => "tutorials/thermal_model.md",
+        "DC Motor with Speed Controller" => "tutorials/dc_motor_pi.md",
     ],
     "About Acausal Connections" => "connectors/connections.md",
     "API" => [

docs/src/API/linear_analysis.md

@@ -79,6 +79,7 @@ L = P*(-C) # Add the minus sign to build the negative feedback into the controll
 
 To obtain the transfer function between two analysis points, we call `linearize`
 ```@example LINEAR_ANALYSIS
+using ModelingToolkit # hide
 matrices_PS = linearize(sys, :plant_input, :plant_output)[1]
 ```
 this particular transfer function should be equivalent to the linear system `P(s)S(s)`, i.e., equivalent to

docs/src/index.md

@@ -17,6 +17,7 @@ import Pkg; Pkg.add("ModelingToolkitStandardLibrary")
 - [RC Circuit](http://mtkstdlib.sciml.ai/dev/tutorials/rc_circuit/)
 - [Custom Component](http://mtkstdlib.sciml.ai/dev/tutorials/custom_component/)
 - [Thermal Model](http://mtkstdlib.sciml.ai/dev/tutorials/thermal_model/)
+- [DC Motor with PI-controller](http://mtkstdlib.sciml.ai/dev/tutorials/dc_motor_pi/)
 
 ## Libraries
 

docs/src/tutorials/dc_motor_pi.md

@@ -0,0 +1,93 @@
+# DC Motor with PI-controller
+In this example a PI-controller is setup for speed control of a DC-motor. First the needed packages 
+are imported and the parameters of the model defined.
+
+```@example dc_motor_pi
+using ModelingToolkit
+using ModelingToolkitStandardLibrary.Electrical
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using OrdinaryDiffEq
+using Plots
+
+@parameters t
+
+R = 0.5 # [Ohm]
+L = 4.5e-3 # [H]
+k = 0.5 # [N.m/A]
+J = 0.02 # [kg.m²]
+f = 0.01 # [N.m.s/rad]
+tau_L_step = -3 # [N.m]
+nothing # hide
+```
+
+The actual model can now be composed.
+
+```@example dc_motor_pi
+@named ground = Ground()
+@named source = Voltage()
+@named ref = Blocks.Step(height = 1, start_time = 0)
+@named pi_controller = Blocks.LimPI(k = 1.1, T = 0.035, u_max = 0.6, Ta = 0.035)
+@named feedback = Blocks.Feedback()
+@named R1 = Resistor(R = R)
+@named L1 = Inductor(L = L)
+@named emf = EMF(k = k)
+@named fixed = Fixed()
+@named load = Torque(use_support = false)
+@named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+@named inertia = Inertia(J = J)
+@named friction = Damper(d = f)
+@named speed_sensor = SpeedSensor()
+
+connections = [connect(fixed.flange, emf.support, friction.flange_b)
+               connect(emf.flange, friction.flange_a, inertia.flange_a)
+               connect(inertia.flange_b, load.flange)
+               connect(inertia.flange_b, speed_sensor.flange)
+               connect(load_step.output, load.tau)
+               connect(ref.output, feedback.input1)
+               connect(speed_sensor.w, feedback.input2)
+               connect(feedback.output, pi_controller.err_input)
+               connect(pi_controller.ctr_output, source.V)
+               connect(source.p, R1.p)
+               connect(R1.n, L1.p)
+               connect(L1.n, emf.p)
+               connect(emf.n, source.n, ground.g)]
+
+@named model = ODESystem(connections, t,
+                         systems = [
+                             ground,
+                             ref,
+                             pi_controller,
+                             feedback,
+                             source,
+                             R1,
+                             L1,
+                             emf,
+                             fixed,
+                             load,
+                             load_step,
+                             inertia,
+                             friction,
+                             speed_sensor,
+                         ])
+nothing # hide
+```
+
+Now the model can be simulated. Typical rotational mechanical systems are described via `DAE` 
+(differential algebraic equations), however in this case ModelingToolkit can simplify the model enough
+so that it can be represented as a system of `ODEs` (ordinary differential equations).
+
+```@example dc_motor_pi
+sys = structural_simplify(model)
+prob = ODEProblem(sys, [], (0, 6.0))
+sol = solve(prob, Rodas4())
+
+p1 = Plots.plot(sol.t, sol[inertia.w], ylabel = "Angular Vel. in rad/s",
+                label = "Measurement", title = "DC Motor with Speed Controller")
+Plots.plot!(sol.t, sol[ref.output.u], label = "Reference")
+p2 = Plots.plot(sol.t, sol[load.tau.u], ylabel = "Disturbance in Nm", label = "")
+Plots.plot(p1, p2, layout = (2, 1))
+nothing # hide
+```
+
+![DC Motor with speed controller](dc_motor_pi.png)

links.mp4

@@ -26,7 +26,7 @@ See [OnePort](@ref)
 - `n` Negative pin
 
 # Parameters:
-- `R`: [`Ω`] Resistance
+- `R`: [`Ohm`] Resistance
 """
 function Resistor(; name, R)
     @named oneport = OnePort()
@@ -70,8 +70,7 @@ end
 Creates an ideal capacitor.
 
 # States:
-- `v(t)`: [`V`]
-  The voltage across the capacitor, given by `D(v) ~ p.i / C`
+- `v(t)`: [`V`] The voltage across the capacitor, given by `D(v) ~ p.i / C`
 
 # Connectors:
 - `p` Positive pin
@@ -160,3 +159,73 @@ function Short(; name)
     eqs = [v ~ 0]
     extend(ODESystem(eqs, t, [], []; name = name), oneport)
 end
+
+"""
+    HeatingResistor(;name, R_ref=1.0, T_ref=300.15, alpha=0)
+
+Temperature dependent electrical resistor
+
+# States
+- See [OnePort](@ref)
+- `R(t)`: [`Ohm`] Temperature dependent resistance `R ~ R_ref*(1 + alpha*(heat_port.T(t) - T_ref))`
+
+# Connectors
+- `p` Positive pin
+- `n` Negative pin
+
+# Parameters: 
+- `R_ref`: [`Ω`] Reference resistance 
+- `T_ref`: [K] Reference temperature
+"""
+function HeatingResistor(; name, R_ref = 1.0, T_ref = 300.15, alpha = 0)
+    @named oneport = OnePort()
+    @unpack v, i = oneport
+    @named heat_port = HeatPort()
+    pars = @parameters begin
+        R_ref = R_ref
+        T_ref = T_ref
+        alpha = alpha
+    end
+    @variables R(t) = R_ref
+    eqs = [R ~ R_ref * (1 + alpha * (heat_port.T - T_ref))
+           heat_port.Q_flow ~ -v * i # -LossPower
+           v ~ i * R]
+    extend(ODESystem(eqs, t, [R], pars; name = name, systems = [heat_port]), oneport)
+end
+
+"""
+    EMF(;name, k)
+
+Electromotoric force (electric/mechanic transformer)
+
+# States
+- `v(t)`: [`V`] The voltage across component `p.v - n.v`
+- `i(t)`: [`A`] The current passing through positive pin
+- `phi`: [`rad`] Rotation angle (=flange.phi - support.phi)
+- `w`: [`rad/s`] Angular velocity (= der(phi))
+
+# Connectors
+- `p` [Pin](@ref) Positive pin
+- `n` [Pin](@ref) Negative pin
+- `flange` [Flange](@ref) Shaft of EMF shaft
+- `support` [Support](@ref) Support/housing of emf shaft
+
+# Parameters: 
+- `k`: [`N⋅m/A`] Transformation coefficient 
+"""
+function EMF(; name, k)
+    @named p = Pin()
+    @named n = Pin()
+    @named flange = Flange()
+    @named support = Support()
+    @parameters k = k
+    @variables v(t)=0.0 i(t)=0.0 phi(t)=0.0 w(t)=0.0
+    eqs = [v ~ p.v - n.v
+           0 ~ p.i + n.i
+           i ~ p.i
+           phi ~ flange.phi - support.phi
+           D(phi) ~ w
+           k * w ~ v
+           flange.tau ~ -k * i]
+    ODESystem(eqs, t, [v, i, phi, w], [k]; name = name, systems = [p, n, flange, support])
+end

src/Electrical/Analog/ideal_components.jl

@@ -6,16 +6,26 @@ module Electrical
 
 using ModelingToolkit, Symbolics, IfElse
 using OffsetArrays
+using ..Thermal: HeatPort
+using ..Mechanical.Rotational: Flange, Support
+using ..Blocks: RealInput, RealOutput
 
 @parameters t
 D = Differential(t)
 
-using ..Blocks: RealInput, RealOutput
-
+export Pin, OnePort
 include("utils.jl")
+
+export Capacitor, Ground, Inductor, Resistor, Conductor, Short, IdealOpAmp, EMF,
+       HeatingResistor
 include("Analog/ideal_components.jl")
+
+export CurrentSensor, PotentialSensor, VoltageSensor, PowerSensor, MultiSensor
 include("Analog/sensors.jl")
+
+export Voltage, Current
 include("Analog/sources.jl")
+
 # include("Digital/components.jl")
 # include("Digital/gates.jl")
 # include("Digital/tables.jl")
@@ -26,22 +36,4 @@ include("Analog/sources.jl")
 # - machines
 # - multi-phase
 
-export #Interface
-       Pin,
-# Analog Components
-       Capacitor, Ground, Inductor, Resistor, Conductor,
-       Short, IdealOpAmp,
-# Analog Sensors
-       CurrentSensor, PotentialSensor, VoltageSensor,
-       PowerSensor, MultiSensor,
-# Analog Sources
-       Voltage, Current
-
-# # Digital Gates
-# And, Or, Not, Xor, Nand, Nor, Xnor,
-# # Digital components
-# HalfAdder, FullAdder, MUX, DEMUX, Encoder, Decoder,
-# # Digital Sources
-# DigitalPin, Pulse, PulseDiff
-
 end

src/Electrical/Electrical.jl

@@ -39,6 +39,7 @@ end
 function Inertia(; name, J, phi_start = 0.0, w_start = 0.0, a_start = 0.0)
     @named flange_a = Flange()
     @named flange_b = Flange()
+    J > 0 || throw(ArgumentError("Expected `J` to be positive"))
     @parameters J = J
     sts = @variables begin
         phi(t) = phi_start
@@ -73,6 +74,7 @@ Linear 1D rotational spring
 function Spring(; name, c, phi_rel0 = 0.0)
     @named partial_comp = PartialCompliant()
     @unpack phi_rel, tau = partial_comp
+    c > 0 || throw(ArgumentError("Expected `c` to be positive"))
     pars = @parameters begin
         c = c
         phi_rel0 = phi_rel0
@@ -102,6 +104,7 @@ Linear 1D rotational damper
 function Damper(; name, d)
     @named partial_comp = PartialCompliantWithRelativeStates()
     @unpack w_rel, tau = partial_comp
+    d > 0 || throw(ArgumentError("Expected `d` to be positive"))
     pars = @parameters d = d
     eqs = [tau ~ d * w_rel]
     extend(ODESystem(eqs, t, [], pars; name = name), partial_comp)
@@ -130,6 +133,7 @@ This element characterizes any type of gear box which is fixed in the ground and
 function IdealGear(; name, ratio, use_support = false)
     @named partial_element = PartialElementaryTwoFlangesAndSupport2(use_support = use_support)
     @unpack phi_support, flange_a, flange_b = partial_element
+    ratio > 0 || throw(ArgumentError("Expected `ratio` to be positive"))
     @parameters ratio = ratio
     sts = @variables phi_a(t)=0.0 phi_b(t)=0.0
     eqs = [phi_a ~ flange_a.phi - phi_support

src/Mechanical/Rotational/components.jl

@@ -1,5 +1,5 @@
 """
-    Torque(;name) 
+    Torque(; name, use_support=false) 
 
 Input signal acting as external torque on a flange
 

src/Mechanical/Rotational/sources.jl

@@ -118,7 +118,7 @@ Partial model for a component with one rotational 1-dim. shaft flange and a supp
 function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
     @named flange = Flange()
     sys = [flange]
-    @variables phi_support(t)
+    @variables phi_support(t) = 0.0
     if use_support
         @named support = Support()
         eqs = [support.phi ~ phi_support

src/Mechanical/Rotational/utils.jl

@@ -2,8 +2,8 @@ module ModelingToolkitStandardLibrary
 
 include("Blocks/Blocks.jl")
 include("Mechanical/Mechanical.jl")
+include("Thermal/Thermal.jl")
 include("Electrical/Electrical.jl")
 include("Magnetic/Magnetic.jl")
-include("Thermal/Thermal.jl")
 
 end

src/ModelingToolkitStandardLibrary.jl

@@ -102,7 +102,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u]≈cosine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3
 
@@ -169,7 +168,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u]≈ramp.(sol.t, offset, height, duration, start_time) atol=1e-3
 
@@ -365,7 +363,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u]≈exp_sine.(sol.t, amplitude, frequency, damping, phase,
                                       start_time) atol=1e-3

test/Blocks/sources.jl

@@ -157,10 +157,10 @@ end
                                  systems = [resistor, capacitor, source, ground, voltage])
         sys = structural_simplify(model)
         prob = ODAEProblem(sys, [capacitor.v => 10.0], (0.0, 10.0))
-        sol = solve(prob, Rodas5())
-        @test sol.retcode == :Success
         sol = solve(prob, Tsit5())
         @test sol.retcode == :Success
+        sol = solve(prob, Rodas4())
+        @test sol.retcode == :Success
 
         # Plots.plot(sol; vars=[voltage.v, capacitor.v])
     end