🟡 Period doubling in Lur'e problem (PD aBS)

• 🟡 Period doubling in Lur'e problem (PD aBS)
• • Periodic orbits with orthogonal collocation
  • Periodic orbits with Parallel Standard Shooting
  • Branch of periodic orbits with finite differences

The following model is an adaptive control system of Lur’e type. It is an example from the MatCont library.

\[\left\{\begin{array}{l} \dot{x}=y \\ \dot{y}=z \\ \dot{z}=-\alpha z-\beta y-x+x^{2} \end{array}\right.\]

The model is interesting because there is a period doubling bifurcation and we want to show the branch switching capabilities of BifurcationKit.jl in this case. We provide 3 different ways to compute this periodic orbits and highlight their pro / cons.

It is easy to encode the ODE as follows

using Revise, Plots
using BifurcationKit
const BK = BifurcationKit

recordFromSolution(x, p; k...) = (u1 = x[1], u2 = x[2])

function lur!(dz, u, p, t = 0)
	(;α, β) = p
	x, y, z = u
	dz[1] = y
	dz[2] =	z
	dz[3] = -α * z - β * y - x + x^2
	dz
end

# bifurcation problem
prob_bif = ODEBifProblem(lur!, zeros(3), (α = -1.0, β = 1.0), (@optic _.α);
    record_from_solution = recordFromSolution)

We first compute the branch of equilibria

# continuation options
opts_br = ContinuationPar(p_min = -1.4, p_max = 1.8, dsmax = 0.01, max_steps = 1000)

# computation of the branch
br = continuation(prob_bif, PALC(), opts_br)

scene = plot(br)

With detailed information:

br

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 200
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter α starts at -1.0, ends at 1.8
 ├─ Algo: PALC
 └─ Special points:

- #  1,     hopf at α ≈ +1.00022831 ∈ (+0.99934442, +1.00022831), |δp|=9e-04, [converged], δ = (-2, -2), step = 142
- #  2, endpoint at α ≈ +1.80000000,                                                                     step = 199

We note the Hopf bifurcation point which we shall investigate now.

Periodic orbits with orthogonal collocation

We compute the branch of periodic orbits from the Hopf bifurcation point. We rely on a the state of the art method for computing periodic orbits of ODE: orthogonal collocation.

We first define a plotting function and a record function which are used for all cases below:

# plotting function
function plotPO(x, p; k...)
	xtt = get_periodic_orbit(p.prob, x, p.p)
	plot!(xtt.t, xtt[1,:]; markersize = 2, k...)
	plot!(xtt.t, xtt[2,:]; k...)
	plot!(xtt.t, xtt[3,:]; legend = false, k...)
end

# record function
function recordPO(x, p; k...)
	xtt = get_periodic_orbit(p.prob, x, p.p)
	period = getperiod(p.prob, x, p.p)
	mn, mx = extrema(xtt[1,:])
	return (;max = mx, min = mn, period)
end

recordPO (generic function with 1 method)

Continuation of periodic orbits from the Hopf point:

# continuation parameters
opts_po_cont = ContinuationPar(opts_br, dsmax = 0.03, ds = 0.01, dsmin = 1e-4, max_steps = 80, tol_stability = 1e-4, plot_every_step = 20)

br_po = continuation(
	br, 1, opts_po_cont,
	PeriodicOrbitOCollProblem(40, 4);
	plot = true,
	record_from_solution = recordPO,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; k...)
		# plot previous branch
		plot!(br, subplot = 1, putbifptlegend = false)
		end,
	normC = norminf)

scene = plot(br, br_po)

Note that you can compute the PD normal form

get_normal_form(br_po, 1)

Period-Doubling bifurcation point of periodic orbit
├─ Period = 6.364079037215513 -> 12.728158074431025
├─ Problem: PeriodicOrbitOCollProblem
├─ α ≈ 0.6302627694057382
├─ type: ├─ Normal form (Iooss):
├	∂τ = 1 + a₀₁⋅δp + a₂⋅ξ²
├	∂ξ = ξ⋅(c₁₁⋅δp + c₃⋅ξ²)
├─── a₀₁ = 0.030826838792648283
├─── a₂  = 0.003295291803080316
├─── c₁₁ = -1.2958259311950344
└─── c₃  = -0.30503757864675174

