[Dev] BK sprint

[jbcaillau]

@rveltz @j-l-s @ocots @gergaud notes to prepare BK sprint (Dec. 15 to 19) and collect input from OptimalControl.jl

• [ ] @gergaud @ocots can you please put a link towards your preliminary tests somewhere, using variable in the flow?
• FTR: use case to be tested with expected rich behaviour https://github.com/control-toolbox/OptimalControl.jl/issues/224 (check Fig. 5 below from this paper - typical swallowtail sing. for the value function)
• FTR: @abavoil 's post on discourse regarding ForwardDiff + ODE

See also: https://github.com/control-toolbox/OptimalControl.jl/issues/655

[ocots]

Discrete homotopy versus BifurcationKit homotopy:

using Pkg
Pkg.activate(".")

using OptimalControl
using Plots
using OrdinaryDiffEq
using NonlinearSolve

# parameters
t0 = 0
tf = 2
x0 = 0
t1 = 1
xf = 1 - exp(t1-tf)

# model
ocp = @def begin
    ε ∈ R, variable
    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control
    x(t0) == x0
    x(tf) == xf
    ẋ(t) == -x(t)+u(t)
    ∫( abs(u(t)) - ε*( log(abs(u(t))) + log( max( 1 - abs(u(t)), 1e-8 ) ) ) ) → min
end

# maximising control
function control(x, p, ε)    
    ψ = -1 + abs(p)
    u = -(2*ε)/(ψ - 2*ε - √(ψ^2+4*ε^2))
    return p < 0 ? -u : u
end

# Hamiltonian flow
f = Flow(ocp, control);

# state projection, p being the costate
π((x, p)) = x

# shooting function
S(p0, ε) = π( f(t0, x0, p0, tf, ε) ) - xf

# discrete homotopy on ε
ε_init  = 1
ε_final = 1e-4
ε_span = logrange(ε_init, ε_final, 10)
p0_span = Float64[]
let 
    p0_guess = 0.2
    for ε ∈ ε_span
        prob = NonlinearProblem(S, p0_guess, ε)
        sol = solve(prob)
        if !(sol.retcode === SciMLBase.ReturnCode.Success)
            error("Newton solver failed")
        end
        p0_guess =  sol.u[1]
        push!(p0_span, p0_guess) 
    end
end

# plot of the solutions
plt_sol = plot(; dpi=200)
for (p0, ε) ∈ zip(p0_span, ε_span)
    sol = f((t0, tf), x0, p0, ε; saveat=range(t0, tf, 200))
    plot!(plt_sol, sol; label="ε = $(ε)", legend=:none, color=1, linewidth=1)
end
ylims!(plt_sol[3], (0, 1))
plt_sol

# plot the homotopy path
plt_path = plot(p0_span, ε_span; label="Discrete homotopy path", xlabel="p₀", ylabel="ε", color=1, linewidth=2, dpi=200)
hline!(plt_path, [ε_init]; color=2, label=:none)
hline!(plt_path, [0]; color=2, label=:none)

# homotopy with BifurcationKit
import BifurcationKit: BifurcationProblem, continuation, PALC, ContinuationPar

F(p0::AbstractVector, λ::AbstractVector) = [S(p0[1], 1-λ[1])]

p0_init = p0_span[1]
λ_init  = 1 - ε_span[1]
λ_final = 1 - ε_span[end]

prob = BifurcationProblem(F, [p0_init], [λ_init], 1; record_from_solution = (x, p; k...) -> x[1])
Δε = abs(ε_final - ε_init)
dsmax = Δε/100
dsmin = Δε/1000
ds_init = dsmax 
options = ContinuationPar(p_min=λ_init, p_max=λ_final, dsmax=dsmax, dsmin=dsmin, ds=ds_init, max_steps=500)
br = continuation(prob, PALC(), options)

p0_span_bk = Float64[]
ε_span_bk  = Float64[]
for sol ∈ br.sol
    push!(p0_span_bk, sol.x[1])
    push!(ε_span_bk, 1 - sol.p[1])
end

plot!(plt_path, p0_span_bk, ε_span_bk; label="BifurcationKit", color=3, linewidth=2)

savefig(plt_path, "plot_path_homotopy.png")
savefig(plt_sol, "plot_sol_homotopy.png")

[ocots]

Optimal control problem with a bifurcation.

using Pkg
Pkg.activate(".")

