Introduction to Model Design in JuliaSim

In this demo we will show how to

  • build a model from scratch using JuliaSim and standard library components.
  • perform a design optimization using JuliaSimModelOptimizer.jl to satisfy a given design specification.

In this example, we will design a stopper for a runaway train. We shall assume the following design spec:

  1. The stopper needs to stop the train within a stroke/clearance of 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 JuliaSim IDE application on JuliaHub. When the app finishes launching, you will land on a code environment with an open Julia terminal. Now, we can paste the following code into the editor, save the file, and execute the following code inline by hitting Shift+Enter.

Building the System

To build the model, we need to design two basic components:

  1. The stopper itself, which we shall abstract as a mass-spring-damper system and,
  2. A model of a train, so we can study the interaction of the train with the stopper.

We shall start with a train.

Building a Train Model

We consider a train to be a series of train cars, connected by a coupler. First, we shall model a train car as a mass, with a velocity v and acceleration a. To validate this component, let us instantiate it with constant velocity, and plot position as a function of time. We should expect to see a straight line with constant slope.

using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t
using ModelingToolkitStandardLibrary.Mechanical.Translational
using DifferentialEquations
using Plots

@mtkmodel TrainCar begin
    @parameters begin
        m
    end
    @variables begin
        x_a(t)
        x_b(t)
        v(t)
        a(t)=0
    end
    @components begin
        port_a = MechanicalPort(;v)
        port_b = MechanicalPort(;v)
        mass = Mass(;m,v)
    end
    @equations begin
        D(x_a) ~ port_a.v
        D(x_b) ~ port_b.v
        v ~ port_a.v
        D(v) ~ a
        connect(port_a, mass.flange, port_b)
    end
end

@mtkbuild car = TrainCar(;m=100, v=1, x_a=0, x_b=0.1)
prob = ODEProblem(car, [], (0,1))
sol = solve(prob)
Plots.plot(sol; idxs=[car.x_a, car.x_b], xlabel="time [s]", ylabel="x [m]", title = "Car Displacement with Constant Velocity")

We can visualize the TrainCar component using the package ModelingToolkitDesigner. As we expect, it is simply a mass with a velocity, acceleration and two ports. 

using ModelingToolkitDesigner
fig = ModelingToolkitDesigner.view(ODESystemDesign(TrainCar(;name=:car), "."), false)

Now that we have one car component, we can make a train with several cars. But first, let's make a coupler that connect the cars. We'll use a spring-damper to "model" a connector.

@mtkmodel TrainCarCoupler begin
    @parameters begin
        k
        d
        v
    end
    @variables begin
    end
    @components begin
        port_a = MechanicalPort(;v)
        port_b = MechanicalPort(;v)
        spring = Spring(;k,delta_s=0)
        damper = Damper(;d,v=0,f=0)
    end
    @equations begin
        connect(port_a, spring.flange_a, damper.flange_a)
        connect(spring.flange_b, damper.flange_b, port_b)
    end
end

Finally, we can make our train model by just connecting a few cars and couplers together, with the connect statement in ModelingToolkit. Finally, we shall perform a quick sanity check on this sytem as well, initializing all the cars with a constant velocity.

@mtkmodel Train begin
    @parameters begin
        v=10
        k=1e5
        d=1e5
    end
    @variables begin 
    end
    @components begin
        car1 = TrainCar(;m=1000, v, x_a=0, x_b=0.9)
        cx1 = TrainCarCoupler(;k,d,v)

        car2 = TrainCar(;m=1000, v, x_a=1.1, x_b=1.9)
        cx2 = TrainCarCoupler(;k,d,v)

        car3 = TrainCar(;m=1000, v, x_a=2.1, x_b=2.9)
        cx3 = TrainCarCoupler(;k,d,v)

        engine = TrainCar(;m=2000, v, x_a=3.1, x_b=3.9)

        port = MechanicalPort(;v,f=0)
    end
    @equations begin
        connect(car1.port_b, cx1.port_a)
        connect(car2.port_a, cx1.port_b)
        connect(car2.port_b, cx2.port_a)
        connect(car3.port_a, cx2.port_b)
        connect(car3.port_b, cx3.port_a)
        connect(engine.port_a, cx3.port_b)
        connect(engine.port_b, port)
    end
end

@mtkbuild train = Train()
prob = ODEProblem(train, [], (0, 1))
sol = solve(prob)
Plots.plot(sol; idxs=[train.car1.x_a, train.car2.x_a, train.car3.x_a, train.engine.x_a], title = "Train displacement with constant velocity")

And finally, visualize the Train model.

fig = ModelingToolkitDesigner.view(ODESystemDesign(Train(;name=:train), "."), false)

Now that we have designed the train, it's time to design the stopper. The physics of the stopper component are straightforward: After the train makes contact, the force $F$ exerted by stopper is given by the following equation

