Skip to content

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:

  1. The stopper needs to stop the train within 2 meters.

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

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

dyad
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

julia
# 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:

julia
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,

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

dyad
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:

julia
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]
julia
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:

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

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

julia
@mtkbuild sys = TrainCarWithStopper()
dcars.car2.s(t)dt=cars.car2.v(t)dcars.car1.s(t)dt=cars.car1.v(t)dcars.car3.s(t)dt=cars.car3.v(t)dcars.engine.s(t)dt=cars.engine.v(t)dcars.cx1.initial_stretch(t)dt=0dcars.cx2.initial_stretch(t)dt=0dcars.cx3.initial_stretch(t)dt=0dcars.car1.v(t)dt=cars.car1.a(t)dcars.car2.v(t)dt=cars.car2.a(t)dcars.car3.v(t)dt=cars.car3.a(t)dcars.engine.v(t)dt=cars.engine.a(t)

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

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

julia
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!

julia
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

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