# Can ContinuousCallbacks in DifferentialEquations.jl be used to watch parameters?

I’m trying to use a `VectorContinuousCallback` to handle friction and end stops. Some of the conditions are not based on the state variables `u` but on a value stuffed into a field of the parameters `p`.

Code:

``````using DifferentialEquations
using Plots
using Parameters
using StaticArrays

@enum State LEFT_STOP SLIDING_RIGHT STOPPED SLIDING_LEFT RIGHT_STOP

@with_kw mutable struct Params
M::Float64 = 1.0 # mass of sliding object in kg
ga::Float64 = 9.81 # gravitational constant, m/s^2
μs::Float64 = 0.8 # coefficient of static friction
μk::Float64 = 0.5 # coefficient of dynammic friction
cor::Float64 = 0.2 # coefficient of restitution, 0=inelastic, 1=elastic
xls::Float64 = -0.5 # location of left stop in m
xrs::Float64 = 0.5 # location of right stop in m
state::State = STOPPED # state of system
netforce::Float64 = 0.0 # to be written by the ODE so callbacks can see force
end

function FricSMCondx(out, u, t, integ)
@unpack_Params integ.p
v,x = u
out[1] = (state==LEFT_STOP)*netforce
out[2] = (state==SLIDING_RIGHT)*(-v)
out[3] = (state==SLIDING_RIGHT)*(x-xrs)
out[4] = (state==STOPPED)*(netforce)
out[5] = (state==STOPPED)*(-netforce)
out[6] = (state==SLIDING_LEFT)*v
out[7] = (state==SLIDING_LEFT)*(-(x-xls))
out[8] = (state==RIGHT_STOP)*(-netforce)
println("condx ",t," ",out)
end
function FricSMaffect!(integ, idx)
println("affect ",idx)
if     idx==1 integ.p.state = SLIDING_RIGHT
elseif idx==2 integ.p.state = STOPPED
elseif idx==3 integ.p.state = RIGHT_STOP
elseif idx==4 integ.p.state = SLIDING_RIGHT
elseif idx==5 integ.p.state = SLIDING_LEFT
elseif idx==6 integ.p.state = STOPPED
elseif idx==7 integ.p.state = LEFT_STOP
elseif idx==8 integ.p.state = SLIDING_LEFT
end
end
cb = VectorContinuousCallback(FricSMCondx, FricSMaffect!, nothing, 8)

force_ext(t) = 1.0*t

function sbode(u,p,t)
v,x = u # unpacking state variables
@unpack_Params p # unpack parameters by name into function scope
if state==LEFT_STOP
p.netforce = force_ext(t)-M*ga*μs # friction to the left
dv = 0.0 # not moving
dx = 0.0 # not moving
elseif state==SLIDING_RIGHT
p.netforce = force_ext(t)-M*ga*μk # friction to the left
dv = M*p.netforce
dx = v
elseif state==STOPPED
fext = force_ext(t)
p.netforce = fext-M*ga*μs*sign(fext)
dv = 0.0
dx = 0.0
elseif state==SLIDING_LEFT
p.netforce = force_ext()+M*ga*μk # friction to the right
dv = M*p.netforce
dx = v
else # state==RIGHT_STOP
Ffriction = M*ga*μs # positive because to the right
p.netforce = force_ext(t)+Ffriction
dv = 0.0
dx = 0.0
end
println("ode ",t," ",netforce, " ",state)
@SVector[dv, dx]
end

function sim()
u0 = @SVector[0.0, 0.0]
tspan = (0.0, 100.0)
p = Params()
prob = ODEProblem(sbode, u0, tspan, p)
sol = solve(prob, callback = cb)
end

sol = sim()
``````

The idea is that when the system is stopped due to friction, the callback will trigger on the force becoming high enough to cause motion, then the state will be changed, changing the model behavior accordingly.

However, what I see from an excerpt of the `println` outputs is the following:

