Merge pull request #313 from ven-k/vkb/tutorials-with-mtkmodel

Refactor all tutorials with `@mtkmodel` and `@mtkbuild`

Author: ChrisRackauckas

docs/src/tutorials/custom_component.md

@@ -114,14 +114,14 @@ Since the initial voltage of the capacitors was already specified via `v` and th
 
 ```@example components
 @mtkbuild sys = ChaoticAttractor()
-prob = ODEProblem(sys, Pair[], (0, 5e4), saveat = 0.01)
-sol = solve(prob)
+prob = ODEProblem(sys, Pair[], (0, 5e4))
+sol = solve(prob; saveat = 1.0)
 
-plot(sol[capacitor1.v], sol[capacitor2.v], title = "Chaotic Attractor", label = "",
+plot(sol[sys.capacitor1.v], sol[sys.capacitor2.v], title = "Chaotic Attractor", label = "",
     ylabel = "C1 Voltage in V", xlabel = "C2 Voltage in V")
 ```
 
 ```@example components
-plot(sol; idxs = [capacitor1.v, capacitor2.v, inductor.i],
+plot(sol; idxs = [sys.capacitor1.v, sys.capacitor2.v, sys.inductor.i],
     labels = ["C1 Voltage in V" "C2 Voltage in V" "Inductor Current in A"])
 ```

docs/src/tutorials/dc_motor_pi.md

@@ -20,51 +20,54 @@ using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq
 using Plots
-
-R = 0.5 # [Ohm] armature resistance
-L = 4.5e-3 # [H] armature inductance
-k = 0.5 # [N.m/A] motor constant
-J = 0.02 # [kg.m²] inertia
-f = 0.01 # [N.m.s/rad] friction factor
-tau_L_step = -0.3 # [N.m] amplitude of the load torque step
-nothing # hide
 ```
 
 The actual model can now be composed.
 
 ```@example dc_motor_pi
-systems = @named begin
-    ground = Ground()
-    source = Voltage()
-    ref = Blocks.Step(height = 1, start_time = 0)
-    pi_controller = Blocks.LimPI(k = 1.1, T = 0.035, u_max = 10, Ta = 0.035)
-    feedback = Blocks.Feedback()
-    R1 = Resistor(R = R)
-    L1 = Inductor(L = L)
-    emf = EMF(k = k)
-    fixed = Fixed()
-    load = Torque()
-    load_step = Blocks.Step(height = tau_L_step, start_time = 3)
-    inertia = Inertia(J = J)
-    friction = Damper(d = f)
-    speed_sensor = SpeedSensor()
+@mtkmodel DCMotor begin
+    @parameters begin
+        R = 0.5, [description = "Armature resistance"] # Ohm
+        L = 4.5e-3, [description = "Armature inductance"] # H
+        k = 0.5, [description = "Motor constant"] # N.m/A
+        J = 0.02, [description = "Inertia"] # kg.m²
+        f = 0.01, [description = "Friction factor"] # N.m.s/rad
+        tau_L_step = -0.3, [description = "Amplitude of the load torque step"] # N.m
+    end
+    @components begin
+        ground = Ground()
+        source = Voltage()
+        ref = Blocks.Step(height = 1, start_time = 0)
+        pi_controller = Blocks.LimPI(k = 1.1, T = 0.035, u_max = 10, Ta = 0.035)
+        feedback = Blocks.Feedback()
+        R1 = Resistor(R = R)
+        L1 = Inductor(L = L)
+        emf = EMF(k = k)
+        fixed = Fixed()
+        load = Torque()
+        load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+        inertia = Inertia(J = J)
+        friction = Damper(d = f)
+        speed_sensor = SpeedSensor()
+    end
+    @equations begin
+        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, :y, feedback.input2)
+        connect(feedback.output, pi_controller.err_input)
+        connect(pi_controller.ctr_output, :u, source.V)
+        connect(source.p, R1.p)
+        connect(R1.n, L1.p)
+        connect(L1.n, emf.p)
+        connect(emf.n, source.n, ground.g)
+    end
 end
 
-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, :y, feedback.input2)
-               connect(feedback.output, pi_controller.err_input)
-               connect(pi_controller.ctr_output, :u, 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)
+@named model = DCMotor()
 nothing # hide
 ```
 
@@ -75,12 +78,12 @@ so that it can be represented as a system of `ODEs` (ordinary differential equat
 ```@example dc_motor_pi
 sys = structural_simplify(model)
 prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0, 6.0))
-sol = solve(prob, Rodas4())
+sol = solve(prob)
 
