Design Optimization
In this demo we will show how to:
build a model from scratch using Dyad and standard library components
perform a design optimization using DyadModelOptimizer.jl to satisfy a given design specification.
In this example, we will design a stopper for a runaway train. We assume the following design specification:
The stopper needs to stop the train within 2 meters.
The stopper needs to do so at an acceleration less than 5*g.
By the end, we should have something that looks like this:
To get started, we shall launch the Dyad Studio
application on JuliaHub. When the app finishes launching, you will land in a VS Code environment. Run the "Dyad: Create Library" command in the command palette (shift + control + p), and create a new library named "DyadTutorial".
Building the System
We will design a stopper for a runaway train. The train is modeled as a series of connected train cars, and the stopper is modeled as a simple damped spring.
We can use standard components from the TranslationalComponents library.
using Pkg
Pkg.add(["TranslationalComponents"])
We will then create a feedback control system as given below:
Train Cars
First, let's design the train cars. This should be written to a file named dyad/train_cars.dyad
.
We model each train car as a mass, coupled via a spring damper on each side to another train car. The mass is modeled as a system with two flanges, mechanical connection points that transmit force.
component TrainCars
parameter v0::Velocity = 10
parameter a0::Acceleration = 0
parameter k::Real = 1e5
parameter d::Real = 1e5
car1 = TranslationalComponents.Mass(m = 1000, L = 0.9)
car2 = TranslationalComponents.Mass(m = 1000, L = 0.9)
car3 = TranslationalComponents.Mass(m = 1000, L = 0.9)
engine = TranslationalComponents.Mass(m = 2000, L = 0.9)
cx1 = TranslationalComponents.SpringDamper(c=k, d=d)
cx2 = TranslationalComponents.SpringDamper(c=k, d=d)
cx3 = TranslationalComponents.SpringDamper(c=k, d=d)
flange = Flange()
relations
initial car1.s = 0.45
initial car1.v = v0
initial car2.s = 1.55
initial car2.v = v0
initial car3.s = 2.55
initial car3.v = v0
initial engine.s = 3.55
initial engine.v = v0
initial engine.a = a0
initial cx1.initial_stretch = 0
initial cx2.initial_stretch = 0
connect(car1.flange_b, cx1.flange_a)
connect(car2.flange_a, cx1.flange_b)
connect(car2.flange_b, cx2.flange_a)
connect(car3.flange_a, cx2.flange_b)
connect(car3.flange_b, cx3.flange_a)
connect(engine.flange_a, cx3.flange_b)
connect(engine.flange_b, flange)
end
analysis TrainCarsTransientAnalysis
extends TransientAnalysis(stop=1)
model = TrainCars()
end
Solving the System
# using <LibraryName>
using DyadInterface
res = TrainCarsTransientAnalysis()
Transient Analysis Solution for TransientAnalysis
retcode: Success
Interpolation: 3rd order Hermite
t: 96-element Vector{Float64}:
0.0
0.011435302256567779
0.016665455302875883
0.025123478202032375
0.03372207804356194
0.0449417097668294
0.05763108530582586
0.07030355641107763
0.0820612394228276
0.09290329090485726
⋮
0.920015156770167
0.9306320935726707
0.9412490294041549
0.9518659650738024
0.9624829013907964
0.973099837626872
0.9837167745912125
0.9943337114746347
1.0
u: 96-element Vector{Vector{Float64}}:
[0.45, 1.55, 2.55, 3.55, 0.0, 0.0, 0.09999999999999963, 10.0, 10.0, 10.0, 10.0]
[0.5652277390613808, 1.6640298333589167, 2.6639741881807306, 3.66426667611368, 0.0, 0.0, 0.09999999999999963, 10.10848413620305, 10.000202119472036, 9.920734754933466, 9.985289494695726]
[0.6181923954043378, 1.7162809403684078, 2.715892707019466, 3.716453361175791, 0.0, 0.0, 0.09999999999999963, 10.140929869326488, 9.986677721054662, 9.92885829754259, 9.971767056038132]
[0.7040889509219107, 1.8007657178953225, 2.7998676984368367, 3.800725771423774, 0.0, 0.0, 0.09999999999999963, 10.166713241540766, 9.993953450977079, 9.926150933941418, 9.956591186770371]
[0.7915823704140935, 1.8867304716727067, 2.8852160886561746, 3.8862874857175607, 0.0, 0.0, 0.09999999999999963, 10.182309859865542, 10.001673349863912, 9.925003242440688, 9.94550677391493]
[0.9058992522528773, 1.9989873651060985, 2.99656941242002, 3.9978147292812363, 0.0, 0.0, 0.09999999999999963, 10.194111145204813, 10.009511232834283, 9.92406631844138, 9.936155651759764]
[1.0353102642190752, 2.1260280819528945, 3.12250742115847, 4.123854248980426, 0.0, 0.0, 0.09999999999999963, 10.200604271145911, 10.015870500779041, 9.923023068020703, 9.930251080027174]
[1.164612696752944, 2.2529430272394864, 3.2482864249245296, 4.24966783581846, 0.0, 0.0, 0.09999999999999963, 10.202052669472257, 10.02234420398979, 9.920309525033804, 9.927646800752076]
[1.2846029627293498, 2.370717414081785, 3.3649955760558505, 4.366373009137197, 0.0, 0.0, 0.09999999999999963, 10.200638592649941, 10.028652615064175, 9.91667390255003, 9.927017444867928]
[1.3952452933958472, 2.4793413725988493, 3.4726145703562064, 4.473981654445979, 0.0, 0.0, 0.09999999999999963, 10.199544060422419, 10.030854869437693, 9.916057712057667, 9.926771679041114]
⋮
[9.782340301092429, 10.762097383210097, 11.702076396148353, 12.70212187902873, 0.0, 0.0, 0.09999999999999963, 10.08675122411148, 10.014907687773274, 9.962336639540872, 9.968002224287186]
[9.889453950802421, 10.868358888249524, 11.80790550710215, 12.807943166239715, 0.0, 0.0, 0.09999999999999963, 10.085789628498004, 10.014808031988522, 9.96268024188463, 9.96836104881442]
[9.996557436514271, 10.974619331555473, 11.913738236898254, 12.913768232619866, 0.0, 0.0, 0.09999999999999963, 10.084838375667735, 10.014709506298804, 9.963020182943318, 9.968715967545071]
[10.103650875596397, 11.080878733185477, 12.019574554913342, 13.019597044997445, 0.0, 0.0, 0.09999999999999963, 10.083897353041305, 10.01461210081856, 9.96335649920377, 9.969067023468181]
[10.210734384239704, 11.187137113061054, 12.125414430937893, 13.125429570650578, 0.0, 0.0, 0.09999999999999963, 10.082966449238707, 10.0145158058079, 9.963689226734868, 9.969414259109262]
[10.31780806195817, 11.293394475599223, 12.231257819829159, 13.231265761978518, 0.0, 0.0, 0.09999999999999963, 10.082045556643319, 10.014420606030686, 9.964018406124502, 9.969757715600746]
[10.424872022629883, 11.399650840454528, 12.337104692190485, 13.337105587142858, 0.0, 0.0, 0.09999999999999963, 10.0811345662309, 10.014326492032932, 9.964344072620534, 9.970097434557816]
[10.531926362667432, 11.505906210971391, 12.442955002874744, 13.442948998609078, 0.0, 0.0, 0.09999999999999963, 10.080233372847834, 10.014233448563994, 9.964666266259323, 9.970433456164425]
[10.589050934005508, 11.562630146555346, 12.499435301384649, 13.499441809027243, 0.0, 0.0, 0.09999999999999963, 10.0819774208939, 10.00904672957471, 9.969361809304838, 9.969807020113278]
Now that we've solved this system, we can also plot it:
using CairoMakie
trains = artifacts(res, :SimplifiedSystem)
f, a, p = lines(
res;
idxs = [trains.car1.s, trains.car2.s, trains.car3.s, trains.engine.s],
)
axislegend(a; position = :lt)
f
This makes sense! In the absence of any external force, the train should continue to move at a constant velocity, which it does.
Stopper
Now, let's design the stopper,
component Stopper
extends TranslationalComponents.PartialTwoFlanges
parameter c::TranslationalSpringConstant
parameter d::DampingCoefficient
variable v_a::Velocity
variable v_b::Velocity
variable s_a::Dyad.Position
variable s_b::Dyad.Position
variable f::Force
relations
flange_a.s = s_a
flange_b.s = s_b
flange_a.f = f
flange_b.f = -f
D(s_a) = v_a
D(s_b) = v_b
f = if s_a >= s_b then
(v_a - v_b)*d + c*(s_a - s_b)
else
0
end
Building a complete system
Finally, we can unify these two components in a TrainCarWithStopper
component.
component TrainCarWithStopper
cars = TrainCars()
stopper = Stopper(c=1e1, d=1e4)
reference = TranslationalComponents.Fixed(s0=10)
relations
connect(cars.flange, stopper.flange_a)
connect(stopper.flange_b, reference.flange)
end
analysis TrainCarWithStopperTransientAnalysis
extends TransientAnalysis(stop=10)
model = TrainCarWithStopper()
end
Let's run this once, without any tuning, to see what happens:
res = TrainCarWithStopperTransientAnalysis()
Transient Analysis Solution for TransientAnalysis
retcode: Success
Interpolation: 3rd order Hermite
t: 946-element Vector{Float64}:
0.0
0.011435302256567779
0.016665455302875883
0.025123478202032375
0.03372207804356194
0.0449417097668294
0.05763108530582586
0.07030355641107763
0.0820612394228276
0.09290329090485726
⋮
9.924693940620172
9.935306588570777
9.945919239271495
9.956531893531185
9.967144551349847
9.977757211918624
9.98836987442866
9.998982540416785
10.0
u: 946-element Vector{Vector{Float64}}:
[1.55, 0.45, 2.55, 3.55, 0.0, 0.0, 0.09999999999999963, 10.0, 10.0, 10.0, 10.0]
[1.6640298333589167, 0.5652277390613808, 2.6639741881807306, 3.66426667611368, 0.0, 0.0, 0.09999999999999963, 10.10848413620305, 10.000202119472036, 9.920734754933466, 9.985289494695726]
[1.7162809403684078, 0.6181923954043378, 2.715892707019466, 3.716453361175791, 0.0, 0.0, 0.09999999999999963, 10.140929869326488, 9.986677721054662, 9.92885829754259, 9.971767056038132]
[1.8007657178953225, 0.7040889509219107, 2.7998676984368367, 3.800725771423774, 0.0, 0.0, 0.09999999999999963, 10.166713241540766, 9.993953450977079, 9.926150933941418, 9.956591186770371]
[1.8867304716727067, 0.7915823704140935, 2.8852160886561746, 3.8862874857175607, 0.0, 0.0, 0.09999999999999963, 10.182309859865542, 10.001673349863912, 9.925003242440688, 9.94550677391493]
[1.9989873651060985, 0.9058992522528773, 2.99656941242002, 3.9978147292812363, 0.0, 0.0, 0.09999999999999963, 10.194111145204813, 10.009511232834283, 9.92406631844138, 9.936155651759764]
[2.1260280819528945, 1.0353102642190752, 3.12250742115847, 4.123854248980426, 0.0, 0.0, 0.09999999999999963, 10.200604271145911, 10.015870500779041, 9.923023068020703, 9.930251080027174]
[2.2529430272394864, 1.164612696752944, 3.2482864249245296, 4.24966783581846, 0.0, 0.0, 0.09999999999999963, 10.202052669472257, 10.02234420398979, 9.920309525033804, 9.927646800752076]
[2.370717414081785, 1.2846029627293498, 3.3649955760558505, 4.366373009137197, 0.0, 0.0, 0.09999999999999963, 10.200638592649941, 10.028652615064175, 9.91667390255003, 9.927017444867928]
[2.4793413725988493, 1.3952452933958472, 3.4726145703562064, 4.473981654445979, 0.0, 0.0, 0.09999999999999963, 10.199544060422419, 10.030854869437693, 9.916057712057667, 9.926771679041114]
⋮
[12.597855537443857, 11.697856538268859, 13.49783870722118, 14.497804279991184, 0.0, 0.0, 0.09999999999999963, -0.005047123750323944, -0.005039703614550496, -0.005027699702341372, -0.004986381891587704]
[12.597802019286874, 11.697802998174641, 13.497785389693739, 14.49775135748227, 0.0, 0.0, 0.09999999999999963, -0.005045963358283595, -0.005038578226912877, -0.005026797558539928, -0.004985908429473178]
[12.59774851300069, 11.697749470310296, 13.497732081669154, 14.497698439959976, 0.0, 0.0, 0.09999999999999963, -0.005044815786498825, -0.0050374651520203855, -0.005025905098497863, -0.004985439752003508]
[12.597695018451239, 11.697695954536506, 13.497678783041241, 14.497645527369768, 0.0, 0.0, 0.09999999999999963, -0.005043680888384841, -0.005036364247527176, -0.005025022214252803, -0.0049849758052263205]
[12.59764153551002, 11.69764245071959, 13.497625493709034, 14.497592619661711, 0.0, 0.0, 0.09999999999999963, -0.00504255851832028, -0.005035275374606878, -0.005024148797504703, -0.004984516536106441]
[12.597588064054081, 11.697588958731476, 13.497572213576762, 14.497539716790467, 0.0, 0.0, 0.09999999999999963, -0.005041448531663794, -0.005034198397837556, -0.005023284739678637, -0.004984061892504511]
[12.597534603961913, 11.697535478445609, 13.497518942549782, 14.497486818711243, 0.0, 0.0, 0.09999999999999963, -0.005040350785395115, -0.005033133183355705, -0.005022429933394453, -0.004983611822853138]
[12.59748115510161, 11.69748200972512, 13.49746568052281, 14.497433925368107, 0.0, 0.0, 0.09999999999999963, -0.005039265140130506, -0.00503207959339233, -0.0050215842768341285, -0.004983166275228778]
[12.597476034644279, 11.69747688277579, 13.497460571808377, 14.497428855128346, 0.0, 0.0, 0.09999999999999963, -0.005038702551396517, -0.005033041787856357, -0.0050205661953714445, -0.004983293371690752]
sol = artifacts(res, :RawSolution)
train_with_stopper = artifacts(res, :SimplifiedSystem)
left_box_idxs = [getproperty(train_with_stopper.cars, Symbol("car", i)).flange_a.s for i in 1:3]
push!(left_box_idxs, train_with_stopper.cars.engine.flange_a.s)
right_box_idxs = [getproperty(train_with_stopper.cars, Symbol("car", i)).flange_b.s for i in 1:3]
push!(right_box_idxs, train_with_stopper.cars.engine.flange_b.s)
t = Observable(0.0)
boxes = lift(t) do t
starts = sol(t; idxs = left_box_idxs)
stops = sol(t; idxs = right_box_idxs)
Rect2f.(tuple.(starts, 0.0), tuple.(stops .- starts, 1.0))
end
f, a, p = poly(boxes)
ylims!(a, -4, 5)
xlims!(a, 0, 14)
f
We can also animate this:
record(f, "train_with_stopper.mp4", LinRange(0, 1, 100)) do ti
t[] = ti
a.title[] = string("t = ", round(ti; digits = 3))
end
"train_with_stopper.mp4"
Optimizing the System
Now, let's optimize the system to satisfy the design specification.
In order to optimize, we will use the DyadModelOptimizer.jl
package, which allows us to optimize a model by specifying an experiment design, constraints, a loss function (running cost), and tunable parameters.
In this case, we will be tuning the spring constant (c
) and damping coefficient (d
) of the stopper.
using DyadModelOptimizer
using OrdinaryDiffEqRosenbrock: Rodas5P
using ModelingToolkit # to manipulate and construct the Dyad system
First, we need to build the system. The @mtkbuild
macro is simply a convenience that will name the system, and structurally simplify it so it can be solved.
@mtkbuild sys = TrainCarWithStopper()
Now, we define our constraints. We want the acceleration of the last car to be less than 5g
, we want the stopper to stop the train before it gets compressed to less than 2 meters, and we want the train to completely stop at the end of the simulation.
In this case, our "loss function", or penalty, is higher the faster the engine is at the end of the simulation period.
We also only want to tune two variables - stopper.c
and stopper.d
. We allow them to range from 0 to 1,000,000,000 (1e9).
g = 9.807
constraints = [
-5*g ≲ sys.cars.car3.a
(sys.stopper.s_a - sys.stopper.s_b) ≲ 2
0 ≲ sys.cars.engine.v[end]
]
running_cost = abs(sys.cars.engine.v[end]);
tunables = [
sys.stopper.c => (0, 1e9),
sys.stopper.d => (0, 1e9),
]
2-element Vector{Pair{Symbolics.Num, Tuple{Int64, Float64}}}:
stopper₊c => (0, 1.0e9)
stopper₊d => (0, 1.0e9)
Now that we've defined our constraints and metric of success, we can wrap all of this into a DesignConfiguration
, which specifies the experiment design.
Then, we specify the algorithm to use for the optimization. In this case, we will use the SingleShooting
algorithm, which is a type of direct collocation method.
Finally. we calibrate the inverse problem and algorithm, and run a simulation.
designconfig = DesignConfiguration(sys;
alg = Rodas5P(),
tspan = (0.0, 2.0),
saveat = 0.0:0.01:2.0,
constraints = constraints,
running_cost = running_cost,
reduction = minimum
)
invprob = InverseProblem(designconfig, tunables)
alg = SingleShooting(maxiters = 1000)
r = calibrate(invprob, alg)
optimized_sol = simulate(designconfig, r)
retcode: Success
Interpolation: 1st order linear
t: 201-element Vector{Float64}:
0.0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
⋮
1.92
1.93
1.94
1.95
1.96
1.97
1.98
1.99
2.0
u: 201-element Vector{Vector{Float64}}:
[1.55, 0.45, 2.55, 3.55, 0.0, 0.0, 0.09999999999999963, 10.0, 10.0, 10.0, 10.0]
[1.6497907403134684, 0.5506685407687827, 2.649652922540951, 3.6499438981884014, 0.0, 0.0, 0.09999999999999963, 10.111889385977737, 9.974288208604568, 9.943031869925562, 9.98539526774607]
[1.749589107955689, 0.6520257809735258, 2.7489991289797104, 3.7496929910455434, 0.0, 0.0, 0.09999999999999963, 10.154369534589858, 9.986249288304899, 9.929377823873674, 9.965001676615788]
[1.8495124112534058, 0.7536918959619073, 2.8482717217850757, 3.849261985499813, 0.0, 0.0, 0.09999999999999963, 10.176789154618033, 9.997855840710779, 9.925843710326307, 9.949755647172438]
[1.949534529663358, 0.8555313523613428, 2.9475235647490967, 3.94870527661311, 0.0, 0.0, 0.09999999999999963, 10.189971340950338, 10.006027993242045, 9.924727672487437, 9.93963649666009]
[2.0496237313053594, 0.9574735022345243, 3.046768884714239, 4.048066940872948, 0.0, 0.0, 0.09999999999999963, 10.197755193829497, 10.011427803156293, 9.92442342412542, 9.933196789444391]
[2.149756745699252, 1.0594751124139954, 3.146013326884671, 4.147377407501051, 0.0, 0.0, 0.09999999999999963, 10.202112262021956, 10.014917201452032, 9.924513844672822, 9.92922834592659]
[2.2499178624577927, 1.161508373210003, 3.2452599314317516, 4.24665691645024, 0.0, 0.0, 0.09999999999999963, 10.204243342093934, 10.01713579759466, 9.924837351856038, 9.926891754227684]
[2.3500966859083325, 1.2635551650146672, 3.3445105589183024, 4.345918795079363, 0.0, 0.0, 0.09999999999999963, 10.204921423024302, 10.018516845865383, 9.925307388772126, 9.925627171169097]
[2.4502863753128032, 1.3656036618564062, 3.4437663873864564, 4.445171787722184, 0.0, 0.0, 0.09999999999999963, 10.204651877761075, 10.019347459964191, 9.925870490506382, 9.92506508588418]
⋮
[9.763393462485444, 8.869866256583167, 10.607087906259364, 11.501778236971804, 0.0, 0.0, 0.09999999999999963, -0.16202995515774607, -0.15529649741408116, -0.09847710413941166, 0.0071478435017350616]
[9.761841856410202, 8.868247358682716, 10.606102232412129, 11.501844386759014, 0.0, 0.0, 0.09999999999999963, -0.1617399133134658, -0.15501552946749367, -0.09864877041609484, 0.006090977463915978]
[9.760293330438927, 8.866631647026562, 10.605115104043302, 11.501900229784312, 0.0, 0.0, 0.09999999999999963, -0.1613932441318509, -0.1546809850293691, -0.09876846687544254, 0.005086100412588114]
[9.758748405729264, 8.865019672435801, 10.604127027671979, 11.501946274688489, 0.0, 0.0, 0.09999999999999963, -0.16099303495366898, -0.15429578261462631, -0.0988388267551703, 0.004130964587329432]
[9.757207574521463, 8.863411955141558, 10.603138483722738, 11.501983007830182, 0.0, 0.0, 0.09999999999999963, -0.1605422351735929, -0.15386271041073565, -0.09886236735083388, 0.003223417815626283]
[9.755671302015266, 8.86180898677155, 10.602149928195953, 11.502010894664085, 0.0, 0.0, 0.09999999999999963, -0.16004366481252338, -0.1533844343747553, -0.0988414971953558, 0.002361397642242054]
[9.754140027828646, 8.860211231894054, 10.601161793965836, 11.502030380811822, 0.0, 0.0, 0.09999999999999963, -0.1595000211782592, -0.15286350452469494, -0.09877852163736887, 0.001542926768070391]
[9.752614166077274, 8.85861912810203, 10.600174490851204, 11.502041892120399, 0.0, 0.0, 0.09999999999999962, -0.15891387922908665, -0.15230235528294986, -0.09867564314576896, 0.0007661128010330205]
[9.751094107100558, 8.857033087839984, 10.599188407151837, 11.502045835930565, 0.0, 0.0, 0.09999999999999963, -0.1582876994620577, -0.1517033129272507, -0.09853496791668122, 2.91428529325355e-5]
We can see from this simulation what the optimal values of c
and d
are.
Let's plot and animate this optimized system!
sol = optimized_sol
left_box_idxs = [getproperty(sys.cars, Symbol("car", i)).flange_a.s for i in 1:3]
push!(left_box_idxs, sys.cars.engine.flange_a.s)
right_box_idxs = [getproperty(sys.cars, Symbol("car", i)).flange_b.s for i in 1:3]
push!(right_box_idxs, sys.cars.engine.flange_b.s)
t = Observable(0.0)
boxes = lift(t) do t
starts = sol(t; idxs = left_box_idxs)
stops = sol(t; idxs = right_box_idxs)
Rect2f.(tuple.(starts, 0.0), tuple.(stops .- starts, 1.0))
end
f, a, p = poly(boxes)
ylims!(a, -4, 5)
xlims!(a, 0, 14)
f
Animating the tuned system
record(f, "optimized_train_with_stopper.mp4", LinRange(0, 1, 100)) do ti
t[] = ti
a.title[] = string("t = ", round(ti; digits = 3))
end
"optimized_train_with_stopper.mp4"
Conclusion
In this tutorial, we have shown how to design a stopper for a train using Dyad.