Component-Based Modeling a Spring-Mass System

In this tutorial we will build a simple component-based model of a spring-mass system. A spring-mass system consists of one or more masses connected by springs. Hooke's law gives the force exerted by a spring when it is extended or compressed by a given distance. This specifies a differential-equation system where the acceleration of the masses is specified using the forces acting on them.

Copy-Paste Example

using ModelingToolkit, Plots, DifferentialEquations, LinearAlgebra
using Symbolics: scalarize

@variables t
D = Differential(t)

function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
    ps = @parameters m=m
    sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u
    eqs = scalarize(D.(pos) .~ v)
    ODESystem(eqs, t, [pos..., v...], ps; name)
end

function Spring(; name, k = 1e4, l = 1.)
    ps = @parameters k=k l=l
    @variables x(t), dir(t)[1:2]
    ODESystem(Equation[], t, [x, dir...], ps; name)
end

function connect_spring(spring, a, b)
    [
        spring.x ~ norm(scalarize(a .- b))
        scalarize(spring.dir .~ scalarize(a .- b))
    ]
end

spring_force(spring) = -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l)  ./ spring.x

m = 1.0
xy = [1., -1.]
k = 1e4
l = 1.
center = [0., 0.]
g = [0., -9.81]
@named mass = Mass(m=m, xy=xy)
@named spring = Spring(k=k, l=l)

eqs = [
    connect_spring(spring, mass.pos, center)
    scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
]

@named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], [])
@named model = compose(_model, mass, spring)
sys = structural_simplify(model)

prob = ODEProblem(sys, [], (0., 3.))
sol = solve(prob, Rosenbrock23())
plot(sol)

Explanation

Building the components

For each component we use a Julia function that returns an ODESystem. At the top, we define the fundamental properties of a Mass: it has a mass m, a position pos and a velocity v. We also define that the velocity is the rate of change of position with respect to time.

function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
    ps = @parameters m=m
    sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u
    eqs = scalarize(D.(pos) .~ v)
    ODESystem(eqs, t, [pos..., v...], ps; name)
end

Mass (generic function with 1 method)

Note that this is an incompletely specified ODESystem. It cannot be simulated on its own since the equations for the velocity v[1:2](t) are unknown. Notice the addition of a name keyword. This allows us to generate different masses with different names. A Mass can now be constructed as:

Mass(name = :mass1)

\[ \begin{align} \mathrm{\frac{d}{d t}}\left( pos(t)_1 \right) =& v(t)_1 \\ \mathrm{\frac{d}{d t}}\left( pos(t)_2 \right) =& v(t)_2 \end{align} \]

Or using the @named helper macro

@named mass1 = Mass()

\[ \begin{align} \mathrm{\frac{d}{d t}}\left( pos(t)_1 \right) =& v(t)_1 \\ \mathrm{\frac{d}{d t}}\left( pos(t)_2 \right) =& v(t)_2 \end{align} \]

Next we build the spring component. It is characterised by the spring constant k and the length l of the spring when no force is applied to it. The state of a spring is defined by its current length and direction.

function Spring(; name, k = 1e4, l = 1.)
    ps = @parameters k=k l=l
    @variables x(t), dir(t)[1:2]
    ODESystem(Equation[], t, [x, dir...], ps; name)
end

Spring (generic function with 1 method)

We now define functions that help construct the equations for a mass-spring system. First, the connect_spring function connects a spring between two positions a and b. Note that a and b can be the pos of a Mass, or just a fixed position such as [0., 0.]. In that sense, the length of the spring x is given by the length of the vector dir joining a and b.

function connect_spring(spring, a, b)
    [
        spring.x ~ norm(scalarize(a .- b))
        scalarize(spring.dir .~ scalarize(a .- b))
    ]
end

connect_spring (generic function with 1 method)

Lastly, we define the spring_force function that takes a spring and returns the force exerted by this spring.

spring_force(spring) = -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l)  ./ spring.x

spring_force (generic function with 1 method)

To create our system, we will first create the components: a mass and a spring. This is done as follows:

m = 1.0
xy = [1., -1.]
k = 1e4
l = 1.
center = [0., 0.]
g = [0., -9.81]
@named mass = Mass(m=m, xy=xy)
@named spring = Spring(k=k, l=l)

\[ \begin{align} \end{align} \]

We can now create the equations describing this system, by connecting spring to mass and a fixed point.

eqs = [
    connect_spring(spring, mass.pos, center)
    scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
]

\begin{align} \mathrm{spring{+}x}\left( t \right) =& \sqrt{\left|mass{+}pos(t)1\right|^{2} + \left|mass{+}pos(t)2\right|^{2}} \ spring{+}dir(t)1 =& mass{+}pos(t)1 \ spring{+}dir(t)2 =& mass{+}pos(t)2 \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)1 \right) =& \frac{ - spring{+}k \left( - spring{+}l + \mathrm{spring{+}x}\left( t \right) \right) spring{+}dir(t)1}{mass{+}m \mathrm{spring{+}x}\left( t \right)} \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)2 \right) =& -9.81 + \frac{ - spring{+}k \left( - spring{+}l + \mathrm{spring{+}x}\left( t \right) \right) spring{+}dir(t)2}{mass{+}m \mathrm{spring_{+}x}\left( t \right)} \end{align}

Finally, we can build the model using these equations and components.

