DifferentialEquations; how to simulate limit stops

Specific Domains Modelling & Simulations
1

I’m trying to simulate a spring mass system in which the mass is constrained by limits stops or walls.

There’s a driving force which, if it pushes the mass to one wall or the other, essentially halts the system until the net force will take the mass off the wall. For MWE purposes, imagine a sinusoidal source of force which is of large enough magnitude to alternately pin the mass against one wall or the other, and the interesting behavior is what happens in between.

I’m having trouble figuring out how to arrange callbacks to achieve this, partly because the condition to get off a wall is not dependent on the state (which in the MWE would only include position and velocity, but not the force acting on the mass).

Any thoughts about how to structure such a problem?

2

Adding a MWE:

using DifferentialEquations, Plots
plotly()

fdrive(t) = t%4.0>2.0 ? 5.0 : -5.0

function spring_mass_ODE!(du,u,p,t)
	ẋ, x = u
	mass, k, kd = p
	du[1] = 1/mass*(-k*x-kd*ẋ+fdrive(t))
	du[2] = ẋ
end

function smsim()
	tspan = [0.0, 10.0]
	u0 = [1.0, 1.0]
	p=[1.0, 1.0, 0.2]
	prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
	sol = solve(prob,Tsit5(), abstol=1e-7,reltol=1e-7)
end

function plotsmsim(sol)
	p1=plot(sol, vars=[1,],ylabel="velocity")
	p2=plot(sol, vars=[2,],ylabel="position")
	plot(p1,p2,layout=(2,1),link=:x,leg=false)
end

If you run smsim() |> plotsmsim you see a position that ranges from -6.something to over 6.something. I want to place walls at say, +/- 3, so that it runs against the wall at -3, then once the force reverses swings up to +3, sticks there a while, and so on.

3

What do you mean? So have a callback that starts the stick, and have a callback that unsticks using? A timer? You just need to specify what that model is.

4

One easy way that doesn’t need callback is to put a very stiff spring at each of your limits:

# useful for defining piecewise functions, 0 for negative argument, 1 for positive
heaviside(x) = 0.5 * (sign(x) + 1)                                               
                                                                                 
function spring_mass_ODE!(du,u,p,t)
    ẋ, x = u
    mass, k, kd, stop_rate = p
    negative_stop_force =  heaviside(-3 - x) * stop_rate
    positive_stop_force = -heaviside(x - 3) * stop_rate
    du[1] = 1/mass*(-k*x-kd*ẋ+negative_stop_force+positive_stop_force+fdrive(t))
    du[2] = ẋ
end

function smsim()
    tspan = [0.0, 10.0]
    u0 = [1.0, 1.0]
    p=[1.0, 1.0, 0.2, 100.0] # new fourth parameter for stop spring strength                                                      
    prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
    sol = solve(prob,Tsit5(), abstol=1e-7,reltol=1e-7)
end

Another solution that would need callbacks is the impulse method. A nice example is given in the documentation, using a limit detection function like

function detect_limit(u, t, integrator)                                          
    return (u[2] - 3)*(u[2] + 3)                                                 
end

Should let you do a two-sided version of that example.

1 Like
5

So here’s an attempt at use of callbacks to change the state of the model and it is almost working.

If you modify the MWE below by changing to MAGNITUDE = 10.0, you’ll see a response in which the limit conditions don’t get hit. If you set MAGNITUDE = 30.0, the limits seem to work, but it seems to get stuck on one of them. What I expect to see is several cycles in which the position saturates at the limits, but as the driving force decays, the end looks the same as when magnitude is started smaller. I suspect this is due to some kind of deadlock condition between the callbacks, but I’m not sure how to resolve it.

BTW, I am not trying to simulate elastic collisions as in the bouncing ball demo, but an inelastic one in which if the mass is pressed against one wall or the other it will sit there until the net force is reversed as appropriate to take it off the wall (I am not conserving momentum).

using DifferentialEquations, Plots
plotly()

@enum state neglimit moving poslimit

NEGLIMIT = -10.0
POSLIMIT = 10.0

MAGNITUDE = 30.0

