Dyad

Appearance

TUTORIAL

Building an Inversion-Based Controller ​

Download as a Dyad projectActiveSuspensionID.zipOpen in Dyad Studio

In this tutorial we will model and simulate an active suspension model using components from our standard libraries.

To get started, create a new component library, as indicated on the Getting started page.

The Active Suspension Model ​

We will use components from our standard libraries: DyadExampleComponents, TranslationalComponents, and BlockComponents. Additionally, we will import DyadControlSystems. In the Julia REPL, run the following commands:

julia
using Pkg
Pkg.add(["DyadExampleComponents", "TranslationalComponents", "BlockComponents", "DyadControlSystems"])

The diagram below is the graphical representation of the ActiveSuspension model. It represents a vehicle moving over a bumpy road. On the vehicle is a PID controller controlling an actuating force which dampens the seat response to the reference position of 1.5. The vehicle consists of a seat, suspension system, and wheel. Each of these three elements is modeled as a MassSpringDamper from the DyadExampleComponents library.

Let's define this model in our library. In the dyad folder, create a new file named ActiveSuspension.dyad and include the code below.

To simulate this model, we create an analysis that extends the TransientAnalysis type and uses our ActiveSuspension model:

dyad
analysis ActiveSuspensionTransient
  extends TransientAnalysis(alg=ODEAlg.Auto(), abstol=10m, reltol=1m, start=0, stop=10, dtmax=10m)
  model = ActiveSuspension()
end

Let's run this analysis and plot two signals: the seat position and force response. Run the following commands in the Julia REPL:

julia
# using <LibraryName>
using ModelingToolkit
using Plots
using DyadInterface

sol1 = ActiveSuspensionTransient()
as = artifacts(sol1, :SimplifiedSystem)

plot(sol1, idxs = [as.seat_pos.s],
  title = "Seat position and actuator force",
  label = "Seat position",
  xlabel = "Time (sec)",
  ylabel = "Position"
)
plot!(sol1, idxs = [as.force.f],
  label = "Force",
  ylabel = "Force (N)"
)

Tuning the Controller ​

The controller in the above model is not tuned. We can specify another analysis to tune the controller. Since the PID autotuning analysis internally opens the loop at the measurement point y, the connection from seat_pos.s to pid.u_m is broken during linearization, and pid.u_m becomes a free input that needs an operating-point value. We also need to pin the algebraic spring-damper forces and the internal PID-block outputs so that the linearization solver does not encounter cyclic guesses. We use Dyad's guess keyword to provide these operating-point values:

dyad
component ActiveSuspensionForAutotuning
  extends ActiveSuspension()
relations
  guess pid.u_m = 1.5
  guess wheel.spring_damper.f_c = 0
  guess car_and_suspension.spring_damper.f_c = 0
  guess seat.spring_damper.f_c = 0
  guess pid.add_pid.y = 0
  guess pid.derivative.y = 0
end

analysis ActiveSuspensionPIDAutotuningAnalysis
  extends DyadControlSystems.PIDAutotuningAnalysis(
    measurement = "y",
    control_input = "u",
    step_input = "u",
    step_output = "y",
    Ts = 0.01,
    duration = 10.0,
    Ms = 1.2,
    Mt = 1.1,
    Mks = 3e4,
    ki_lb = 2000,
    wl = 1e-2,
    wu = 1e3
  )
  model = ActiveSuspensionForAutotuning()
end

Once we specify the analysis declaratively, we can run it imperatively:

julia
using DyadInterface
using Plots
sol = ActiveSuspensionPIDAutotuningAnalysis()
PID Autotuning Analysis solution for PIDAutotuningAnalysis
Optimized parameters (parallel form): [763.1073446297678, 5431.9891078208675, 714.3126666827111, 0.026504673820710346]
Cost: 9.001960976281153e-8

Disk-based margins:
	Gain margins: 0.2, 5.3
	Phase margin: 68.7°
	Delay margin: 0.07277380397091959s

We may plot some of the available artifacts:

julia
figs  = artifacts(sol, :SensitivityFunctions)
fign  = artifacts(sol, :NoiseSensitivityAndController)
figr  = artifacts(sol, :OptimizedResponse)
figny = artifacts(sol, :NyquistPlot)

plot!(figs, ylims=(1e-2, Inf))
plot!(figny, legend=:bottomleft)

plot(figs, fign, figr, figny,
    link = :none,
    size = (600, 600),
)

We may also ask for the optimized controller parameters:

julia
artifacts(sol, :OptimizedParameters)
1×7 DataFrame
Rowkp_parallelki_parallelkd_parallelKp_standardTi_standardTd_standardTf
Float64Float64Float64Float64Float64Float64Float64
1763.1075431.99714.313763.1070.1404840.9360580.0265047

Inverting the Model ​

Say we need to know the force requirement of the actuator if the control system were to keep the seat position precisely at the reference value of 1.5 as the vehicle moves over the bumpy road. In many workflows, we would first need to autotune the PID controller to determine the appropriate controller gain values before we are able to calculate the corresponding force signal. In a acausal modeling framework like Dyad, however, we are able to solve for the force signal directly without having to tune the PID controller.

To do this, we will create a component, ReverseCausality, that extends from the SI2SO component from the BlockComponents library. ReverseCausality will constrain its inputs to be equal to one another. This will allow us to constrain the seat position to the reference value in our model. When we simulate, the model will then calculate all other signals such that this constraint is met. Another way to think of it: rather than having the PID controller drive the actuating force based on the PID gains, we will have the ReverseCausality block drive the actuating force based on the constraint that the two block inputs must be equal to one another.

In the dyad folder, create a new file named InvertedActiveSuspension.dyad and include the code below.

We will modify our model to remove the PID controller and incorporate the ReverseCausality component. Below is the Dyad implementation for this InvertedActiveSuspension model — include it in the InvertedActiveSuspension.dyad file:

Similarly to before, we create an analysis that extends the TransientAnalysis type and uses our InvertedActiveSuspension model:

dyad
analysis InvertedActiveSuspensionTransient
  extends TransientAnalysis(alg=ODEAlg.Auto(), abstol=10m, reltol=1m, start=0, stop=10, dtmax=10m)
  model = InvertedActiveSuspension()
end

Finally, we can compare the seat positions and force responses between the nominal and the inverted models. First, run both analyses:

julia
using ModelingToolkit
using Plots
using DyadInterface

sol1 = ActiveSuspensionTransient()
as = artifacts(sol1, :SimplifiedSystem)

sol2 = InvertedActiveSuspensionTransient()
ias = artifacts(sol2, :SimplifiedSystem)
┌ Warning: compile_lustre: no continuous time partition found; discrete-only systems are not currently supported
└ @ SynchToolkit ~/.julia/packages/SynchToolkit/Xf0NA/src/SynchToolkit.jl:294
┌ Warning: compile_lustre: no continuous time partition found; discrete-only systems are not currently supported
└ @ SynchToolkit ~/.julia/packages/SynchToolkit/Xf0NA/src/SynchToolkit.jl:294

Compare the seat positions:

julia
plot(sol1, idxs = [as.seat_pos.s], label = "Baseline Model",
  title = "Seat position", xlabel = "Time (sec)", ylabel = "Position")
plot!(sol2, idxs = [ias.seat_pos.s], label = "Inverted Model")

Compare the actuator force:

julia
plot(sol1, idxs = [as.force.f], label = "Baseline Model",
  title = "Force supplied by actuator", xlabel = "Time (sec)", ylabel = "Force (N)")
plot!(sol2, idxs = [ias.force.f], label = "Inverted Model")