@named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], [])
@named model = compose(_model, mass, spring)

\[ \begin{align} \mathrm{spring_{+}x}\left( t \right) =& \sqrt{\left|mass_{+}pos(t)_1\right|^{2} + \left|mass_{+}pos(t)_2\right|^{2}} \\ spring_{+}dir(t)_1 =& mass_{+}pos(t)_1 \\ spring_{+}dir(t)_2 =& mass_{+}pos(t)_2 \\ \mathrm{\frac{d}{d t}}\left( mass_{+}v(t)_1 \right) =& \frac{ - spring_{+}k \left( - spring_{+}l + \mathrm{spring_{+}x}\left( t \right) \right) spring_{+}dir(t)_1}{mass_{+}m \mathrm{spring_{+}x}\left( t \right)} \\ \mathrm{\frac{d}{d t}}\left( mass_{+}v(t)_2 \right) =& -9.81 + \frac{ - spring_{+}k \left( - spring_{+}l + \mathrm{spring_{+}x}\left( t \right) \right) spring_{+}dir(t)_2}{mass_{+}m \mathrm{spring_{+}x}\left( t \right)} \\ \mathrm{\frac{d}{d t}}\left( mass_{+}pos(t)_1 \right) =& mass_{+}v(t)_1 \\ \mathrm{\frac{d}{d t}}\left( mass_{+}pos(t)_2 \right) =& mass_{+}v(t)_2 \end{align} \]

We can take a look at the equations in the model using the equations function.

equations(model)

\begin{align} \mathrm{spring{+}x}\left( t \right) =& \sqrt{\left|mass{+}pos(t)1\right|^{2} + \left|mass{+}pos(t)2\right|^{2}} \ spring{+}dir(t)1 =& mass{+}pos(t)1 \ spring{+}dir(t)2 =& mass{+}pos(t)2 \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)1 \right) =& \frac{ - spring{+}k \left( - spring{+}l + \mathrm{spring{+}x}\left( t \right) \right) spring{+}dir(t)1}{mass{+}m \mathrm{spring{+}x}\left( t \right)} \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)2 \right) =& -9.81 + \frac{ - spring{+}k \left( - spring{+}l + \mathrm{spring{+}x}\left( t \right) \right) spring{+}dir(t)2}{mass{+}m \mathrm{spring{+}x}\left( t \right)} \ \mathrm{\frac{d}{d t}}\left( mass{+}pos(t)1 \right) =& mass{+}v(t)1 \ \mathrm{\frac{d}{d t}}\left( mass{+}pos(t)2 \right) =& mass{+}v(t)_2 \end{align}

The states of this model are:

states(model)

7-element Vector{Term{Real}}:
 spring₊x(t)
 (spring₊dir(t))[1]
 (spring₊dir(t))[2]
 (mass₊pos(t))[1]
 (mass₊pos(t))[2]
 (mass₊v(t))[1]
 (mass₊v(t))[2]

And the parameters of this model are:

parameters(model)

3-element Vector{Any}:
 mass₊m
 spring₊k
 spring₊l

Simplifying and solving this system

This system can be solved directly as a DAE using one of the DAE solvers from DifferentialEquations.jl. However, we can symbolically simplify the system first beforehand. Running structural_simplify eliminates unnecessary variables from the model to give the leanest numerical representation of the system.

sys = structural_simplify(model)
equations(sys)

\begin{align} \mathrm{\frac{d}{d t}}\left( mass{+}pos(t)1 \right) =& mass{+}v(t)1 \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)1 \right) =& mass{+}pos(t)[1]ˍtt \ \mathrm{\frac{d}{d t}}\left( mass{+}pos(t)2 \right) =& mass{+}v(t)2 \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)2 \right) =& mass{+}pos(t)[2]ˍtt \end{align}

We are left with only 4 equations involving 4 state variables (mass.pos[1], mass.pos[2], mass.v[1], mass.v[2]). We can solve the system by converting it to an ODEProblem. Some observed variables are not expanded by default. To view the complete equations, one can do

full_equations(sys)

\begin{align} \mathrm{\frac{d}{d t}}\left( mass{+}pos(t)1 \right) =& mass{+}v(t)1 \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)1 \right) =& - \frac{spring{+}k \left( - spring{+}l + \sqrt{\left|mass{+}pos(t)1\right|^{2} + \left|mass{+}pos(t)2\right|^{2}} \right) mass{+}pos(t)1}{mass{+}m \sqrt{\left|mass{+}pos(t)1\right|^{2} + \left|mass{+}pos(t)2\right|^{2}}} \ \mathrm{\frac{d}{d t}}\left( mass{+}pos(t)2 \right) =& mass{+}v(t)2 \ \mathrm{\frac{d}{d t}}\left( mass{+}v(t)2 \right) =& -9.81 + \frac{ - spring{+}k \left( - spring{+}l + \sqrt{\left|mass{+}pos(t)1\right|^{2} + \left|mass{+}pos(t)2\right|^{2}} \right) mass{+}pos(t)2}{mass{+}m \sqrt{\left|mass{+}pos(t)1\right|^{2} + \left|mass{+}pos(t)2\right|^{2}}} \end{align}

This is done as follows:

prob = ODEProblem(sys, [], (0., 3.))
sol = solve(prob, Rosenbrock23())
plot(sol)