ModelingToolkit.jl from the ground up: a DC motor under speed control

The problem: a DC motor that must hold 100 rad/s

A small permanent-magnet DC motor has to hold its load at 100 rad/s, about 955 rpm. The winding has 1 Ω and 1 mH, the motor constant is 0.05 V·s/rad, and rotor and load have an inertia of \(10^{-4}\) kg·m² with a little viscous friction. A PI controller sets the voltage. The setpoint jumps from 0 to 100 rad/s at the start, and at 0.2 s a load torque of 0.03 N·m switches on. ModelingToolkit.jl builds this model from the parts of the circuit and derives its equations itself.

Written by hand, the model is three equations in three states, the current \(i\), the speed \(\omega\), and the controller's integrated error:

\[L\frac{di}{dt} = u - R\,i - k\,\omega, \qquad J\frac{d\omega}{dt} = k\,i - d\,\omega - \tau_\text{load}, \qquad u = k_p\Big(e + \frac{1}{T}\int e\,dt\Big).\]

Here \(e = \omega_\text{ref} - \omega\) is the speed error, and the controller's voltage \(u\) is proportional to the error plus its integral, with gain \(k_p\) and integral time \(T\). Add a gear or a voltage limit, and you rewrite these by hand, signs included.

ModelingToolkit takes the other route. Each component carries its own equations, the resistor Ohm's law and the inertia Newton's law, and every connection between ports adds Kirchhoff's laws or the balance of torques. You draw the circuit in code, and the library collects and simplifies the equations.

Speed in rad/s (top) and armature current in A (bottom) of a DC motor over 0.4 s for two PI gains. With kp = 0.1 the speed settles at 100 rad/s in 78 ms at a 9.1 A current peak; with kp = 0.02 it is still 14 rad/s short when the load switches on at 0.2 s.

The low gain has not reached the setpoint when the load arrives. Five times the gain settles in 78 ms and pays with a current peak of 9.1 A. Both curves come from one compiled model, solved twice.

Setup

The model's components come from ModelingToolkitStandardLibrary, imported by name so that you see which domain each one belongs to. As in DifferentialEquations.jl from the ground up, the solvers come from OrdinaryDiffEq, the part of the DifferentialEquations.jl umbrella that holds them; using DifferentialEquations runs every line below unchanged. Install once with import Pkg; Pkg.add(["ModelingToolkit", "ModelingToolkitStandardLibrary", "OrdinaryDiffEq", "CairoMakie"]). The code assumes ModelingToolkit 11; older versions use other names for the same functions, listed in the pitfalls. Step 1 explains t and D. Be patient with the first run: the whole notebook took 223 s on the machine that built this page, most of it Julia compiling ModelingToolkit's functions for this model. A retuned solve in the same session takes seconds.

using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq, Printf, CairoMakie
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary.Electrical: OnePort, Voltage, Resistor, Inductor, Ground
using ModelingToolkitStandardLibrary.Mechanical.Rotational: Flange, Inertia, Damper, Fixed, Torque, SpeedSensor
using ModelingToolkitStandardLibrary.Blocks: Step, Feedback, PI

const R, L, k = 1.0, 1e-3, 0.05            # Ω, H, V·s/rad
const J, d = 1e-4, 1e-5                    # kg·m², N·m·s/rad
const ω_ref, τ_load, t_load = 100.0, 0.03, 0.2   # rad/s, N·m, s

const INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"
set_theme!(Theme(                          # the look of every figure below
    size = (770, 396), fontsize = 17,
    palette = (color = [INK, ACCENT, SECOND, MUTED],),
    Axis = (topspinevisible = false, rightspinevisible = false, xgridvisible = true, ygridvisible = true),
    Lines = (linewidth = 2.5,),
))

println("ModelingToolkit ", pkgversion(ModelingToolkit),
        ", StandardLibrary ", pkgversion(ModelingToolkitStandardLibrary))

ModelingToolkit 11.45.3, StandardLibrary 2.29.8

Step 1: Write one component, the EMF that couples current and torque

The one component you write yourself is the motor's coupling element: it turns current into torque and speed into back-EMF. Step 2 plugs it into the motor next to the library's parts. Start by reading a library component, the base of every two-terminal element:

@named oneport = OnePort()
equations(oneport)

