OK here is my go of making a PID controller. Notes:

- I was completely arbitrary in my choice of A B C D matrices and PID parameters
- Defaults for each component could be added
- Look how nice ModelingToolkit.jl makes it to access and plot particular variables by name in the solution

Any way this could be made cleaner, let me know.

```
using ModelingToolkit, OrdinaryDiffEq
## helper functions. reserving u_i and o_i as names of the ith input/output for each system
function get_input(sys, i)
return getproperty(sys, Symbol(:u, Symbol('₀' + i)))
end
get_inputs(sys, i) = get_input.((sys,), 1:i)
function get_output(sys, i)
return getproperty(sys, Symbol(:o, Symbol('₀' + i)))
end
get_outputs(sys, i) = get_output.((sys,), 1:i)
function linear_plant(A,B,C,D; name=:lp)
num_inputs = size(B,2)
num_states = size(A,1)
num_outputs = size(C,1)
@parameters t
dt = Differential(t)
@variables x[1:num_states](t)
@variables u[1:num_inputs](t)
@variables o[1:num_outputs](t)
eqs = vcat(
dt.(x) .~ A*x .+ B*u,
o .~ C*x .+ D*u
)
return ODESystem(eqs;
name=name
)
end
function PID_controller(Pgain, Igain, Dgain, RC, num_inputs; name=:pid)
@parameters t
D = Differential(t)
@variables u[1:num_inputs](t)
@variables o[1:num_inputs](t)
@variables int[1:num_inputs](t)
@variables mv[1:num_inputs](t)
eqs = vcat(
D.(int) .~ Igain.*u, #integrator
D.(mv) .~ u - RC.*(mv), # low pass filter
o .~ Dgain .* (u - RC.*(mv)) + Igain.*int + Pgain.*u
)
return ODESystem(eqs; name = name)
end
A = [0 1 0 0; -2 -2 0 0; 0 0 0 1; 0 0 -2 -2]
B = reshape([0 0 1. 0], 4,1)
C = reshape([1 0 0 0],1,4)
D = reshape([0], 1,1)
lp = linear_plant(A,B,C,D)
pid = PID_controller(1., 2.,1.5, 0.5, 1)
function connect(plant, ref, pid)
@parameters t
eqs = vcat(
get_input(pid,1) ~ ref(t) - get_output(plant,1),
get_output(pid,1) ~ get_input(plant, 1),
)
return ODESystem(eqs; systems=[plant,pid])
end
all = connect(lp,sin,pid)
final_sys = structural_simplify(all)
u0 = final_sys.states .=> zeros(length(final_sys.states))
prob = ODEProblem(final_sys, u0, (0.,100.), jac=true, sparse=true)
sol = solve(prob, Tsit5())
using Plots
plot(sol, vars = [pid.int₁])
sol[lp.u₁]
```