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.
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:
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:
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:
# 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:
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:
using DyadInterface
using Plots
sol = ActiveSuspensionPIDAutotuningAnalysis()
PID Autotuning Analysis solution for PIDAutotuningAnalysis
Optimized parameters (parallel form): [763.1058525857778, 5432.011716657712, 714.3127705596897, 0.026504675367611308]
Cost: 9.001959047177901e-8
Disk-based margins:
Gain margins: 0.2, 5.3
Phase margin: 68.7°
Delay margin: 0.07277361380002935s
We may plot some of the available artifacts:
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:
artifacts(sol, :OptimizedParameters)
Row | kp_parallel | ki_parallel | kd_parallel | Kp_standard | Ti_standard | Td_standard | Tf |
---|---|---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | 763.106 | 5432.01 | 714.313 | 763.106 | 0.140483 | 0.93606 | 0.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:
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:
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:
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:
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)"
)