Control Library

The control library implements functionality for the design, simulation and optimization of control systems. Current functionality includes

Example: PID control with ModelingToolkit

This example will demonstrate a simple workflow for simulation of a control system.

The system model

We will consider a system consisting of two masses connected by a flexible element.

The code belows defines the system model, which also includes a damper to model dissipation of energy by the flexible element.

using ModelingToolkit, JuliaSimControls, OrdinaryDiffEq, Plots
@parameters t
Dₜ = Differential(t)

@parameters u(t) [input=true]  # Indicate that this is a controlled input
@parameters y(t) [output=true] # Indicate that this is a measured output

function Mass(; name, m = 1.0, p = 0, v = 0)
    ps = @parameters m=m
    sts = @variables pos(t)=p vel(t)=v
    eqs = Dₜ(pos) ~ vel
    ODESystem(eqs, t, [pos, vel], ps; name)
end

function Spring(; name, k = 1e4)
    ps = @parameters k=k
    @variables x(t)=0 # Spring deflection
    ODESystem(Equation[], t, [x], ps; name)
end

function Damper(; name, c = 10)
    ps = @parameters c=c
    @variables vel(t)=0
    ODESystem(Equation[], t, [vel], ps; name)
end

function SpringDamper(; name, k=false, c=false)
    spring = Spring(; name=:spring, k)
    damper = Damper(; name=:damper, c)
    compose(
        ODESystem(Equation[], t; name),
        spring, damper)
end


connect_sd(sd, m1, m2) = [sd.spring.x ~ m1.pos - m2.pos, sd.damper.vel ~ m1.vel - m2.vel]
sd_force(sd) = -sd.spring.k * sd.spring.x - sd.damper.c * sd.damper.vel

# Parameters
m1 = 1
m2 = 1
k = 1000 # Spring stiffness
c = 10   # Damping coefficient

@named mass1 = Mass(; m=m1)
@named mass2 = Mass(; m=m2)
@named sd = SpringDamper(; k, c)

function Model(u, d=0)
    eqs = [
        connect_sd(sd, mass1, mass2)
        Dₜ(mass1.vel) ~ ( sd_force(sd) + u) / mass1.m
        Dₜ(mass2.vel) ~ (-sd_force(sd) + d) / mass2.m
    ]
    @named _model = ODESystem(eqs, t; observed=[y~mass2.pos])
    @named model = compose(_model, mass1, mass2, sd)
end
Model (generic function with 2 methods)

We are now ready to simulate the system. We simulate the response of the system to a sinusoidal input and plot how the states evolve in time

model = Model(sin(30t))
sys = structural_simplify(model)
prob = ODEProblem(sys, Pair[], (0., 1.))
sol = solve(prob, Rosenbrock23())
plot(sol)

The next step is to introduce a controller and close the loop. We will look at how the closed-loop system behaves with a step on the reference. To handle the discontinuous reference step, we tell the solver to continue even though the step size gets small using force_dtmin=true. We also show how to connect a reference filter to shape the reference (feedforward). The closed-loop system will thus look like this

                                  d│                   
                                   │                   
    ┌───────┐        ┌────────┐    │   ┌────────┐      
 r  │       │      e │        │ u  ▼   │        │        y
────► filter├────+───►  PID   ├────+───►   P    ├───────┬───►
    │       │   -▲   │        │        │        │       │
    └───────┘    │   └────────┘        └────────┘       │
                 │                                      │
                 └──────────────────────────────────────┘
using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks: LimPID, FirstOrder
@variables u(t)=0 [input=true]
model = Model(u)
@named pid = LimPID(k=400, Ti=1, Td=1)
@named filt = FirstOrder(T=0.1)
connections = [
    filt.input.u  ~ t >= 1 # a reference step at t = 1
    pid.reference.u ~ filt.output.u
    pid.measurement.u ~ model.mass1.pos # we measure the position of the mass # TODO: this should use the declared output y above but this is currently not working
    pid.ctr_output.u   ~ model.u
]
closed_loop = structural_simplify(ODESystem(connections, t, systems=[model, pid, filt], name=:closed_loop))
prob = ODEProblem(closed_loop, Pair[], (0., 4.))
sol = solve(prob, Rosenbrock23(), dtmin=1e-15, force_dtmin=true);
plot(plot(sol, vars=[filt.y, model.mass1.pos, model.mass2.pos]), plot(sol, vars=[pid.measurement.u], title="Control signal"), legend=:bottomright)

We can also simulate the effects of actuator limitations, the simulation below introduces limits on the control signal of ±30 through the keyword argument u_max.

model = Model(u)
@named pid = LimPID(k=400, Ti=1, Td=1, u_max=30)
connections = [
    filt.input.u  ~ t >= 1 # a reference step at t = 1
    pid.reference.u ~ filt.output.u
    pid.measurement.u ~ model.mass1.pos   # we measure the position of the mass
    pid.ctr_output.u   ~ model.u
]
closed_loop = ODESystem(connections, t, systems=[model, pid, filt], name=:closed_loop)
closed_loop = structural_simplify(closed_loop)
prob = ODEProblem(closed_loop, Pair[], (0., 4.))
sol = solve(prob, Rosenbrock23(), dtmin=1e-15, force_dtmin=true);
plot(plot(sol, vars=[filt.y, model.mass1.pos, model.mass2.pos]), plot(sol, vars=[pid.measurement.u], title="Control signal"), legend=:bottomright)

Example: State feedback using ModelingToolkit

In this example we will instead simulate a state-feedback interconnection. We assume for now that all states are measurable and that a suitable feedback gain is given. The scenario this time is disturbance rejection, i.e., we simulate that a unit disturbance acts as a force on the second mass between 1 < t < 3.

using ModelingToolkitStandardLibrary.Blocks: MatrixGain
@variables u(t)=0 [input=true]
model = Model(u, (1 <= t) & (t <= 3)) # Model with load disturbance
lin, xs = JuliaSimControls.linearize(model, model.u, model.y)  # Linearize model
xs = filter(!ModelingToolkit.isinput, xs) # remove inputs from the state vector
# using ControlSystems
# L = lqr(ControlSystems.Continuous, lin.A, lin.B, 100I(4), I(1))[:]
L = [9.19578  6.82981  15.9409  -1.79881] # state-feedback gain
@named state_feedback = MatrixGain(L)
connections = [
    Symbolics.scalarize(state_feedback.input.u .~ xs)
    Symbolics.scalarize(model.u .~ -state_feedback.output.u) # negative feedback
]
closed_loop = ODESystem(connections, t, systems=[model, state_feedback], name=:closed_loop)
closed_loop = structural_simplify(closed_loop)
prob = ODEProblem(closed_loop, Pair[], (0., 10.))
sol = solve(prob, Rosenbrock23(), dtmin=1e-15, force_dtmin=true);
plot(plot(sol, vars=[model.mass1.pos, model.mass2.pos]), plot(sol, vars=[state_feedback.output.u], title="Control signal"))