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

Specific Domains Modelling & Simulations
1

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?

2

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.

3

Thanks. I’ll do some refactoring.

4

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.

image

Thanks again.

5

Cool. That’s a pretty wild model.

6

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.

7

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.

8

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.

9

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

3 Likes
10

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

11

Hi, thanks for reading.

I have a couple of real world systems in which I really needed the end stop behavior and adding static and dynamic friction seemed like a good idea while I was at it. I’m hoping this is a decent example of how to use callbacks and state machines to handle what can be very confusing cases of varying system behavior.

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.