-p1 = plot(sol.t, sol[inertia.w], ylabel = "Angular Vel. in rad/s",
+p1 = plot(sol.t, sol[sys.inertia.w], ylabel = "Angular Vel. in rad/s",
     label = "Measurement", title = "DC Motor with Speed Controller")
-plot!(sol.t, sol[ref.output.u], label = "Reference")
-p2 = plot(sol.t, sol[load.tau.u], ylabel = "Disturbance in Nm", label = "")
+plot!(sol.t, sol[sys.ref.output.u], label = "Reference")
+p2 = plot(sol.t, sol[sys.load.tau.u], ylabel = "Disturbance in Nm", label = "")
 plot(p1, p2, layout = (2, 1))
 ```
 
@@ -89,8 +92,8 @@ plot(p1, p2, layout = (2, 1))
 When implementing and tuning a control system in simulation, it is a good practice to analyze the closed-loop properties and verify robustness of the closed-loop with respect to, e.g., modeling errors. To facilitate this, we added two analysis points to the set of connections above, more specifically, we added the analysis points named `:y` and `:u` to the connections (for more details on analysis points, see [Linear Analysis](@ref))
 
 ```julia
-connect(speed_sensor.w, :y, feedback.input2)
-connect(pi_controller.ctr_output, :u, source.V)
+connect(sys.speed_sensor.w, :y, feedback.input2)
+connect(sys.pi_controller.ctr_output, :u, source.V)
 ```
 
 one at the plant output (`:y`) and one at the plant input (`:u`). We may use these analysis points to calculate, e.g., sensitivity functions, illustrated below. Here, we calculate the sensitivity function $S(s)$ and the complimentary sensitivity function $T(s) = I - S(s)$, defined as
@@ -108,7 +111,7 @@ matrices_S, simplified_sys = Blocks.get_sensitivity(
     model, :y, op = Dict(unknowns(sys) .=> 0.0))
 So = ss(matrices_S...) |> minreal # The output-sensitivity function as a StateSpace system
 matrices_T, simplified_sys = Blocks.get_comp_sensitivity(
-    model, :y, op = Dict(inertia.phi => 0.0, inertia.w => 0.0))
+    model, :y, op = Dict(model.inertia.phi => 0.0, model.inertia.w => 0.0))
 To = ss(matrices_T...)# The output complementary sensitivity function as a StateSpace system
 bodeplot([So, To], label = ["S" "T"], plot_title = "Sensitivity functions",
     plotphase = false)

docs/src/tutorials/rc_circuit.md

@@ -12,27 +12,32 @@ using ModelingToolkitStandardLibrary.Electrical
 using ModelingToolkitStandardLibrary.Blocks: Constant
 using ModelingToolkit: t_nounits as t
 
-R = 1.0
-C = 1.0
-V = 1.0
-systems = @named begin
-    resistor = Resistor(R = R)
-    capacitor = Capacitor(C = C, v = 0.0)
-    source = Voltage()
-    constant = Constant(k = V)
-    ground = Ground()
+@mtkmodel RC begin
+    @parameters begin
+        R = 1.0
+        C = 1.0
+        V = 1.0
+    end
+    @components begin
+        resistor = Resistor(R = R)
+        capacitor = Capacitor(C = C, v = 0.0)
+        source = Voltage()
+        constant = Constant(k = V)
+        ground = Ground()
+    end
+    @equations begin
+        connect(constant.output, source.V)
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n, ground.g)
+    end
 end
 
-rc_eqs = [connect(constant.output, source.V)
-          connect(source.p, resistor.p)
-          connect(resistor.n, capacitor.p)
-          connect(capacitor.n, source.n, ground.g)]
-
-@named rc_model = ODESystem(rc_eqs, t; systems)
-sys = structural_simplify(rc_model)
+@mtkbuild sys = RC()
 prob = ODEProblem(sys, Pair[], (0, 10.0))
-sol = solve(prob, Tsit5())
-plot(sol, idxs = [capacitor.v, resistor.i],
+sol = solve(prob)
+
+plot(sol, idxs = [sys.capacitor.v, sys.resistor.i],
     title = "RC Circuit Demonstration",
     labels = ["Capacitor Voltage" "Resistor Current"])
 ```

docs/src/tutorials/thermal_model.md

@@ -10,31 +10,34 @@ from dividing the total initial energy in the system by the sum of the heat capa
 using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Plots
 using ModelingToolkit: t_nounits as t
 
-C1 = 15
-C2 = 15
-systems = @named begin
-    mass1 = HeatCapacitor(C = C1, T = 373.15)
-    mass2 = HeatCapacitor(C = C2, T = 273.15)
-    conduction = ThermalConductor(G = 10)
-    Tsensor1 = TemperatureSensor()
-    Tsensor2 = TemperatureSensor()
+@mtkmodel HeatConductionModel begin
+    @parameters begin
+        C1 = 15
+        C2 = 15
+    end
+    @components begin
+        mass1 = HeatCapacitor(C = C1, T = 373.15)
+        mass2 = HeatCapacitor(C = C2, T = 273.15)
+        conduction = ThermalConductor(G = 10)
+        Tsensor1 = TemperatureSensor()
+        Tsensor2 = TemperatureSensor()
+    end
+    @equations begin
+        connect(mass1.port, conduction.port_a)
+        connect(conduction.port_b, mass2.port)
+        connect(mass1.port, Tsensor1.port)
+        connect(mass2.port, Tsensor2.port)
+    end
 end
 
-connections = [
-    connect(mass1.port, conduction.port_a),
-    connect(conduction.port_b, mass2.port),
-    connect(mass1.port, Tsensor1.port),
-    connect(mass2.port, Tsensor2.port)
-]
-
-@named model = ODESystem(connections, t; systems)
-sys = structural_simplify(model)
+@mtkbuild sys = HeatConductionModel()
 prob = ODEProblem(sys, Pair[], (0, 5.0))
-sol = solve(prob, Tsit5())
+sol = solve(prob)
 
-T_final_K = sol[(mass1.T * C1 + mass2.T * C2) / (C1 + C2)]
+T_final_K = sol[(sys.mass1.T * sys.C1 + sys.mass2.T * sys.C2) / (sys.C1 + sys.C2)]
 
 plot(title = "Thermal Conduction Demonstration")
-plot!(sol, idxs = [mass1.T, mass2.T], labels = ["Mass 1 Temperature" "Mass 2 Temperature"])
+plot!(sol, idxs = [sys.mass1.T, sys.mass2.T],
+    labels = ["Mass 1 Temperature" "Mass 2 Temperature"])
 plot!(sol.t, T_final_K, label = "Steady-State Temperature")
 ```