We provide Automatic Branch Switching from the PD point and computing the bifurcated branch is as simple as:

# aBS from PD
br_po_pd = continuation(deepcopy(br_po), 1, ContinuationPar(br_po.contparams, max_steps = 100, dsmax = 0.02, plot_every_step = 10, ds = 0.005);
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; k...)
		## add previous branch
		plot!(br_po; subplot = 1)
	end,
	record_from_solution = recordPO,
	normC = norminf,
	)

scene = plot(br_po, br_po_pd, title = "Collocation based")

Periodic orbits with Parallel Standard Shooting

We use a different method to compute periodic orbits: we rely on a fixed point of the flow. To compute the flow, we use DifferentialEquations.jl. This way of computing periodic orbits should be more precise than the previous one. We use a particular instance called multiple shooting which is computed in parallel. This is an additional advantage compared to the previous method. Finally, please note the close similarity to the code of the previous part. As before, we first rely on Hopf aBS.

import OrdinaryDiffEq as ODE

# ODE problem for using DifferentialEquations
prob_ode = ODE.ODEProblem(lur!, prob_bif.u0, (0., 1.), prob_bif.params; abstol = 1e-12, reltol = 1e-10)

# continuation parameters
# we decrease a bit the newton tolerance to help automatic branch switching from PD point
opts_po_cont = ContinuationPar(opts_br, dsmax = 0.03, newton_options = NewtonPar(tol = 1e-10), tol_stability = 1e-5, n_inversion = 8, nev = 3, max_steps = 100)

br_po = continuation(
	br, 1, opts_po_cont,
	# parallel shooting functional with 10 sections
	ShootingProblem(10, prob_ode, ODE.Vern9(); parallel = true);
	# plot = true,
	record_from_solution = recordPO,
	plot_solution = plotPO,
	# limit the residual, useful to help DifferentialEquations
	callback_newton = BK.cbMaxNorm(10),
	normC = norminf)

plot(br_po)

Note that you can compute the PD normal form

get_normal_form(br_po, 1; autodiff = true)

Period-Doubling bifurcation point of periodic orbit
├─ Period = 6.364071562915145 -> 12.72814312583029
├─ Problem: ShootingProblem
SuperCritical - Period-Doubling ┌─ Normal form:
├	 x ─▶ x⋅(a⋅δp - 1 + c⋅x²)
├─ a = -698.5786918693161
└─ c = 25.662290190964843

We provide Automatic Branch Switching from the PD point and computing the bifurcated branch is as simple as:

# aBS from PD
br_po_pd = continuation(deepcopy(br_po), 1,
	ContinuationPar(br_po.contparams, max_steps = 6, ds = -0.01);
	plot = true, verbosity = 2,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; k...)
		# add previous branch
		plot!(br_po; subplot = 1)
	end,
	record_from_solution = recordPO,
	autodiff_nf = true,
	normC = norminf,
	callback_newton = BK.cbMaxNorm(10),
	)

scene = plot(br, br_po, br_po_pd, title = "Shooting based")

Branch of periodic orbits with finite differences

We use finite differences to discretize the problem for finding periodic orbits. We appeal to automatic branch switching from the Hopf point as follows

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.02, dsmin = 1e-4, p_max = 1.1, max_steps = 80, tol_stability = 1e-4)

br_po = continuation(
	br, 1, opts_po_cont,
	PeriodicOrbitTrapProblem(M = 120);
	record_from_solution = recordPO,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; k...)
		## plot previous branch
		plot!(br, subplot=1, putbifptlegend = false)
		end,
	normC = norminf)

scene = plot(br, br_po)

Two period doubling bifurcations were detected. We shall now compute the branch of periodic orbits from these PD points. We do not provide Automatic Branch Switching as we do not have the PD normal form computed for PeriodicOrbitTrapProblem. Hence, it takes some trial and error to find the ampfactor of the PD branch.

# aBS from PD
br_po_pd = continuation(deepcopy(br_po), 1, ContinuationPar(br_po.contparams, max_steps = 70);
	# plot = true,
	ampfactor = .2, δp = -0.005,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; k...)
		# add previous branch
		plot!(br_po; legend=false, subplot=1)
	end,
	record_from_solution = recordPO,
	normC = norminf
	)
Scene = title!("")

plot(br, br_po, br_po_pd, title = "Trapezoid based")