using OptimalControl
using Plots
using OrdinaryDiffEq
import BifurcationKit: BifurcationProblem, continuation, PALC, ContinuationPar

# parameters
t0 =  0
x0 = -2
tf =  1

# model
ocp = @def begin
    λ ∈ R, variable
    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control
    x(t0) == x0
    x(tf)^2 - λ^2 == 0
    ẋ(t) == u(t)
    ∫( u(t)^2 ) → min
end

# maximising control
u(x, p, λ) = p/2

# Hamiltonian flow
f = Flow(ocp, u);

# state projection
π((x, p)) = x

# shooting function
S(p0::Number, λ::Number) =  π( f(t0, x0, p0, tf, λ) )^2 - λ^2
S(p0::AbstractVector, λ::AbstractVector) = [S(p0[1], λ[1])]

# continuation: branch 1
p0_init =  2
λ_init  = -1
λ_final =  1
prob = BifurcationProblem(S, Float64[p0_init], Float64[λ_init], 1; 
    record_from_solution = (x,p; k...) -> x[1])
br = continuation(prob, PALC(), ContinuationPar(p_min = Float64(λ_init), p_max = Float64(λ_final)))
plt_branch_1 = plot(br; dpi=200)

savefig(plt_branch_1, "plot_branch_1.png")

# continuation: branch 2
p0_init =  6
λ_init  = -1
λ_final =  1
prob = BifurcationProblem(S, Float64[p0_init], Float64[λ_init], 1; 
    record_from_solution = (x,p; k...) -> x[1])
br = continuation(prob, PALC(), ContinuationPar(p_min = Float64(λ_init), p_max = Float64(λ_final)))
plt_branch_2 = plot(br; dpi=200)

savefig(plt_branch_2, "plot_branch_2.png")

Bifurcation Analysis of Parameterized Optimal Control Problem

1. Problem Formulation