\[\begin{aligned} F &= v*d + k*x \end{aligned}\]

where:

  • k is the spring constant.
  • d is the damping coefficient.
  • v is the relative velocity.
  • x is the relative position.
@mtkmodel Stopper begin
    @parameters begin
        k
        d
    end
    @variables begin
        v_a(t)
        v_b(t)
        x_a(t)
        x_b(t)
        f(t)=0
    end
    @components begin
        port_a = MechanicalPort(;v=v_a,f=+f)
        port_b = MechanicalPort(;v=v_b,f=-f)
    end
    @equations begin
        # differentials
        D(x_a) ~ v_a
        D(x_b) ~ v_b

        # connectors
        port_a.v ~ v_a
        port_b.v ~ v_b
        port_a.f ~ +f
        port_b.f ~ -f

        # physics
        f ~ ifelse(x_a >= x_b, (v_a - v_b)*d + k*(x_a - x_b), 0)
    end
end

Now, having built the stopper, we can assemble the system, and examine the behaviour of the system. 

@mtkmodel System begin
    @parameters begin
        v=10
    end
    @variables begin
    end
    @components begin
        train = Train(;v)
        stopper = Stopper(;k=1e1, d=1e4, v_a=v, v_b=0, x_a=3.9, x_b=10, f=0)
        reference = Fixed()
    end
    @equations begin
        connect(train.port, stopper.port_a)
        connect(stopper.port_b, reference.flange)
    end
end

@mtkbuild sys = System()
prob = ODEProblem(sys, [], (0, 5))
sol = solve(prob)

function plot_sys(sol)
    Plots.plot(sol; idxs=(sys.stopper.x_a - sys.stopper.x_b)*(sys.stopper.x_a > sys.stopper.x_b) , 
        label = "Stopper Compression")
    Plots.plot!(sol; idxs=sys.train.engine.v, label= "Engine Car Velocity")
    Plots.plot!(sol; idxs=sys.train.car1.a/9.807, label= "Car1 acc. / g")
    Plots.plot!(sol; idxs=sys.train.car2.a/9.807, label= "Car2 acc. / g")
    Plots.plot!(sol; idxs=sys.train.car3.a/9.807, label= "Car3 acc. / g")
    Plots.hline!([2], linestyle = :dash, label = "2", lw = 0.7)
    Plots.hline!([-5], linestyle = :dash, label = "-5", lw = 0.7)
    Plots.plot!(title= "Dynamics on Contact with Stopper")
end
plot_sys(sol)

Let us unpack the plot above. The red line indicates the velocity of the engine car which, as we expect, comes down to zero after contact with the stopper. The acceleration of the train cars, denoted by the green, pink and gold lines are well within 5*g. However, the relative position of the stopper, denoted by the blue plot, far overshoots our desired spec of 2 meters.

To understand what's happening, let us visualize the dynamic behaviour of the system.

Additionally, let us visualize the entire system architecture.

fig = ModelingToolkitDesigner.view(ODESystemDesign(System(;name=:sys), "."), false)

So how do we deal with this mismatch from our design? We can try and tune the model by performing a design optimization routine. We do this using JuliaSimModelOptimizer. Using, this package we shall define a DesignConfiguration which takes in constraints, which we define symbolically as follows.

We contain the decelerration of Car3 to be within 5g's (given that car 3 has the highest decelerration), specify that the stopper compression should be within 2 metres, and additionally ensure that the train comes to a halt by the end of the timespan.

using JuliaSimModelOptimizer

g = 9.807
cnst = [
    -5*g ≲ sys.train.car3.a
    (sys.stopper.x_a - sys.stopper.x_b) ≲ 2
    0 ≲ sys.train.engine.v[end] 
]

designconfig = DesignConfiguration(sys;
    alg = Rodas4(),
    abstol = 1e-6,
    reltol = 1e-6,
    tspan = (0.0, 2.0),
    saveat = 0.0:0.01:2.0,
    constraints=cnst,
    running_cost = abs(sys.train.engine.v[end]),
    reduction = minimum
)

Our design variables are the contact stiffness $k$ and damping coefficient $d$ of the stopper. We give the package a large search space like follows and call calibrate. 

tunables =[
    sys.stopper.k => (0, 1e9),
    sys.stopper.d => (0, 1e9),
]
invprob = InverseProblem(designconfig, tunables)
alg = SingleShooting(maxiters = 1000, optimizer = IpoptOptimizer(verbose = true))
r = calibrate(invprob, alg)
optimized_sol = simulate(designconfig, r)
plot_sys(optimized_sol)

After optimizing our system, we can see that the acceleration has reduced, and the clearance is now at 2 meters. We have now satisfied our spec, and finished our model design.