Simulation of a pendulum with dry friction using DifferentialEquations.jl and JuMP.jl

PendulumWithDryFriction.jl

using DifferentialEquations
using JuMP, PATHSolver
using Plots

g = 9.81 # gravitational acceleration
l = 1 # pendulum's length
m = 1 # pendulum's mass
τₛ = 1 # static friction torque

model = Model(PATHSolver.Optimizer)
@variable(model, -τₛ <= τ <= τₛ)
@variable(model, bias)
@constraint(model, - bias + τ ⟂ τ)
set_silent(model)

function dynamics!(du, u, p, t)
	du[1] = u[2]

	c = g * m * l * sin(u[1])

	if u[2] == 0
		fix(bias, c)
		optimize!(model)
		du[2] = (- c + value(τ)) / (m * l^2)
	else
		du[2] = (- c - τₛ * sign(u[2])) / (m * l^2)
	end

	return nothing
end

condition(u, t, integrator) = u[2]
affect!(integrator) = nothing
velocityCallback = ContinuousCallback(condition, affect!)

f! = ODEFunction(dynamics!, syms=[:φ, :ω])
prob = ODEProblem(f!, [pi/4, 0], [0., 5.])
sol = solve(prob, callback=velocityCallback)

plot(sol)