Dyad

Appearance

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)"
)