We consider the following parameterized optimal control problem:

    (P)\left\{
    \begin{array}{l}
        \min \int_{0}^{1}u^2\,dt \\[0.5em]
        \dot{x}=u \\[0.3em]
        x(0)=-2 \\[0.3em]
        x(1)^2=\lambda^2
    \end{array}
    \right.

where:

• $x \in \mathbb{R}$ is the state variable
• $u \in \mathbb{R}$ is the control variable
• $\lambda \in \mathbb{R}$ is a parameter
• The initial condition is fixed at $x(0) = -2$
• The final condition is parameterized: the state must reach $\pm\lambda$ at time $t=1$

This is a minimum energy problem where we seek to transfer the system from $x(0)=-2$ to a final state satisfying $x(1)^2 = \lambda^2$, minimizing the control effort.

2. Pontryagin Maximum Principle (PMP)

2.1 Maximizing Control (Feedback Form)

The Hamiltonian of the system is:

    H(t,x,p,u) = -u^2 + pu

where $p$ is the costate variable. Applying the maximization condition from PMP:

    \frac{\partial H}{\partial u} = -2u + p = 0

yields the optimal control in feedback form:

    u(x,p) = \frac{p}{2}

2.2 True Hamiltonian

Substituting the optimal control into the Hamiltonian gives the true Hamiltonian:

    H_r(t,x,p) = H(t,x,p,u(x,p)) = -\frac{p^2}{4} + p \cdot \frac{p}{2} = \frac{p^2}{4}

2.3 Hamiltonian Dynamics

The state and costate dynamics are governed by:

    \left\{
    \begin{array}{l}
        \dot{x} = \frac{\partial H_r}{\partial p} = \frac{p}{2} \\[0.5em]
        \dot{p} = -\frac{\partial H_r}{\partial x} = 0
    \end{array}
    \right.

3. Extremal Solutions

3.1 General Solution

Since $\dot{p} = 0$, the costate is constant: $p(t) = a \in \mathbb{R}$ for all $t \in [0,1]$.

Integrating the state equation with initial condition $x(0) = -2$:

    \left\{
    \begin{array}{ll}
        x(t) = \frac{at}{2} - 2 & t\in[0,1]\\[0.5em]
        p(t) = a  & t\in[0,1]
    \end{array}
    \right.

3.2 Determination of the Costate Parameter

The final condition $x(1)^2 = \lambda^2$ determines the value of $a$ as a function of $\lambda$:

    \begin{array}{rcl}
        x(1)^2 = \lambda^2  & \Leftrightarrow &  \left(\frac{a}{2} - 2\right)^2 = \lambda^2 \\[0.5em]
        & \Leftrightarrow & \frac{a^2}{4} - 2a + 4-\lambda^2 = 0 \\[0.5em]
        & \Leftrightarrow & a^2 - 8a + 16 - 4\lambda^2 = 0
    \end{array}

Solving this quadratic equation:

    a = \frac{8 \pm \sqrt{64 - 64 + 16\lambda^2}}{2} = \frac{8 \pm 4\lambda}{2} = 4 \pm 2\lambda

This yields two branches:

    \begin{array}{l}
        a^1_\lambda = 4 - 2\lambda \\[0.3em]
        a^2_\lambda = 4 + 2\lambda
    \end{array}

4. Parametrization of the Solution Branches

For $i \in {1, 2}$, we define the initial costate as a function of the parameter:

    \begin{array}{lllll}
        p^i_0 & : & [-1,1]  & \longrightarrow & \mathbb{R} \\  
            &   & \lambda & \longmapsto     & p^i_0(\lambda) = a^i_\lambda
    \end{array}

We represent each branch as a curve in the $(p_0, \lambda)$ plane:

    \begin{array}{lllll}
        c^i & : & [-1,1]  & \longrightarrow & \mathbb{R}^2 \\  
            &   & \lambda & \longmapsto     & c^i(\lambda) = (p^i_0(\lambda), \lambda)
    \end{array}

Explicitly:

• Branch 1: $c^1(\lambda) = (4 - 2\lambda, \lambda)$
• Branch 2: $c^2(\lambda) = (4 + 2\lambda, \lambda)$

5. Arc-Length Parametrization

5.1 Computation of Arc Length

To perform bifurcation analysis, we reparametrize using arc length $s$. Computing the derivatives:

    \frac{dc^1}{d\lambda} = (-2, 1), \quad \frac{dc^2}{d\lambda} = (2, 1)

Both have the same norm: $\left|\frac{dc^i}{d\lambda}\right| = \sqrt{4 + 1} = \sqrt{5}$

The arc length from $\lambda = -1$ to $\lambda$ is:

    s(\lambda) = \int_{-1}^{\lambda} \left\|\frac{dc}{d\mu}\right\| d\mu = \int_{-1}^{\lambda} \sqrt{5}\, d\mu = \sqrt{5}(\lambda+1)

The total arc length for $\lambda \in [-1, 1]$ is $s_{\text{max}} = 2\sqrt{5}$.

5.2 Tangent Vectors in Arc-Length Parametrization

The tangent vectors with respect to arc length are:

    \begin{array}{lllll}
        \dot{c}^i & : & [0, 2\sqrt{5}]  & \longrightarrow & \mathbb{R}^2 \\  
            &   & s & \longmapsto     & \dot{c}^i(s) = \frac{1}{\sqrt{5}}\frac{dc^i}{d\lambda}
    \end{array}

Explicitly:

• Branch 1: $\dot{c}^1(s) = \left(-\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right)$
• Branch 2: $\dot{c}^2(s) = \left(\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right)$

5.3 Curves in Arc-Length Parametrization

Converting to arc-length parameter where $s = \sqrt{5}(\lambda + 1)$, i.e., $\lambda = \frac{s}{\sqrt{5}} - 1$:

    \begin{array}{lllll}
        c^i & : & [0, 2\sqrt{5}]  & \longrightarrow & \mathbb{R}^2 \\  
            &   & s & \longmapsto     & c^i(s) = \left((-1)^i\frac{2s}{\sqrt{5}} + a^i_{-1}, \frac{s}{\sqrt{5}} - 1\right)
    \end{array}

where $a^1_{-1} = p^1_0(-1) = 4 - 2(-1) = 6$ and $a^2_{-1} = p^2_0(-1) = 4 + 2(-1) = 2$.

Branch 1:

    c^1(s) = \left(-\frac{2s}{\sqrt{5}} + 6, \frac{s}{\sqrt{5}} - 1\right)

Branch 2:

    c^2(s) = \left(\frac{2s}{\sqrt{5}} + 2, \frac{s}{\sqrt{5}} - 1\right)

6. Bifurcation Analysis

6.1 Bifurcation Point

The two branches meet when $c^1(s) = c^2(s)$:

From the $\lambda$-component: both give $\lambda = \frac{s}{\sqrt{5}} - 1$, so we need:

    -\frac{2s}{\sqrt{5}} + 6 = \frac{2s}{\sqrt{5}} + 2

Solving:

    4 = \frac{4s}{\sqrt{5}} \Rightarrow s = \sqrt{5}

At $s = \sqrt{5}$ (corresponding to $\lambda = 0$), the bifurcation point is:

    c^1(\sqrt{5}) = c^2(\sqrt{5}) = (4, 0)

This means:

• $p_0 = 4$ (initial costate)
• $\lambda = 0$ (parameter value)

6.2 Tangent Vectors at the Bifurcation Point

At the bifurcation point $s = \sqrt{5}$, the two tangent vectors are:

    \begin{array}{l}
        \dot{c}^1(\sqrt{5}) = \left(-\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right) \\[0.5em]
        \dot{c}^2(\sqrt{5}) = \left(\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right)
    \end{array}

Observation: The two tangent vectors:

• Have the same $\lambda$-component: $\frac{1}{\sqrt{5}}$ (both branches increase in $\lambda$)
• Have opposite $p_0$-components: $\pm\frac{2}{\sqrt{5}}$ (branches diverge symmetrically)
• Are not collinear, confirming a genuine bifurcation

The angle between the tangent vectors can be computed:

    \cos\theta = \frac{\dot{c}^1 \cdot \dot{c}^2}{\|\dot{c}^1\| \|\dot{c}^2\|} = \frac{-\frac{4}{5} + \frac{1}{5}}{1 \cdot 1} = -\frac{3}{5}

giving $\theta = \arccos(-3/5) \approx 126.87°$.

---

Summary

This parameterized optimal control problem exhibits a transcritical bifurcation at $(p_0, \lambda) = (4, 0)$. For each value of the parameter $\lambda \in [-1, 1] \setminus {0}$, there exist exactly two distinct extremal solutions characterized by different initial costates. At $\lambda = 0$, these two branches merge into a single solution, and the two tangent vectors at this bifurcation point have distinct directions, confirming the bifurcation structure.

[gergaud]

👍

Le 3 oct. 2025 à 20:40, Olivier Cots @.***> a écrit :

ocots left a comment (control-toolbox/OptimalControl.jl#654) https://github.com/control-toolbox/OptimalControl.jl/issues/654#issuecomment-3366824571 Optimal control problem with a bifurcation.

using Pkg Pkg.activate(".")

using OptimalControl using Plots using OrdinaryDiffEq import BifurcationKit: BifurcationProblem, continuation, PALC, ContinuationPar

parameters

t0 = 0 x0 = -2 tf = 1

model

ocp = @def begin λ ∈ R, variable t ∈ [t0, tf], time x ∈ R, state u ∈ R, control x(t0) == x0 x(tf)^2 - λ^2 == 0 ẋ(t) == u(t) ∫( u(t)^2 ) → min end

maximising control

u(x, p, λ) = p/2

Hamiltonian flow

f = Flow(ocp, u);

state projection

π((x, p)) = x

shooting function

S(p0::Number, λ::Number) = π( f(t0, x0, p0, tf, λ) )^2 - λ^2 S(p0::AbstractVector, λ::AbstractVector) = [S(p0[1], λ[1])]

continuation: branch 1

p0_init = 2 λ_init = -1 λ_final = 1 prob = BifurcationProblem(S, Float64[p0_init], Float64[λ_init], 1; record_from_solution = (x,p; k...) -> x[1]) br = continuation(prob, PALC(), ContinuationPar(p_min = Float64(λ_init), p_max = Float64(λ_final))) plt_branch_1 = plot(br; dpi=200)

savefig(plt_branch_1, "plot_branch_1.png")

continuation: branch 2

p0_init = 6 λ_init = -1 λ_final = 1 prob = BifurcationProblem(S, Float64[p0_init], Float64[λ_init], 1; record_from_solution = (x,p; k...) -> x[1]) br = continuation(prob, PALC(), ContinuationPar(p_min = Float64(λ_init), p_max = Float64(λ_final))) plt_branch_2 = plot(br; dpi=200)

savefig(plt_branch_2, "plot_branch_2.png") plot_branch_1.png (view on web) https://github.com/user-attachments/assets/b14ab292-3fe8-4cc8-986c-84d72fc35650 plot_branch_2.png (view on web) https://github.com/user-attachments/assets/1100d507-0569-4e16-9f4a-3c812db5f99c Bifurcation Analysis of Parameterized Optimal Control Problem

1. Problem Formulation

We consider the following parameterized optimal control problem:

$$(P)\left{ \begin{array}{l} \min \int_{0}^{1}u^2,dt \[0.5em] \dot{x}=u \[0.3em] x(0)=-2 \[0.3em] x(1)^2=\lambda^2 \end{array} \right.$$ where:

$x \in \mathbb{R}$ is the state variable $u \in \mathbb{R}$ is the control variable $\lambda \in \mathbb{R}$ is a parameter The initial condition is fixed at $x(0) = -2$ The final condition is parameterized: the state must reach $\pm\lambda$ at time $t=1$ This is a minimum energy problem where we seek to transfer the system from $x(0)=-2$ to a final state satisfying $x(1)^2 = \lambda^2$, minimizing the control effort.

2. Pontryagin Maximum Principle (PMP)

2.1 Maximizing Control (Feedback Form)

The Hamiltonian of the system is:

$$H(t,x,p,u) = -u^2 + pu$$ where $p$ is the costate variable. Applying the maximization condition from PMP:

$$\frac{\partial H}{\partial u} = -2u + p = 0$$ yields the optimal control in feedback form:

$$u(x,p) = \frac{p}{2}$$ 2.2 True Hamiltonian

Substituting the optimal control into the Hamiltonian gives the true Hamiltonian:

$$H_r(t,x,p) = H(t,x,p,u(x,p)) = -\frac{p^2}{4} + p \cdot \frac{p}{2} = \frac{p^2}{4}$$ 2.3 Hamiltonian Dynamics

The state and costate dynamics are governed by:

$$\left{ \begin{array}{l} \dot{x} = \frac{\partial H_r}{\partial p} = \frac{p}{2} \[0.5em] \dot{p} = -\frac{\partial H_r}{\partial x} = 0 \end{array} \right.$$ 3. Extremal Solutions

3.1 General Solution

Since $\dot{p} = 0$, the costate is constant: $p(t) = a \in \mathbb{R}$ for all $t \in [0,1]$.

Integrating the state equation with initial condition $x(0) = -2$:

$$\left{ \begin{array}{ll} x(t) = \frac{at}{2} - 2 & t\in[0,1]\[0.5em] p(t) = a & t\in[0,1] \end{array} \right.$$ 3.2 Determination of the Costate Parameter

The final condition $x(1)^2 = \lambda^2$ determines the value of $a$ as a function of $\lambda$:

$$\begin{array}{rcl} x(1)^2 = \lambda^2 & \Leftrightarrow & \left(\frac{a}{2} - 2\right)^2 = \lambda^2 \[0.5em] & \Leftrightarrow & \frac{a^2}{4} - 2a + 4-\lambda^2 = 0 \[0.5em] & \Leftrightarrow & a^2 - 8a + 16 - 4\lambda^2 = 0 \end{array}$$ Solving this quadratic equation:

$$a = \frac{8 \pm \sqrt{64 - 64 + 16\lambda^2}}{2} = \frac{8 \pm 4\lambda}{2} = 4 \pm 2\lambda$$ This yields two branches:

$$\begin{array}{l} a^1_\lambda = 4 - 2\lambda \[0.3em] a^2_\lambda = 4 + 2\lambda \end{array}$$ 4. Parametrization of the Solution Branches

For $i \in {1, 2}$, we define the initial costate as a function of the parameter:

$$\begin{array}{lllll} p^i_0 & : & [-1,1] & \longrightarrow & \mathbb{R} \ & & \lambda & \longmapsto & p^i_0(\lambda) = a^i_\lambda \end{array}$$ We represent each branch as a curve in the $(p_0, \lambda)$ plane:

$$\begin{array}{lllll} c^i & : & [-1,1] & \longrightarrow & \mathbb{R}^2 \ & & \lambda & \longmapsto & c^i(\lambda) = (p^i_0(\lambda), \lambda) \end{array}$$ Explicitly:

Branch 1: $c^1(\lambda) = (4 - 2\lambda, \lambda)$ Branch 2: $c^2(\lambda) = (4 + 2\lambda, \lambda)$ 5. Arc-Length Parametrization

5.1 Computation of Arc Length

To perform bifurcation analysis, we reparametrize using arc length $s$. Computing the derivatives:

$$\frac{dc^1}{d\lambda} = (-2, 1), \quad \frac{dc^2}{d\lambda} = (2, 1)$$ Both have the same norm: $\left|\frac{dc^i}{d\lambda}\right| = \sqrt{4 + 1} = \sqrt{5}$

The arc length from $\lambda = -1$ to $\lambda$ is:

$$s(\lambda) = \int_{-1}^{\lambda} \left|\frac{dc}{d\mu}\right| d\mu = \int_{-1}^{\lambda} \sqrt{5}, d\mu = \sqrt{5}(\lambda+1)$$ The total arc length for $\lambda \in [-1, 1]$ is $s_{\text{max}} = 2\sqrt{5}$.

5.2 Tangent Vectors in Arc-Length Parametrization

The tangent vectors with respect to arc length are:

$$\begin{array}{lllll} \dot{c}^i & : & [0, 2\sqrt{5}] & \longrightarrow & \mathbb{R}^2 \ & & s & \longmapsto & \dot{c}^i(s) = \frac{1}{\sqrt{5}}\frac{dc^i}{d\lambda} \end{array}$$ Explicitly:

Branch 1: $\dot{c}^1(s) = \left(-\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right)$ Branch 2: $\dot{c}^2(s) = \left(\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right)$ 5.3 Curves in Arc-Length Parametrization

Converting to arc-length parameter where $s = \sqrt{5}(\lambda + 1)$, i.e., $\lambda = \frac{s}{\sqrt{5}} - 1$:

$$\begin{array}{lllll} c^i & : & [0, 2\sqrt{5}] & \longrightarrow & \mathbb{R}^2 \ & & s & \longmapsto & c^i(s) = \left((-1)^i\frac{2s}{\sqrt{5}} + a^i_{-1}, \frac{s}{\sqrt{5}} - 1\right) \end{array}$$ where $a^1_{-1} = p^1_0(-1) = 4 - 2(-1) = 6$ and $a^2_{-1} = p^2_0(-1) = 4 + 2(-1) = 2$.

Branch 1:

$$c^1(s) = \left(-\frac{2s}{\sqrt{5}} + 6, \frac{s}{\sqrt{5}} - 1\right)$$ Branch 2:

$$c^2(s) = \left(\frac{2s}{\sqrt{5}} + 2, \frac{s}{\sqrt{5}} - 1\right)$$ 6. Bifurcation Analysis

6.1 Bifurcation Point

The two branches meet when $c^1(s) = c^2(s)$:

From the $\lambda$-component: both give $\lambda = \frac{s}{\sqrt{5}} - 1$, so we need:

$$-\frac{2s}{\sqrt{5}} + 6 = \frac{2s}{\sqrt{5}} + 2$$ Solving:

$$4 = \frac{4s}{\sqrt{5}} \Rightarrow s = \sqrt{5}$$ At $s = \sqrt{5}$ (corresponding to $\lambda = 0$), the bifurcation point is:

$$c^1(\sqrt{5}) = c^2(\sqrt{5}) = (4, 0)$$ This means:

$p_0 = 4$ (initial costate) $\lambda = 0$ (parameter value) 6.2 Tangent Vectors at the Bifurcation Point

At the bifurcation point $s = \sqrt{5}$, the two tangent vectors are:

$$\begin{array}{l} \dot{c}^1(\sqrt{5}) = \left(-\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right) \[0.5em] \dot{c}^2(\sqrt{5}) = \left(\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}}\right) \end{array}$$ Observation: The two tangent vectors:

Have the same $\lambda$-component: $\frac{1}{\sqrt{5}}$ (both branches increase in $\lambda$) Have opposite $p_0$-components: $\pm\frac{2}{\sqrt{5}}$ (branches diverge symmetrically) Are not collinear, confirming a genuine bifurcation The angle between the tangent vectors can be computed:

$$\cos\theta = \frac{\dot{c}^1 \cdot \dot{c}^2}{|\dot{c}^1| |\dot{c}^2|} = \frac{-\frac{4}{5} + \frac{1}{5}}{1 \cdot 1} = -\frac{3}{5}$$ giving $\theta = \arccos(-3/5) \approx 126.87°$.

Summary

This parameterized optimal control problem exhibits a transcritical bifurcation at $(p_0, \lambda) = (4, 0)$. For each value of the parameter $\lambda \in [-1, 1] \setminus {0}$, there exist exactly two distinct extremal solutions characterized by different initial costates. At $\lambda = 0$, these two branches merge into a single solution, and the two tangent vectors at this bifurcation point have distinct directions, confirming the bifurcation structure.

— Reply to this email directly, view it on GitHub https://github.com/control-toolbox/OptimalControl.jl/issues/654#issuecomment-3366824571, or unsubscribe https://github.com/notifications/unsubscribe-auth/AP2IVRRURG3KADBN544BUP33V27LHAVCNFSM6AAAAACIHH2N6KVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTGNRWHAZDINJXGE. You are receiving this because you were mentioned.

[gergaud]

Olivier,

Le control est calculer pour le coût sans le 0.5 : Il faut donc mettre :

∫( abs(u(t)) - ε*( log(abs(u(t))) + log( max( 1 - abs(u(t)), 1e-8 ) ) ) ) → min

L’état est alors correct

La solution est plus jolie (la commutation est vers 1.3.

Pour avoir une commutation est en t1 il faut prendre xf = 1 - exp(t1-tf)

Le 3 oct. 2025 à 20:30, Olivier Cots @.***> a écrit :

ocots left a comment (control-toolbox/OptimalControl.jl#654) https://github.com/control-toolbox/OptimalControl.jl/issues/654#issuecomment-3366789711 Discrete homotopy versus BifurcationKit homotopy:

using Pkg Pkg.activate(".")

using OptimalControl using Plots using OrdinaryDiffEq using NonlinearSolve

parameters

t0 = 0 tf = 2 x0 = 0 xf = 0.5 umax = 1

model

ocp = @def begin ε ∈ R, variable t ∈ [t0, tf], time x ∈ R, state u ∈ R, control x(t0) == x0 x(tf) == xf ẋ(t) == -x(t)+u(t) -umax ≤ u(t) ≤ umax 0.5∫( abs(u(t)) - ε*( log(abs(u(t))) + log( max( 1 - abs(u(t)), 1e-8 ) ) ) ) → min end

maximising control

function control(x, p, ε)
ψ = -1 + abs(p) u = -(2ε)/(ψ - 2ε - √(ψ^2+4*ε^2)) return p < 0 ? -u : u end

Hamiltonian flow

f = Flow(ocp, control);

state projection, p being the costate

π((x, p)) = x

shooting function

S(p0, ε) = π( f(t0, x0, p0, tf, ε) ) - xf

discrete homotopy on ε

ε_init = 1 ε_final = 1e-4 ε_span = logrange(ε_init, ε_final, 10) p0_span = Float64[] let p0_guess = 0.2 for ε ∈ ε_span prob = NonlinearProblem(S, p0_guess, ε) sol = solve(prob) if !(sol.retcode === SciMLBase.ReturnCode.Success) error("Newton solver failed") end p0_guess = sol.u[1] push!(p0_span, p0_guess) end end

plot of the solutions

plt_sol = plot(; dpi=200) for (p0, ε) ∈ zip(p0_span, ε_span) sol = f((t0, tf), x0, p0, ε; saveat=range(t0, tf, 200)) plot!(plt_sol, sol; label="ε = $(ε)", legend=:none, color=1) end ylims!(plt_sol[3], (0, 1)) plt_sol

plot the homotopy path

plt_path = plot(p0_span, ε_span; label="Discrete homotopy path", xlabel="p₀", ylabel="ε", color=1, linewidth=2, dpi=200) hline!(plt_path, [ε_init]; color=2, label=:none) hline!(plt_path, [ε_final]; color=2, label=:none)

homotopy with BifurcationKit

import BifurcationKit: BifurcationProblem, continuation, PALC, ContinuationPar

F(p0::AbstractVector, λ::AbstractVector) = [S(p0[1], 1-λ[1])]

p0_init = p0_span[1] λ_init = 1 - ε_span[1] λ_final = 1 - ε_span[end]

prob = BifurcationProblem(F, [p0_init], [λ_init], 1; record_from_solution = (x, p; k...) -> x[1]) Δε = abs(ε_final - ε_init) dsmax = Δε/100 dsmin = Δε/1000 ds_init = dsmax options = ContinuationPar(p_min=λ_init, p_max=λ_final, dsmax=dsmax, dsmin=dsmin, ds=ds_init, max_steps=500) br = continuation(prob, PALC(), options)

p0_span_bk = Float64[] ε_span_bk = Float64[] for sol ∈ br.sol push!(p0_span_bk, sol.x[1]) push!(ε_span_bk, 1 - sol.p[1]) end

plot!(plt_path, p0_span_bk, ε_span_bk; label="BifurcationKit", color=3, linewidth=2)

savefig(plt_path, "plot_path_homotopy.png") savefig(plt_sol, "plot_sol_homotopy.png") plot_path_homotopy.png (view on web) https://github.com/user-attachments/assets/283c3ad7-793e-4468-8bab-f5992eef3145 plot_sol_homotopy.png (view on web) https://github.com/user-attachments/assets/e2a7a67a-16d5-48e0-a0c3-d999cec8814a — Reply to this email directly, view it on GitHub https://github.com/control-toolbox/OptimalControl.jl/issues/654#issuecomment-3366789711, or unsubscribe https://github.com/notifications/unsubscribe-auth/AP2IVRVPMY6BAS43CJQXG7L3V26EFAVCNFSM6AAAAACIHH2N6KVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTGNRWG44DSNZRGE. You are receiving this because you were mentioned.

[ocots]

Merci c'est mis à jour.

Olivier,

Le control est calculer pour le coût sans le 0.5 : Il faut donc mettre : ∫( abs(u(t)) - ε*( log(abs(u(t))) + log( max( 1 - abs(u(t)), 1e-8 ) ) ) ) → min L’état est alors correct

La solution est plus jolie (la commutation est vers 1.3.

Pour avoir une commutation est en t1 il faut prendre xf = 1 - exp(t1-tf) …

[jbcaillau]

@ocots c mis à jour mais où 🥲 ?

[jbcaillau]

@gergaud @ocots dans le pb https://github.com/control-toolbox/OptimalControl.jl/issues/654#issuecomment-3366789711

je vois contrôle dans un segment (donc signe quelconque) mais une pénalisation log qui interdit le changement de signe : qu’en est-il ? idem sur terme de pénalisation pour le bord, umax = 1 ou pas ?

[ocots]

@ocots c mis à jour mais où 🥲 ?

Directement dans le commentaire. J'ai changé le critère en principe et les figures.

[ocots]

En fait je ne sais pas à quoi sert la contrainte sur le contrôle j'ai juste repris le problème de Joseph.

Pour l'autre question c'est pas clair. Oui umax vaut 1.

En fait oui il faudrait enlever la contrainte il me semble. Je l'ai enlevé.

@gergaud @ocots dans le pb https://github.com/control-toolbox/OptimalControl.jl/issues/654#issuecomment-3366789711

je vois contrôle dans un segment (donc signe quelconque) mais une pénalisation log qui interdit le changement de signe : qu’en est-il ? idem sur terme de pénalisation pour le bord, umax = 1 ou pas ?

[ocots]

@rveltz Comment on fait pour avoir directement tous les chemins pour le problème avec une bifurcation ? Là j'ai fait 2 calculs.

[joseph-gergaud]

Au départ cette homotopie a été choisie pour que u(x,p) soit lisse contrairement à \lambda|u(t) + (1-\lambda)|u|^2 par exemple et pour le problème de transfert optimal : u est dans R^3 et la contrainte est 0 <= |u| <= 1. On pénalise les 2 bornes de cette contrainte par les barrières log pour être lisse Pour le cas où u est dans R 0 < |u| < 1 n’est plus connexe et il faut que le p0 initial ait le bon signe car on ne peut franchir le 0.

[jbcaillau]

👍🏽 @gergaud OK

• on avait en effet fait ça avec Bilel en dim > 1 (domaine de contrôle pas simplement connexe, mais connexe)
• cf. section 4.2 https://caillau.perso.math.cnrs.fr/research/cele-2012.pdf
• mettre umax = 1 (ou remplacer 1 par umax) dans le code

Au départ cette homotopie a été choisie pour que u(x,p) soit lisse contrairement à \lambda|u(t) + (1-\lambda)|u|^2 par exemple et pour le problème de transfert optimal : u est dans R^3 et la contrainte est 0 <= |u| <= 1. On pénalise les 2 bornes de cette contrainte par les barrières log pour être lisse Pour le cas où u est dans R 0 < |u| < 1 n’est plus connexe et il faut que le p0 initial ait le bon signe car on ne peut franchir le 0.

[ocots]

@jbcaillau En fait j'ai déjà enlevé la contrainte dans le commentaire.