Skip to content

Commit

Permalink
Merge pull request #313 from ven-k/vkb/tutorials-with-mtkmodel
Browse files Browse the repository at this point in the history
Refactor all tutorials with `@mtkmodel` and `@mtkbuild`
  • Loading branch information
ChrisRackauckas authored Jul 29, 2024
2 parents 7e3f5b8 + 802c415 commit 31114a0
Show file tree
Hide file tree
Showing 4 changed files with 98 additions and 87 deletions.
8 changes: 4 additions & 4 deletions docs/src/tutorials/custom_component.md
Original file line number Diff line number Diff line change
Expand Up @@ -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"])
```
93 changes: 48 additions & 45 deletions docs/src/tutorials/dc_motor_pi.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
```

Expand All @@ -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))
```

Expand All @@ -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
Expand All @@ -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)
Expand Down
41 changes: 23 additions & 18 deletions docs/src/tutorials/rc_circuit.md
Original file line number Diff line number Diff line change
Expand Up @@ -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"])
```
43 changes: 23 additions & 20 deletions docs/src/tutorials/thermal_model.md
Original file line number Diff line number Diff line change
Expand Up @@ -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")
```

0 comments on commit 31114a0

Please sign in to comment.