Can OrdinaryDiffEq.jl solve Fillippov system?

Fillippov systems are piecewise smooth ODEs and the solutions can slide along the switch surfaces. Here is an excellent illustration of such systems, taken from this page [Events and Discontinuities in Differential Equations—Wolfram Language Documentation]


Here is a typical Fillippov system:

\begin{aligned} \dot{x}=&y, \\ \dot{y}=&-cy-kx+F(y-v_{s})+\epsilon \cos(\omega t),\\ \end{aligned}

where F is given by

F(x)=-A\ \text{Sign}(x)(\frac{\alpha}{1+\gamma|x|}+\lambda+\eta x^{2}).

So I use OrdinaryDiffEq.jl to solve this ODE. Here is my code:

using StaticArrays, OrdinaryDiffEq, Plots

function F(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A*sign(x)*(α/(1+γ*abs(x))+λ+η*x^2)
end

function vec(x,p,t)
    SVector(x[2],-p[1]*x[2]-p[2]*x[1]+F(x[2]-p[3])+p[4]*cos(p[5]*t)) # p=[c,k,v_s,ϵ,ω]
end

prob=ODEProblem(vec, SVector(1.0,3.0),(0,20),SVector(0.2,1,1,0.2,1.067))

sol=solve(prob, Vern6())

plot(sol, idxs = (1, 2))

Here is the solution:


As you can see from the picture, in the switch surface y=1, the solution is oscillating slightly and the integrator spent lots of time on the switch surface. It’s clearly not what I want.

I have tried using the event handler, but it just gives the wrong solution:

using StaticArrays, OrdinaryDiffEq, Plots

function F(x,s)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A*s*(α/(1+γ*abs(x))+λ+η*x^2)
end

function vec(x,p,t)
    SVector(x[2],-0.2*x[2]-x[1]+F(x[2]-1, p[1])+0.2*cos(1.067*t)) 
end

function h(u, t, integrator)
    u[2]-integrator.p[1]
end

function affect!(integrator)
    integrator.p=-integrator.p
end

cb=ContinuousCallback(h, affect!)

prob=ODEProblem(vec, SVector(1.0,3.0),(0,20),[1])

sol=solve(prob, Vern6(), callback=cb, dtmax=0.1)

plot(sol, idxs = (1, 2))


So, does any solution available?

Those aren’t the same equation. In the first case you have -A*sign(x)* .... In the second case, you have -A *s*..., so presumably s = sign(x) where s = p[1]. So does p[1] = sign(x)? In the second case you have the condition u[2]-integrator.p[1], so when u[2] = integrator.p[1] you switch the parameters. When you start, p[1] = 1, so when u[2] = 1 you swap p[1] = -1. Now next time around, the sign change would occur when u[2] = -1. In your original model, the point where the sign changed did not move after the first crossing, so these are not the same model.

If I do:

function F(x)
    @show sign(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A*sign(x)*(α/(1+γ*abs(x))+λ+η*x^2)
end

you see:

sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0

Those oscillations are real. They are your system switching over and over and over. It seems that there are just “infinitely many” events going on in the small interval.

Thank you for pointing out why my second code was wrong! Yes, the solution was oscillating around u[2]=1, since the vector fields in the two sides of u[2]=1 are squeezing towards u[2]=1. Just like this picture:
image
In non-smooth nonlinear dynamics, this phenomenon is called the sliding mode. My question is whether the existing solvers can handle this sliding mode more accurately. For example, detect when the solution enters the sliding mode or when it exits, and the solution can just slide along the switch surface, without oscillating.

This seems more like something for a bifurcation analysis tool to me. Have you considered something like BifurcationKit?

https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/

1 Like

For this type of thing, you have to simulate differently at the switching surface to eliminate infinite switch. This series of lectures explain this https://www.youtube.com/watch?v=MNbc69pcio0&list=PLp8ijpvp8iCvFDYdcXqqYU5Ibl_aOqwjr&index=37

1 Like

OK. I will check that.

I am afraid not. The parameters of the system do not change during the simulation. I am just looking for a suitable solver for this type of system.

I only have a small amount of experience with these kinds of problems.

I would try to use x[1] as a parameter (instead of a dynamical variable) in the bifurcation analysis. It should then yield this ‘sliding mode’ as a fix point manifold.

Hi, Dear Chris @ChrisRackauckas. I am designing a function to switch the sliding mode manually. In the document of DifferentialEquation.jl, the function ContinuousCallback has a keyword repeat_nudge to avoid repeating find zeros. So, if I have an initial time t0, can I set a fixed time, say 1//100? Only after t0+1//100, the event handler starts to work.

You can change the parameter. It’s just a keyword argument. But this has so many nearby events that I would suggest @josuagrw 's approach as the issue is going to be down to floating point error

1 Like

Hello everyone, thank you for all your advice. Now I have a solution for this type of problem. I will just give the code, and refer to this paper for basic knowledge about the Fillippov systems.

using StaticArrays, OrdinaryDiffEq, Plots, LinearAlgebra

struct ODESol
    solution # ODE solution solved by OrdinaryDiffEq.jl
    state::Symbol # only two symbols are used: :regularmode and :slidemode
end

function fillippov_solve(f1, f2, h, ∇h, v0, p0, t0, tend)
    # it is assumed that f1 is at h(x)<0 and f2 is at h(x)>0
    # f1, f2, h, ∇h are all type of f(u, p, t), the output is a vector
    # initial condition: v0 
    # time interval: t0,tend
    # parameter p0
    # 
    # first define the vector field on the switch surface
    function sv(x, p, t)
        α = dot(∇h(x, p, t), f1(x, p, t)) / dot(∇h(x, p, t), f1(x, p, t) - f2(x, p, t))
        (1 - α) * f1(x, p, t) + α * f2(x, p, t)
    end

    # the hyper surface 
    function newh(x, t, integrator)
        h(x, integrator.p, t)
    end
    # stop the integrator
    function stop(integrator)
        if (integrator.t-integrator.sol.t[1])< 1e-2 # to avoid the integrator stopping at the beginning when it exits the sliding surface
            nothing
        else
            terminate!(integrator)
        end
    end

    # exit the sliding mode when slide_exit=0
    function slide_exit(x, t, integrator)
        dot(∇h(x, integrator.p, t), f2(x, integrator.p, t)) * dot(∇h(x, integrator.p, t), f1(x, integrator.p, t))
    end

    # main data consist of different solutions at h>0, h<0 and h=0
    sols = Vector{ODESol}(undef, 1)

    if abs(h(v0, p0, t0)) < 1e-10 # if this happen, we assume that the initial condition is on h=0
        prob = ODEProblem(sv, v0, (t0, tend), p0)
        cb = ContinuousCallback(slide_exit, stop)
        firstsol = ODESol(solve(prob, Vern6(), callback=cb), :slidemode)
        sols[1] = firstsol
    elseif h(v0, p0, t0) > 0
        prob = ODEProblem(f2, v0, (t0, tend), p0)
        cb = ContinuousCallback(newh, stop)
        firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
        sols[1] = firstsol
    else
        prob = ODEProblem(f1, v0, (t0, tend), p0)
        cb = ContinuousCallback(newh, stop)
        firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
        sols[1] = firstsol
    end
    while sols[end].solution.t[end] != tend # if the last of sols's time ≠ tend
        v1 = sols[end].solution.u[end]
        t1 = sols[end].solution.t[end]
        if sols[end].state == :slidemode # if the last solution is in slidemode, we have to determine which vector field it enter
            if abs(dot(∇h(v1, p0, t1), f2(v1, p0, t1))) < 1e-8 && dot(∇h(v1, p0, t1), f1(v1, p0, t1)) > 1e-5
                prob = ODEProblem(f2, v1, (t1, tend), p0)
                cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                append!(sols, [firstsol])
            elseif abs(dot(∇h(v1, p0, t1), f1(v1, p0, t1))) < 1e-8 && dot(∇h(v1, p0, t1), f2(v1, p0, t1)) < -1e-5
                prob = ODEProblem(f1, v1, (t1, tend), p0)
                cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                append!(sols, [firstsol])
            else
                break
                print("The solution is not defined on state $v1 and time $t1")
            end
        else
            if dot(∇h(v1, p0, t1), f1(v1, p0, t1)) * dot(∇h(v1, p0, t1), f2(v1, p0, t1)) > 0 # crossing h=0
                if dot(∇h(v1, p0, t1), f1(v1, p0, t1)) > 0
                    prob = ODEProblem(f2, v1, (t1, tend), p0)
                    cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                    firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                    append!(sols, [firstsol])
                else
                    prob = ODEProblem(f1, v1, (t1, tend), p0)
                    cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                    firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                    append!(sols, [firstsol])
                end
            else # entering h=0
                prob = ODEProblem(sv, v1, (t1, tend), p0)
                cb = ContinuousCallback(slide_exit, stop)
                firstsol = ODESol(solve(prob, Vern6(), callback=cb), :slidemode)
                append!(sols, [firstsol])
            end
        end
    end
    sols
end

The main function is fillippov_solve, and its output is a vector consist of ODESol.
Now we can test this function using the system I gave before. First, define two vector fields:

function F1(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    A * (α / (1 + γ * abs(x)) + λ + η * x^2)
end

function F2(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A * (α / (1 + γ * abs(x)) + λ + η * x^2)
end

function vec1(x, p, t)
    SVector(x[2], -p[1] * x[2] - p[2] * x[1] + F1(x[2] - p[3]) + p[4] * cos(p[5] * t)) # p=[c,k,vs,ϵ,ω]
end

function vec2(x, p, t)
    SVector(x[2], -p[1] * x[2] - p[2] * x[1] + F2(x[2] - p[3]) + p[4] * cos(p[5] * t)) # p=[c,k,vs,ϵ,ω]
end

And then define the hypersurface and its grad:

function H(x, p, t)
    x[2] - p[3]
end

function ∇H(x, p, t)
    SVector(0.0, 1.0)
end

Note that we have to make sure that when H>0 the vector field is vec2 and when H<0 the vector field is vec1.
Now we can use fillippov_solve to solve this ODE. Just run

sol = fillippov_solve(vec1, vec2, H, ∇H, SVector(1.0, 3.0), SVector(0.2, 1, 1, 0.2, 1.067), 0.0, 20.0)

You will get a vector.

sol[2].state

gives

:slidemode
sol[2].solution

gives a ODE solution

retcode: Terminated
Interpolation: specialized 
6th order lazy interpolation
t: 3-element Vector{Float64}:
 0.5146478247426896
 2.117126414261962
 2.117126414261962
u: 3-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [2.07049529563982, 1.0000000000000895]
 [3.6729738851570852, 1.000000000000007]
 [3.6729738851570852, 1.000000000000007]

Since we define a new vector field sv on the switch surface when the vector fields are squeezing towards the switch surface, the solution is now stable.
We can define a simple function to plot all the solutions:

function myplot(sols)
    plotlyjs()
    figure = plot()
    for i in sols
        plot!(i.solution, idxs=(1,2))
    end
    figure
end

Now we run

myplot(sol)

it gives this picture:


Since I also consider the case of crossing, the codes should work on the usual piecewise smooth case.

2 Likes