fdrive(t) = MAGNITUDE*sin(2π*0.1*t)*exp(-t/15.0)
fspring(x) = -1.0*x
fdamping(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# if we are running into the negative travel limit
neglimit_cond(u,t,integrator) = NEGLIMIT-u[2]
function neglimit_affect!(integrator)
	integrator.p[1] = neglimit
	integrator.u[1] = 0.0
end
nlccb = ContinuousCallback(neglimit_cond, neglimit_affect!, nothing)

# if we are running into the positive travel limit
poslimit_cond(u,t,integrator) = u[2]-POSLIMIT
function poslimit_affect!(integrator)
	integrator.p[1] = poslimit
	integrator.u[1] = 0.0
end
plccb = ContinuousCallback(poslimit_cond, poslimit_affect!, nothing)

# if we are trying to get off of either travel limit condition
function moving_cond(u,t,integrator)
	ẋ,x = u
	state, = integrator.p
	state == moving   ? 0.0 : force(t,x,ẋ)
end
function moving_affect!(integrator)
	integrator.p[1] = moving
end
mvccb = ContinuousCallback(moving_cond, moving_affect!)
	
cbs = CallbackSet(nlccb,plccb,mvccb)
	
function spring_mass_ODE!(du,u,p,t)
	ẋ, x = u
	state, mass = p
	if state == moving
		du[1] = 1/mass*force(t,x,ẋ)
		du[2] = ẋ
	else
		du[1] = 0.0
		du[2] = 0.0
	end
end

function smsim()
	tspan = [0.0, 50.0]
	u0 = [0.0, 0.0]
	p=[moving, 1.0 ]
	prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
	sol = solve(prob, Tsit5(), callback=cbs, abstol=1e-7,reltol=1e-7)
end

function plotsmsim(sol)
	p1=plot(sol, vars=[1,],ylabel="velocity")
	p2=plot(sol, vars=[2,],ylabel="position")
	plot(p1,p2,layout=(2,1),link=:x,leg=false)
end
	
sol = smsim()
plotsmsim(sol)

MAGNITUDE = 10.0 response:
mag10

MAGNITUDE = 30.0 response:
mag30

6

I think the problem here, if I understand your system correctly, is that you need a different condition to leave each limit stop. The force must be positive to leave the negative stop and negative to leave the positive stop. moving_cond just waits for force to cross 0, which is not necessarily correct, and the reason it gets stuck at the negative limit is because the spring force dominates and the force stays positive.

I tried using a DiscreteCallback to fix this:

function start_moving(u, t, integrator)                                          
    ẋ, x = u                                                                     
    state, = integrator.p                                                        
    ((state == neglimit) && (force(t, x, ẋ) > 0) ||                              
     (state == poslimit) && (force(t, x, ẋ) < 0))                                
end                                                                              
mvdcb = DiscreteCallback(start_moving, moving_affect!)
cbs = CallbackSet(nlccb,plccb,mvdcb)

This works differently, it get stuck early but then unsticks, and I’m not sure why. I added force to the plot, it goes negative, but the systems remains stuck at the positive limit.

function plotsmsim(sol)
    p1=plot(sol, vars=[1,],ylabel="velocity")
    p2=plot(sol, vars=[2,],ylabel="position")
    tspan = sol.t[1]:0.01:sol.t[end]
    p3=plot(tspan, force.(tspan, sol(tspan, idxs=2), sol(tspan, idxs=1)), ylabel="force")
    plot(p1,p2,p3,layout=(3,1),link=:x,leg=false)
end

Setting dtmax=0.1 in the call to the solver does fix things, but this seems like a heavy hammer, I’m not sure how to give a gentler hint.

7

Thanks for the help. This is one of several things I’m working on right now and so I’m only able to get to it intermittently. I will do some troubleshooting by adding logging to the callbacks.

8

I now have two solutions but neither is quite right.

First, combining ContinuousCallbacks and a DiscreteCallback works but I have to set dtmax to a small value to get accuracy (is there an equivalent of rootfind=true for DiscreteCallback?)

using DifferentialEquations, Plots
plotly()

@enum state neglimit moving poslimit

NEGLIMIT = -10.0
POSLIMIT = 10.0

MAGNITUDE = 25.0

fdrive(t) = MAGNITUDE*sin(2π*0.1*t)*exp(-t/15.0)
fspring(x) = -1.0*x
fdamping(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# if we are running into the negative travel limit
neglimit_cond(u,t,integrator) = NEGLIMIT-u[2]
function neglimit_affect!(integrator)
    integrator.p[1] = neglimit
    integrator.p[2] = neglimit
    integrator.u[1] = 0.0
end
nlccb = ContinuousCallback(neglimit_cond, neglimit_affect!, nothing)

# if we are running into the positive travel limit
poslimit_cond(u,t,integrator) = u[2]-POSLIMIT
function poslimit_affect!(integrator)
    integrator.p[1] = poslimit
    integrator.p[2] = poslimit
    integrator.u[1] = 0.0
end
plccb = ContinuousCallback(poslimit_cond, poslimit_affect!, nothing)

# if we are trying to get off of neglimit travel limit condition
function moving_cond(u,t,integrator)
    ẋ,x = u
    state, callback_prev_state = integrator.p
    if callback_prev_state == neglimit
        out = force(t,x,ẋ)>0
    elseif callback_prev_state == poslimit
        out = force(t,x,ẋ)<0
    else
        out = false
    end
    out
end
function moving_affect!(integrator)
    integrator.p[1] = moving
end
mvccb = DiscreteCallback(moving_cond, moving_affect!)
    
cbs = CallbackSet(nlccb,plccb,mvccb)
    
function spring_mass_ODE!(du,u,p,t)
    ẋ, x = u
    state, calback_prev_force, mass = p
    if state == moving
        du[1] = 1/mass*force(t,x,ẋ)
        du[2] = ẋ
    else
        du[1] = 0.0
        du[2] = 0.0
    end
end

function smsim()
    tspan = [0.0, 50.0]
    u0 = [0.0, 0.0]
    p=[moving, moving, 1.0 ]
    prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
    sol = solve(prob, Tsit5(), callback=cbs, abstol=1e-7,reltol=1e-7,dtmax=0.1)
end

function plotsmsim(sol)
    p1=plot(sol, vars=[1,],ylabel="velocity")
    p2=plot(sol, vars=[2,],ylabel="position")
    plot(p1,p2,layout=(2,1),link=:x,leg=false)
end

function simplot()
    sol = smsim()
    plotsmsim(sol)
end

The second uses two DiscreteCallbacks but has an issue in which the solutions “buzzes” when it’s against the limit stops.

using DifferentialEquations, Plots
plotly()

NEGLIMIT = -10.0
POSLIMIT = 10.0

MAGNITUDE = 25.0

fdrive(t) = MAGNITUDE*sin(2π*0.1*t)*exp(-t/15.0)
fspring(x) = -1.0*x
fdamping(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# if we are running into the negative travel limit
neglimit_cond(u,t,integrator) = NEGLIMIT-u[2]≥0
function neglimit_affect!(integrator)
    ẋ, x = integrator.u
    integrator.u[1] = 0.0
    integrator.u[2] = NEGLIMIT
end
nlccb = DiscreteCallback(neglimit_cond, neglimit_affect!)

# if we are running into the positive travel limit
poslimit_cond(u,t,integrator) = u[2]-POSLIMIT≥0
function poslimit_affect!(integrator)
    ẋ, x = integrator.u
    integrator.u[1] = 0.0
    integrator.u[2] = POSLIMIT
end
plccb = DiscreteCallback(poslimit_cond, poslimit_affect!)

cbs = CallbackSet(nlccb,plccb)
    
function spring_mass_ODE!(du,u,p,t)
    ẋ, x = u
    mass, = p
    du[1] = 1/mass*force(t,x,ẋ)
    du[2] = ẋ
end

function smsim()
    tspan = [0.0, 50.0]
    u0 = [0.0, 0.0]
    p=[1.0]
    prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
    sol = solve(prob, Tsit5(), callback=cbs, abstol=1e-9,reltol=1e-9)
end

function plotsmsim(sol)
    p1=plot(sol, vars=[1,],ylabel="velocity")
    p2=plot(sol, vars=[2,],ylabel="position")
    plot(p1,p2,layout=(2,1),link=:x,leg=false)
end

function simplot()
    sol = smsim()
    plotsmsim(sol)
end
9

Attaching plots:
From the first simulation (with the buzzing):
sol1
and from the second:
sol2.

10

It looks like that should be directional ContinuousCallbacks?

11

I tried that but it didn’t catch all cases.

When there is a collision, there are two possible scenarios:

  • the driving force can still be in the direction of the wall and can hold it there until the driving force changes sign (the first three collisions visible in the plots), or
  • the driving force was already of opposite sign at time of collision (fourth collision in plot), so the mass coasts into the wall and should immediately bounce off, with exactly one point at the limit position.

I first wrote similar code using ContinuousCallbacks, but they would not find the second scenario as there was no sign change (at least after the last callback event?) to be found.

[EDIT - writing this reply gave me another idea for how to attack it. Will post again in a while.]