[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 Test BifurcationKit.jl for orbit transfer  #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: #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) <#654 (comment)> 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 <#654 (comment)>, 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) <#654 (comment)> 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 <#654 (comment)>, 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
#654 (comment)

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
#654 (comment)

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.