``````ode 0.2720999999999999 -7.736900000000001 STOPPED
ode 0.4380999999999998 -7.575900000000001 STOPPED
ode 1.0110999999999997 -7.409900000000001 STOPPED
ode 1.091125540904509 -6.836900000000001 STOPPED
ode 1.1110999999999995 -6.756874459095492 STOPPED
ode 1.1110999999999995 -6.736900000000001 STOPPED
condx 0.11109999999999996 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 1.1110999999999995 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.22221111111111103 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.3333222222222221 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.4444333333333331 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.5555444444444442 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.6666555555555553 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.7777666666666664 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.8888777777777774 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 0.9999888888888885 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
condx 1.1110999999999995 [-0.0, -0.0, -0.0, -6.736900000000001, 6.736900000000001, 0.0, -0.0, 0.0]
ode 2.7210999999999985 -6.736900000000001 STOPPED
ode 4.381099999999998 -5.126900000000003 STOPPED
ode 10.111099999999995 -3.4669000000000025 STOPPED
ode 10.911355409045091 2.2630999999999943 STOPPED
ode 11.111099999999993 3.0633554090450907 STOPPED
ode 11.111099999999993 3.2630999999999926 STOPPED
condx 1.1110999999999995 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 11.111099999999993 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 2.22221111111111 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 3.3333222222222205 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 4.444433333333331 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 5.5555444444444415 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 6.666655555555551 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 7.777766666666662 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 8.888877777777772 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 9.999988888888883 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
condx 11.111099999999993 [0.0, -0.0, -0.0, 3.2630999999999926, -3.2630999999999926, 0.0, -0.0, -0.0]
``````

In this test case, `out[4]` of the conditional crosses zero in a rising direction between 1.11 s and 11.1 s (see the final iterations of the ode at those times) but when the conditional is tested at different time points it has fixed values of the `p.netforce` field and the affect function is never fired. Is there a way to structure this so the conditional can see (and interpolate on to find the zero cross) the parameter value?

You need to be evaluating the force function in the condition in order for that to work. It cannot necessarily use the evaluations of `f` directly since those are at different times from the zero function’s bisection.

Thanks. I’ll do some refactoring.

I think it’s working now - I had to move that force calculation like you said and also set a `dtmax` because in the stopped states there’s essentially no ODE left and it gets too aggressive about skipping forward in time. Here you can see as the sinusoidal force grows big enough to overcome the friction it starts to move, stopping in a couple of places on the first cycle, then with the coefficient of restitution there’s some bouncing occurring when it hits the end stops. I need to do some code cleanup and some sanity checking but I’m happy to get this working.

Thanks again.

Cool. That’s a pretty wild model.

That is neat, I like this idea. I had been trying to do this with conditionals in the model but it looks like events work better.

I’ve visited versions of this problem a number of times over the years in embedded systems - this could be a printer carriage or a proportional valve or a pneumatic device. This is a primitive friction model (there are many levels of friction models out there) but as they get more complicated it’s more work to characterize your system to the fancier model and this is often enough to reveal control system weaknesses.

I tried to do it a while back and didn’t get it to work to my satisfaction. I have a real world problem that could use this and thought up this structure this morning. I’m going to write this up as a blog post in the next few days after I clean a few things up and add an illustration or two.

Here’s a somewhat cleaned up version: Using state machines and callbacks to model friction and limit stops.

3 Likes

EDIT: My comments are relevant to viscous friction, which upon reading further I realize is not the case. The kinetic friction model is just a constant force, not one decaying with speed, so everything below doesn’t apply.

This is very cool, and contains a lot of functionality in very compact code. I was surprised by the absence of simulation parameters, particularly thresholds. A threshold could decide when velocity is close enough to zero to trigger a change in friction state. In principle, it would take forever for velocity to exponentially decay to exactly zero under zero force. Upon closer inspection, I suspect there is indeed an implicit threshold, perhaps from Float64’s limited representation. Plus the non-zero forcing function plays a big role. But I would expect that the simulation could take a long time to actually come to rest under zero force.

I would therefore suggest including an explicit threshold for “zero” velocity. First, real-world friction probably causes a sliding mass to stop at speeds higher than smallest Float64 (is that called ULP?). Second, it would make the simulation more independent of the particulars of computation. If someone switched to Float128, they might be surprised to find the mass taking longer to slow to rest, and worse for arbitrary precision arithmetic. (It looks like a threshold is not needed for initiating motion, because static friction is probably a bigger effect, and the forcing function changes quickly.)

I don’t intend it to be the best or complete friction model, in fact it’s pretty crude. There’s are a lot of friction models - many try to avoid discontinuities by using `atan` or similar shaped functions, and that seems to be as much about finding a model that can be simulated as about finding an accurate model. The callback system means we don’t have to force fit something smooth if it doesn’t make sense.