Skip to content

Building an Inversion-Based Controller

In this tutorial we are going to assemeble an active suspension model using components from our standard libraries.

The Active Suspension Model

First, we shall make a very simple model of a car, using mass spring damper models. We can pull a mass spring model from the DyadExampleComponents library, along with a few other components we will need. Additionally, we will make use of components from the BlockComponents and TranslationalComponents libraries.

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

We will then create a feedback control system as given below:

Schematic in Dyad Builder

The above schematic is generated by the following Dyad code:

dyad
component ActiveSuspension
  wheel = DyadExampleComponents.MassSpringDamper(m=wheel_mass, d=wheel_damping, c=wheel_stiffness, g=-10, theta=pi/2, s0=wheel_initial_position)
  car_and_suspension = DyadExampleComponents.MassSpringDamper(m=car_mass, d=suspension_damping, c=suspension_stiffness, g=-10, theta=pi/2, s0=suspension_initial_position) 
  seat = DyadExampleComponents.MassSpringDamper(m=human_and_seat_mass, d=seat_damping, c=seat_stiffness, g=-10, theta=pi/2, s0=seat_initial_position)
  road_data = DyadExampleComponents.RoadData() 
  road = DyadExampleComponents.SimplePosition() 
  force = TranslationalComponents.Force()
  pid = BlockComponents.LimPID(k=Kp, Ti=Ti, Td=Td, Nd=Nd)
  set_point = BlockComponents.Constant(k=1.5)
  seat_pos = TranslationalComponents.PositionSensor()
  parameter wheel_mass::Dyad.Mass = 25
  parameter wheel_stiffness::TranslationalSpringConstant = 1e2
  parameter wheel_damping::DampingCoefficient = 1e4
  parameter car_mass::Dyad.Mass = 1000
  parameter suspension_stiffness::TranslationalSpringConstant = 1e4
  parameter suspension_damping::DampingCoefficient = 10
  parameter human_and_seat_mass::Dyad.Mass = 100
  parameter seat_stiffness::TranslationalSpringConstant = 1000
  parameter seat_damping::DampingCoefficient = 1
  parameter wheel_initial_position::Dyad.Position = 0.5
  parameter suspension_initial_position::Dyad.Position = 1
  parameter seat_initial_position::Dyad.Position = 1.5
  parameter Kp::Real = 20
  parameter Ti::Real = 5
  parameter Td::Real = 1
  parameter Nd::Real = 1
relations
  u: analysis_point(pid.y, force.f)
  y: analysis_point(seat_pos.s, pid.u_m)
  pid.u_ff = 0
  connect(road.s, road_data.y)
  connect(road.flange, wheel.flange_sd)
  connect(wheel.flange_m, car_and_suspension.flange_sd)
  connect(car_and_suspension.flange_m, seat.flange_sd, force.flange_a)
  connect(seat.flange_m, force.flange_b, seat_pos.flange)
  connect(pid.u_s, set_point.y)
  connect(pid.y, force.f)
  connect(seat_pos.s, pid.u_m)
end

In order to run this model, we should create an analysis:

dyad
analysis ActiveSuspensionTransient
  extends TransientAnalysis(alg="auto", abstol=10m, reltol=1m, start=0, stop=10)
  model = ActiveSuspension()
end

And we can interact with model in julia:

julia
# using <LibraryName>
using ModelingToolkit
using Plots
using DyadInterface
sol1 = ActiveSuspensionTransient()
as = artifacts(sol1, :SimplifiedSystem)
plot(sol1, idxs = [as.force.f],
  title = "Force supplied by actuator",
  xlabel = "Time (sec)",
  ylabel = "Force (N)"
)

Tuning the Controller

The controller in the above model is not tuned. We can specify another analysis to tune the controller:

dyad
analysis ActiveSuspensionPIDAutotuningAnalysis
  extends DyadExampleComponents.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 = ActiveSuspension()
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.1022575477714, 5432.062394914062, 714.3129261386894, 0.026504679064896338]
Cost: 9.001957474176607e-8

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

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.1025432.06714.313763.1020.1404810.9360640.0265047

Inverting the Model

Sometimes, we may want to ask "what if" questions of the model, for example when solving sizing problems. Consider a case where the control system perfectly cancels out the disturbance from the road. What would be force required from the actuator in order to do so? This knowledge could inform the design requirements of the actuator element.

In this case, we remove the controller, and instead add a simple block:

dyad
component ReverseCausality
    extends BlockComponents.SI2SO
relations
    u1 = u2
end

All this component does is constrain its inputs to be equal to each other. Why does this invert the model? We can constrain the seat position signal and the reference to be strictly equal to each other using this block. This forces the model to then calculate all other signals as a result of this constraint.