3-element Vector{Equation}:
 v(t) ~ -n₊v(t) + p₊v(t)
 [34m0[39m ~ n₊i(t) + p₊i(t)
 i(t) ~ p₊i(t)

Three equations: the voltage \(v\) across the element is the difference of the potentials at its pins p and n, the currents into the two pins sum to zero, and \(i\) is the current entering at p. Every component is a function with a name keyword, and @named oneport = OnePort() is shorthand for oneport = OnePort(; name = :oneport). The name becomes the prefix of every variable inside, which is why the pins print as p₊v and n₊i.

Now the EMF:

@component function MotorEMF(; name, k)
    @named oneport = OnePort()
    @unpack v, i = oneport
    @named flange = Flange()
    @variables phi(t) w(t)
    @parameters k = k
    eqs = [phi ~ flange.phi,            # the shaft angle at the flange
           D(phi) ~ w,
           v ~ k * w,                   # back-EMF
           flange.tau ~ -k * i]         # torque into the EMF; the shaft gets +k*i
    return extend(System(eqs, t, [phi, w], [k]; name, systems = [flange]), oneport)
end

MotorEMF (generic function with 1 method)

t is the symbolic time and D the operator \(d/dt\), so D(phi) ~ w reads \(d\varphi/dt = \omega\). The ~ writes an equation, not an assignment. @variables declares quantities that change in time, @parameters constants that you can retune later without rebuilding. @unpack binds the v and i of the subsystem oneport to local names.

A flange is a connector, like a pin: it carries a potential, the angle phi, and a flow, the torque tau. The back-EMF needs the speed and the flange only gives the angle, hence phi and its derivative w. A flange torque counts positive into the component: the shaft pushes on the EMF with \(-k\,i\), so the EMF drives the shaft with \(+k\,i\). extend merges the pins and the three equations of oneport into the new system, and @component makes @named treat the function like a library component:

@named emf = MotorEMF(k = k)
equations(emf)

7-element Vector{Equation}:
 v(t) ~ -n₊v(t) + p₊v(t)
 [34m0[39m ~ n₊i(t) + p₊i(t)
 i(t) ~ p₊i(t)
 phi(t) ~ flange₊phi(t)
 Differential(t, 1)(phi(t)) ~ w(t)
 v(t) ~ k*w(t)
 flange₊tau(t) ~ -i(t)*k

Seven equations, three inherited and four written here. The library has its own EMF with a second flange for the housing; use it when the housing moves.

Step 2: Connect the motor and its controller from library components

Here is the model you are about to write, each box named as in the code:

Show code

fig = Figure(size = (770, 484))
ax = Axis(fig[1, 1], limits = (-0.6, 19.0, -0.6, 9.2))
hidedecorations!(ax); hidespines!(ax)
boxes = Dict{String, NTuple{4, Float64}}()              # name => (x, y, half width, half height)
function block!(name, x, y; color = INK)
    hw, hh = 0.1 * length(name) + 0.3, 0.45
    poly!(ax, Rect2f(x - hw, y - hh, 2hw, 2hh), color = :white, strokecolor = color, strokewidth = 2)
    text!(ax, x, y, text = name, align = (:center, :center), fontsize = 16, color = INK)
    boxes[name] = (x, y, hw, hh)
end
function side(name, s)
    x, y, hw, hh = boxes[name]
    s == :l ? Point2f(x - hw, y) : s == :r ? Point2f(x + hw, y) : s == :t ? Point2f(x, y + hh) : Point2f(x, y - hh)
end
wire!(pts...) = lines!(ax, collect(Point2f.(pts)), color = INK, linewidth = 1.8)
function signal!(pts...)
    p = collect(Point2f.(pts))
    lines!(ax, p, color = SECOND, linewidth = 1.6)
    u = (p[end] - p[end-1]) / sqrt(sum(abs2, p[end] - p[end-1]))
    arrows2d!(ax, [p[end] - 0.35u], [0.35u], color = SECOND, shaftwidth = 0, tipwidth = 10, tiplength = 10)
end
port!(x, y, s; align = (:left, :bottom)) = text!(ax, x, y, text = s, fontsize = 14, color = MUTED, align = align)

block!("source", 1.2, 5.0); block!("resistor", 3.8, 7.6); block!("inductor", 6.8, 7.6)
block!("emf", 9.0, 5.0; color = ACCENT); block!("ground", 5.0, 2.7)
wire!(side("source", :t), (1.2, 7.6), side("resistor", :l))
wire!(side("resistor", :r), side("inductor", :l))
wire!(side("inductor", :r), (9.0, 7.6), side("emf", :t))
wire!(side("emf", :b), (9.0, 3.6), (1.2, 3.6), side("source", :b))
wire!((5.0, 3.6), side("ground", :t))
port!(1.35, 5.55, "p"); port!(1.35, 4.45, "n"; align = (:left, :top))
port!(9.15, 5.55, "p"); port!(9.15, 4.45, "n"; align = (:left, :top))

block!("inertia", 12.0, 5.0); block!("damper", 15.0, 7.0); block!("fixed", 17.6, 7.0)
block!("load", 15.0, 5.0); block!("load_step", 17.6, 5.0); block!("speed_sensor", 15.7, 3.0)
wire!(side("emf", :r), side("inertia", :l))
wire!(side("inertia", :r), side("load", :l))
wire!((13.5, 5.0), (13.5, 7.0), side("damper", :l))
wire!((13.5, 5.0), (13.5, 3.0), side("speed_sensor", :l))
wire!(side("damper", :r), side("fixed", :l))
scatter!(ax, [13.5], [5.0], color = INK, markersize = 7)
port!(9.72, 5.12, "flange")

block!("controller", 4.0, 1.0); block!("feedback", 8.4, 1.0); block!("setpoint", 12.4, 1.0)
signal!(side("setpoint", :l), side("feedback", :r))
signal!(side("speed_sensor", :b), (15.7, 0.0), (8.4, 0.0), side("feedback", :b))
signal!(side("feedback", :l), side("controller", :r))
signal!(side("controller", :l), (-0.2, 1.0), (-0.2, 5.0), side("source", :l))
signal!(side("load_step", :l), side("load", :r))
port!(15.85, 2.0, "w"; align = (:left, :center)); port!(9.65, 1.12, "input1")
port!(8.25, 0.15, "input2"; align = (:right, :bottom)); port!(-0.1, 5.15, "V")
text!(ax, 4.0, 8.6, text = "electrical", color = MUTED, fontsize = 16)
text!(ax, 14.3, 8.6, text = "rotational", color = MUTED, fontsize = 16)
text!(ax, 0.6, -0.5, text = "signal", color = MUTED, fontsize = 16)
fig

Diagram of the motor model. Electrical loop: source, resistor, inductor, emf, and ground. Shaft: emf drives inertia, linked to a damper to the fixed housing, a load torque, and a speed sensor. Signal loop, with arrows: speed sensor and setpoint into feedback, then controller, then the source.

The electrical loop runs on the left, the shaft on the right, and emf joins them. Below, the signal loop closes from the speed sensor through feedback, which subtracts the measured speed from the setpoint, and the controller back to the source's voltage input V. Physical connections have no arrowheads because they have no direction: current and torque flow whichever way the equations say. Signals do.

Each box becomes one @named line, and each line of the diagram an argument of connect:

@named source = Voltage()
@named resistor = Resistor(R = R)
@named inductor = Inductor(L = L)
@named ground = Ground()
@named inertia = Inertia(J = J)
@named damper = Damper(d = d)
@named fixed = Fixed()
@named load = Torque()
@named speed_sensor = SpeedSensor()
@named setpoint = Step(height = ω_ref, start_time = 0.0)
@named load_step = Step(height = -τ_load, start_time = t_load)   # a torque against the motion
@named feedback = Feedback()
@named controller = PI(k = 0.02, T = 0.04)

connections = [
    # electrical
    connect(source.p, resistor.p),
    connect(resistor.n, inductor.p),
    connect(inductor.n, emf.p),
    connect(emf.n, source.n, ground.g),
    # rotational
    connect(emf.flange, inertia.flange_a),
    connect(inertia.flange_b, damper.flange_a, load.flange, speed_sensor.flange),
    connect(damper.flange_b, fixed.flange),
    # signal
    connect(setpoint.output, feedback.input1),
    connect(speed_sensor.w, feedback.input2),
    connect(feedback.output, controller.err_input),
    connect(controller.ctr_output, source.V),
    connect(load_step.output, load.tau),
]
components = [source, resistor, inductor, emf, ground, inertia, damper, fixed,
              load, speed_sensor, setpoint, load_step, feedback, controller]
@named model = System(connections, t; systems = components)

[0m[1mModel model:[22m
[0m[1mSubsystems (14):[22m see hierarchy(model)
  source
  resistor
  inductor
  emf
  ground
  inertia
  ⋮
[0m[1mEquations (91):[22m
  64 standard: see equations(model)
  27 connecting: see equations(expand_connections(model))
[0m[1mUnknowns (74):[22m see unknowns(model)
  source₊v(t)
  source₊i(t)
  source₊p₊v(t)
  source₊p₊i(t)
  source₊n₊v(t)
  source₊n₊i(t)
[0m  ⋮
[0m[1mParameters (26):[22m see parameters(model)
  resistor₊R: Reference resistance
  resistor₊T_ref: Reference temperature
  resistor₊alpha: Temperature coefficient of resistance
  inductor₊L: Inductance
  emf₊k
  inertia₊J: Moment of inertia
[0m  ⋮

In PI(k = 0.02, T = 0.04) the k is the controller's gain \(k_p\), not the motor constant, and T the integral time; Step 4 shows why 0.04 s. ?PI lists a component's parameters and port names, such as err_input and ctr_output.

On physical ports, connect sets the potentials equal and makes the flows sum to zero. On signal ports it copies the output into the input. equations(model) would still show the connect statements; expand_connections turns them into plain equations:

flat = equations(expand_connections(model))
println(length(flat), " equations")
foreach(println, filter(eq -> occursin("ground", string(eq)), flat))

74 equations
ground₊g₊v(t) ~ 0
emf₊n₊v(t) ~ ground₊g₊v(t)
0 ~ emf₊n₊i(t) + ground₊g₊i(t) + source₊n₊i(t)

Fourteen components and twelve connect lines make 74 equations. The three that mention the ground show what a connection becomes. The ground fixes its own potential at zero. The node where emf.n, source.n, and ground.g meet gets one voltage for all three pins; the line shown is one of the two equalities. The currents into the node sum to zero, which is Kirchhoff's current law.

Step 3: Let mtkcompile reduce the model and read its equations

mtkcompile removes every variable that another equation already fixes and keeps as few unknowns as it can. unknowns lists what the solver integrates, observed the equations for the rest. equations(sys) would write the derivatives through those observed variables, inductor₊v / inductor₊L; full_equations substitutes them, so each right-hand side holds only parameters and unknowns:

sys = mtkcompile(model)
println(unknowns(sys))
println(length(observed(sys)), " observed variables")
foreach(println, full_equations(sys))

SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}[emf₊i(t), emf₊flange₊phi(t), controller₊addPI₊input2₊u(t), inertia₊w(t)]
76 observed variables
Differential(t, 1)(emf₊i(t)) ~ (-emf₊i(t)*resistor₊R - emf₊k*inertia₊w(t) + (controller₊addPI₊input2₊u(t)*controller₊addPI₊k2 + (-inertia₊w(t) + setpoint₊offset + (0.5 + atan((-setpoint₊start_time + t) / 1.0e-5) / π)*setpoint₊height)*controller₊addPI₊k1)*controller₊k) / inductor₊L
Differential(t, 1)(emf₊flange₊phi(t)) ~ inertia₊w(t)
Differential(t, 1)(controller₊addPI₊input2₊u(t)) ~ (-inertia₊w(t) + setpoint₊offset + (0.5 + atan((-setpoint₊start_time + t) / 1.0e-5) / π)*setpoint₊height) / controller₊T
Differential(t, 1)(inertia₊w(t)) ~ (load_step₊offset + (0.5 + atan((-load_step₊start_time + t) / 1.0e-5) / π)*load_step₊height - damper₊d*inertia₊w(t) + emf₊i(t)*emf₊k) / inertia₊J

Four unknowns are left of the 74. emf₊i is the armature current: resistor, inductor, and EMF carry the same current, and the simplification kept one name for it. inertia₊w is the speed, and controller₊addPI₊input2₊u is the integrator's output, the integrated error divided by \(T\), named after the wire it feeds. Inside, PI is three blocks: int integrates the error with gain \(1/T\), addPI adds error and integral with weights k1 and k2, both 1, and gainPI multiplies by k; setpoint₊offset is zero. With these values the first and the last equation are the armature and shaft equations from the top, signs included, with the load step entering as a negative torque. The setpoint is a step smoothed by an arctangent over about 10 μs, so the solver never sees a jump.

The fourth unknown, emf₊flange₊phi, is the shaft angle. Rotational ports carry the angle as their potential, and MotorEMF and Inertia both define the speed as its derivative, so the angle is integrated along. Nothing in this model depends on it. The other 70 variables became observed equations, plus six angle derivatives such as damper₊phi_relˍt created on the way: 76 in all, evaluated from the four unknowns only when you ask for them.

Step 4: Simulate a speed step and a load step

An ODEProblem takes the compiled system, one map of start values, and the time span. Every unknown without a default needs a start value, under any of its alias names: sys.inductor.i is emf₊i, sys.inertia.phi is emf₊flange₊phi. The angle needs one too, or initialization stops as incomplete; the PI integrator defaults to zero. The tolerances are the ones the prerequisite recommends for anything you plot:

prob = ODEProblem(sys, [sys.inductor.i => 0.0, sys.inertia.w => 0.0, sys.inertia.phi => 0.0], (0.0, 0.4))
sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-10)
println(sol.retcode, ", ", length(sol.t), " steps")

Success, 284 steps

You index the solution by symbol, not by position. sol[sys.inertia.w] is the speed at every step, and sol[sys.resistor.i] the current through the resistor, which is not an unknown but comes back anyway: the solution evaluates the observed equations on demand. sol(ts; idxs = sys.inertia.w) interpolates one variable at the times ts. metrics samples every 0.1 ms and takes the last entry into the 2 % band around the setpoint, before the load and after it, as the settling and recovery times:

ts = 0:1e-4:0.4

function metrics(sol)
    w = sol(ts; idxs = sys.inertia.w).u
    outside = abs.(w .- ω_ref) .> 0.02ω_ref       # outside the 2 % band
    before = ts .< t_load
    k_out = findlast(outside .& before)
    t_settle = k_out == findlast(before) ? NaN : ts[k_out + 1]
    k_back = findlast(outside .& .!before)
    t_back = k_back == length(ts) ? NaN : ts[k_back + 1] - t_load
    return (t_settle = t_settle, i_peak = maximum(sol(ts; idxs = sys.resistor.i).u),
            w_load = sol(t_load; idxs = sys.inertia.w), dip = ω_ref - minimum(w[.!before]), t_back = t_back)
end

m1 = metrics(sol)
println(isnan(m1.t_settle) ? "not in the 2 % band before the load step" : @sprintf("settled at %.1f ms", 1e3m1.t_settle))
@printf("current peak %.2f A\n", m1.i_peak)
@printf("speed at %.1f s: %.1f rad/s = %.1f %% of the setpoint; 1 - exp(-2) = %.1f %%\n",
        t_load, m1.w_load, 100m1.w_load / ω_ref, 100(1 - exp(-2)))

not in the 2 % band before the load step
current peak 1.97 A
speed at 0.2 s: 86.4 rad/s = 86.4 % of the setpoint; 1 - exp(-2) = 86.5 %

The low gain has not settled when the load arrives, and a rule predicts its speed at 0.2 s. \(JR/k^2 = 0.04\) s is the time constant of the motor's speed. The back-EMF brakes it through the resistor with \(k^2/R = 2.5\times10^{-3}\) N·m·s/rad, 250 times the friction \(d\), which is why \(J/d = 10\) s plays no part. Neglect the inductance and the motor is a single lag with this time constant. An integral time \(T\) equal to 0.04 s cancels that lag, and the controlled speed rises like one exponential with time constant \(T k / k_p\), 100 ms here. At 0.2 s, two time constants, that is \(1 - e^{-2} = 86.5\) % of the setpoint, and the solver says 86.4 %.

function motor_figure()
    fig = Figure(size = (770, 484))
    ax_w = Axis(fig[1, 1], ylabel = "ω / (rad/s)", limits = (nothing, (-5, 118)), yticklabelspace = 40.0)
    ax_i = Axis(fig[2, 1], xlabel = "t / s", ylabel = "i / A", yticklabelspace = 40.0)   # aligned y labels
    linkxaxes!(ax_w, ax_i); hidexdecorations!(ax_w, grid = false)
    hlines!(ax_w, [ω_ref], color = SECOND, linewidth = 1.2, linestyle = :dash)
    vlines!(ax_w, [t_load], color = MUTED, linewidth = 1.2, linestyle = :dash)
    vlines!(ax_i, [t_load], color = MUTED, linewidth = 1.2, linestyle = :dash)
    text!(ax_w, 0.3, 102, text = "setpoint", color = SECOND, align = (:left, :bottom))
    text!(ax_w, t_load + 0.005, 2, text = "load on", color = MUTED, align = (:left, :bottom))
    return fig, ax_w, ax_i
end

fig, ax_w, ax_i = motor_figure()
lines!(ax_w, ts, sol(ts; idxs = sys.inertia.w).u, color = INK)
lines!(ax_i, ts, sol(ts; idxs = sys.resistor.i).u, color = INK)
fig

Speed in rad/s (top) and armature current in A (bottom) over 0.4 s for the PI gain kp = 0.02. The speed creeps toward the dashed setpoint at 100 rad/s and is at 86 rad/s when the load switches on at 0.2 s; the current peaks near 2 A.

Step 5: Retune the controller with remake

The rule predicts the result before you solve: five times the gain gives a fifth of the time constant, 20 ms, and the speed comes within 2 % of the setpoint when \(e^{-t/20\,\text{ms}} = 0.02\), after \(\ln 50 \times 20\) ms. remake copies the problem with new parameter values and reuses the compiled functions:

@printf("predicted settling time %.1f ms\n", 1e3 * log(50) * 0.04k / 0.1)
prob2 = remake(prob; p = [sys.controller.k => 0.1])
sol2 = solve(prob2, Tsit5(); reltol = 1e-8, abstol = 1e-10)
m2 = metrics(sol2)
@printf("settled at %.1f ms; current peak %.2f A; peak voltage %.1f V\n",
        1e3m2.t_settle, m2.i_peak, maximum(sol2(ts; idxs = sys.source.V.u).u))
@printf("load step: dip %.1f rad/s, back in the band after %.1f ms; current with load %.2f A\n",
        m2.dip, 1e3m2.t_back, sol2(0.4; idxs = sys.resistor.i))

predicted settling time 78.2 ms
settled at 78.4 ms; current peak 9.07 A; peak voltage 10.0 V
load step: dip 3.1 rad/s, back in the band after 61.5 ms; current with load 0.62 A

The solver lands 0.2 ms from the prediction. The current peak is 9.07 A, 4.6 times the 1.97 A of the low gain, at a controller voltage of 10 V. When the load arrives the speed dips by 3.1 rad/s and is back in the band 61.5 ms later. The steady current with load is 0.62 A, which is the load torque plus friction, 0.031 N·m, divided by \(k\).

fig, ax_w, ax_i = motor_figure()
for (s, color) in [(sol, INK), (sol2, ACCENT)]
    lines!(ax_w, ts, s(ts; idxs = sys.inertia.w).u, color = color)
    lines!(ax_i, ts, s(ts; idxs = sys.resistor.i).u, color = color)
end
w_settle = sol2(m2.t_settle; idxs = sys.inertia.w)
scatter!(ax_w, [m2.t_settle], [w_settle], color = ACCENT, markersize = 10)
text!(ax_w, m2.t_settle + 0.005, w_settle - 3, text = @sprintf("settled in %.0f ms", 1e3m2.t_settle),
      color = ACCENT, align = (:left, :top))
k_peak = argmax(sol2(ts; idxs = sys.resistor.i).u)
scatter!(ax_i, [ts[k_peak]], [m2.i_peak], color = ACCENT, markersize = 10)
text!(ax_i, ts[k_peak] + 0.008, m2.i_peak, text = @sprintf("%.1f A", m2.i_peak),
      color = ACCENT, align = (:left, :center))
text!(ax_w, 0.25, 80, text = "kp = 0.02", color = INK, align = (:left, :top))
text!(ax_w, 0.1, 102, text = "kp = 0.1", color = ACCENT, align = (:left, :bottom))
fig

Speed in rad/s (top) and armature current in A (bottom) of a DC motor over 0.4 s for two PI gains. With kp = 0.1 the speed settles at 100 rad/s in 78 ms at a 9.1 A current peak; with kp = 0.02 it is still 14 rad/s short when the load switches on at 0.2 s.

Pitfalls

A missing ground. Leave the ground out of the circuit and mtkcompile refuses:

no_ground = [connections[1:3]; connect(emf.n, source.n); connections[5:end]]
@named broken = System(no_ground, t; systems = components)
try
    mtkcompile(broken)
catch e
    println(first(split(sprint(showerror, e), '\n')))
end

ExtraVariablesSystemException: The system is unbalanced. There are 24 highest order derivative variables and 23 equations.

mtkcompile counts after its first pass has removed the plain equalities, hence 24 variables against 23 equations instead of numbers near 74. Without a ground, every equation of the circuit is about differences of potential, and nothing fixes their level: the system is one equation short. Every electrical circuit needs one Ground, every rotational chain a Fixed somewhere (here at the damper's far end), and every signal input something connected to it.

Reading a variable by position. sol[2, :] looks like the speed if you think of the hand-written model, where the speed came second:

println(length(sol2.u[end]), " numbers in the state")
@printf("sol2[2, end] = %.1f, the shaft angle in rad; speed %.1f rad/s; EMF torque %.4f N·m\n",
        sol2[2, end], sol2[sys.inertia.w][end], sol2[sys.emf.flange.tau][end])

4 numbers in the state
sol2[2, end] = 37.8, the shaft angle in rad; speed 99.9 rad/s; EMF torque -0.0312 N·m

The state holds four numbers in the order mtkcompile chose, and the second one is the angle, 37.8 rad. That order can change with the next version of the library or the next component you add. The EMF torque, −0.0312 N·m into the EMF and so +0.0312 N·m on the shaft, is in no position at all. Index by symbol, and observed(sys) lists everything you can ask for.

API names from another ModelingToolkit version. Search results and older tutorials use ODESystem, structural_simplify, @mtkbuild, @mtkmodel, and ODEProblem(sys, u0map, tspan, pmap). Version 10 renamed the first three to System, mtkcompile, and @mtkcompile, and merged start values and parameters into one map, ODEProblem(sys, op, tspan). Version 11 deprecated @mtkmodel and moved it to the SciCompDSL package. The symptom is an UndefVarError or a deprecation warning on code copied from elsewhere. Compare pkgversion(ModelingToolkit), printed in the setup, with the version of the page you are reading.

Variations

A gear and a flexible shaft. Put an IdealGear(ratio = 10) between emf and a SpringDamper that drives a second Inertia. One more component per part and a few connect lines; now the shaft angle matters, because the spring's torque depends on it.
The motor heats up. Resistor(R = R, T_dep = true, alpha = 0.0039) adds a heat_port. Connect it to a HeatCapacitor and through a ThermalConductor to the ambient temperature, and the winding resistance rises as the motor works.
A supply limit. Add LimPI to the Blocks import and replace PI by LimPI(k = 0.1, T = 0.04, u_max = 6, Ta = 0.04). The controller asks for 10 V at the start and gets 6 V, so through 1 Ω the current stays under 6 A instead of 9.07 A. The 5.6 V the loaded motor needs at 100 rad/s is still within reach. The anti-windup keeps the integrator from running away while the voltage is pinned.
A Bode plot. Name the feedback signal as an analysis point, connect(speed_sensor.w, :w, feedback.input2). Then get_looptransfer(model, model.w; op = Dict(unknowns(sys) .=> 0.0)) returns the matrices A, B, C, D of the loop linearized at rest, with four states, and ControlSystemsBase turns them into a frequency response and stability margins. Leave out the operating point op and the call stops on a missing value for emf₊flange₊phi(t).

Cheat sheet

@component function Comp(; name, p)              # a component constructor
    @variables x(t); @parameters p = p           # changes in time; retunable constant
    System([D(x) ~ -p * x], t, [x], [p]; name, systems = [])   # extend(sys, base) inherits base
end
@named c = Comp(p = 1.0)                         # = Comp(; name = :c, p = 1.0)
@named model = System([connect(a.p, b.n, ground.g)], t; systems = [a, b, ground])   # potentials equal, flows sum to zero
sys = mtkcompile(model)                          # unknowns(sys), observed(sys), full_equations(sys)
prob = ODEProblem(sys, [sys.c.x => 1.0], (0.0, 1.0))   # start values and parameters in one map
sol = solve(prob, Tsit5()); sol[sys.a.v]         # index by symbol, observed variables too
prob2 = remake(prob; p = [sys.c.p => 2.0])       # retune without recompiling