Schematic in Dyad Builder

The result is the following model:

corresponding to the following Dyad code:

dyad
component InvertedActiveSuspension
  wheel = DyadExampleComponents.MassSpringDamper(m=wheel_mass, d=wheel_damping, c=wheel_stiffness, g=-10, theta=pi/2, s0=wheel_initial_position)
  car_and_suspension = DyadExampleComponents.MassSpringDamper(m=car_mass, d=suspension_damping, c=suspension_stiffness, g=-10, theta=pi/2, s0=suspension_initial_position)
  seat = DyadExampleComponents.MassSpringDamper(m=human_and_seat_mass, d=seat_damping, c=seat_stiffness, g=-10, theta=pi/2, s0=seat_initial_position) 
  road_data = DyadExampleComponents.RoadData() 
  rev_causality = BlockComponents.ReverseCausality() 
  road = DyadExampleComponents.SimplePosition() 
  force = TranslationalComponents.Force()
  set_point = BlockComponents.Constant(k=1.5) 
  seat_pos = TranslationalComponents.PositionSensor()
  parameter wheel_mass::Dyad.Mass = 25
  parameter wheel_stiffness::TranslationalSpringConstant = 1e2
  parameter wheel_damping::DampingCoefficient = 1e4
  parameter car_mass::Dyad.Mass = 1000
  parameter suspension_stiffness::TranslationalSpringConstant = 1e4
  parameter suspension_damping::DampingCoefficient = 10
  parameter human_and_seat_mass::Dyad.Mass = 100
  parameter seat_stiffness::TranslationalSpringConstant = 1000
  parameter seat_damping::DampingCoefficient = 1
  parameter wheel_initial_position::Dyad.Position = 0.5
  parameter suspension_initial_position::Dyad.Position = 1
  parameter seat_initial_position::Dyad.Position = 1.5
  parameter kd::Real = 20
  parameter ki::Real = 0.2
  parameter kp::Real = 1
relations
  connect(road.s, road_data.y) 
  connect(road.flange, wheel.flange_sd) 
  connect(wheel.flange_m, car_and_suspension.flange_sd)
  connect(car_and_suspension.flange_m, seat.flange_sd, force.flange_a)
  connect(seat.flange_m, force.flange_b, seat_pos.flange)
  connect(seat_pos.s, rev_causality.u1) 
  connect(set_point.y, rev_causality.u2)
  connect(rev_causality.y, force.f) 
end

Similarly, we can run a TransientAnalysis for this new system too:

dyad
analysis InvertedActiveSuspensionTransient
  extends TransientAnalysis(alg="auto", abstol=10m, reltol=1m, start=0, stop=10)
  model = InvertedActiveSuspension()
end

And finally, we can plot the calculated value of the force:

julia
using ModelingToolkit
using Plots
using DyadInterface
sol1 = ActiveSuspensionTransient()
as = artifacts(sol1, :SimplifiedSystem)
sol2 = InvertedActiveSuspensionTransient()
ias = artifacts(sol2, :SimplifiedSystem)
plot(sol2, idxs = [ias.force.f],
  title = "Force supplied by actuator",
  label = "Inverted Model",
  xlabel = "Time (sec)",
  ylabel = "Force (N)"
)
plot!(sol1, idxs = [as.force.f],
  label = "Baseline Model",
  xlabel = "Time (sec)",
  ylabel = "Force (N)"
)

Layout Switch

Adjust the layout style of VitePress to adapt to different reading needs and screens.

Expand all
The sidebar and content area occupy the entire width of the screen.
Expand sidebar with adjustable values
Expand sidebar width and add a new slider for user to choose and customize their desired width of the maximum width of sidebar can go, but the content area width will remain the same.
Expand all with adjustable values
Expand sidebar width and add a new slider for user to choose and customize their desired width of the maximum width of sidebar can go, but the content area width will remain the same.
Original width
The original layout width of VitePress

Page Layout Max Width

Adjust the exact value of the page width of VitePress layout to adapt to different reading needs and screens.

Adjust the maximum width of the page layout
A ranged slider for user to choose and customize their desired width of the maximum width of the page layout can go.

Content Layout Max Width

Adjust the exact value of the document content width of VitePress layout to adapt to different reading needs and screens.

Adjust the maximum width of the content layout
A ranged slider for user to choose and customize their desired width of the maximum width of the content layout can go.

Spotlight

Highlight the line where the mouse is currently hovering in the content to optimize for users who may have reading and focusing difficulties.

ONOn
Turn on Spotlight.
OFFOff
Turn off Spotlight.