sprint #10 / use cases: simple + advanced tests

[jbcaillau]

Extension of BifurcationKit.jl for General Boundary Value Problems from Optimal Control

Context and Motivation

BifurcationKit.jl (BK.jl) currently specializes in periodic boundary value problems. This issue proposes extending the framework to support general boundary value problems (BVPs) arising from optimal control via the indirect method (Pontryagin Maximum Principle + shooting).

The goal is to leverage BK.jl's powerful numerical methods not only for path following and bifurcation detection but also for solving the shooting equations themselves by providing the Hamiltonian system's right-hand side directly rather than a pre-assembled shooting function.

Proposed Use Cases

Use Case 1: Single Arc - Simple Case ✅ READY FOR TESTING

Reference: control-toolbox/OptimalControl.jl#654 (comment)

Problem: Double integrator time minimization

Minimize: tf
Subject to:
    ẋ₁ = x₂
    ẋ₂ = u,  u ∈ [-1, 1]
    x(0) = (-1, 0)
    x(tf) = (0, 0)

Hamiltonian formulation: The Hamiltonian of the system is:

H(x, p, u) = p₁x₂ + p₂u

BVP Structure:

• State dimension: 2
• Costate dimension: 2
• Total system dimension: 4 (x₁, x₂, p₁, p₂)
• Unknowns: Initial costate p(0) and final time tf (3 unknowns)
• Boundary conditions:
  • x(0) = (-1, 0) — fixed initial state (2 equations)
  • x(tf) = (0, 0) — fixed final state (2 equations)
  • H(tf) = 0 — free final time transversality condition (1 equation)
• Control via Maximum Principle: u(t) = -sign(p₂(t))

Current Status:

• ✅ Hamiltonian RHS is fully ready
• ✅ Can be immediately tested with BK.jl
• 🎯 Objective: Use BK.jl to solve the shooting equations by providing only the Hamiltonian vector field

Expected BK.jl capabilities to test:

• Solving the BVP with Newton-Krylov methods
• Handling split boundary conditions (initial + final)
• Managing the free time parameter (tf as variable)

---

Use Case 2: Single Arc - Advanced Case 🔧 REQUIRES HAMILTONIAN GENERATION

Reference: https://control-toolbox.org/Kepler.jl/stable/#Shooting-(2/2),-Tmax-0.7-Newtons

Problem: Minimum time orbit transfer with low thrust

Minimize: tf
Subject to:
    ẍ = -μ(x/|x|³) + (u/m)
    ṁ = -β|u|
    |u| ≤ Tmax
    Initial and final Keplerian orbits (6D state)

Using Keplerian elements representation:

State: x = (P, ex, ey, hx, hy, L) ∈ R⁶
Control: u ∈ R³ with |u| ≤ Tmax
Dynamics: ẋ = F₀(x) + (T/mass)[u₁F₁(x) + u₂F₂(x) + u₃F₃(x)]

BVP Structure:

• State dimension: 6 (Keplerian orbital elements)
• Costate dimension: 6
• Total system dimension: 12
• Parameter: Tmax (thrust magnitude limit)
• Unknowns: Initial costate p(0), final time tf, final longitude Lf
• Boundary conditions:
  • Initial: all 6 components fixed
  • Final: 5 components fixed (P, ex, ey, hx, hy), free longitude L(tf)

Current Status:

• ⚠️ Hamiltonian symplectic gradient needs to be generated
• 📊 Expected bifurcations: Rich behavior when following solution path w.r.t. parameter Tmax
• 🔍 Reference: Swallowtail singularity in value function (see ECC 2021 paper, Fig. 5)

Expected BK.jl capabilities:

• Continuation in parameter Tmax
• Detection of fold/turning points in solution branches
• Handling singular control structures (bang-bang arcs)
• Bifurcation analysis of optimal trajectories

---

Use Case 3: Two Arcs - Simple Case 🏗️ MULTI-ARC STRUCTURE

Reference: https://control-toolbox.org/OptimalControl.jl/stable/example-double-integrator-time.html#Indirect-method

Problem: Double integrator with control switching

Same dynamics as Use Case 1, but with:
- Arc 1: u = -1 on [0, τ]
- Arc 2: u = +1 on [τ, tf]
- Unknown switching time τ

BVP Structure:

• Arc 1:
  • Dynamics: ẋ = [x₂, -1], t ∈ [0, τ]
  • Hamiltonian: H₁ = p₁x₂ - p₂
• Arc 2:
  • Dynamics: ẋ = [x₂, +1], t ∈ [τ, tf]
  • Hamiltonian: H₂ = p₁x₂ + p₂
• Junction conditions at t = τ:
  • State continuity: x₁(τ⁻) = x₂(τ⁺)
  • Costate continuity: p₁(τ⁻) = p₂(τ⁺)
  • Switching function: p₂(τ) = 0

Unknowns:

• Initial costates for both arcs: p₁(0), p₂(0)
• Junction time τ
• Final time tf

Current Status:

• 🏗️ More advanced: Concatenation of two different vector fields
• 🔗 Interior point conditions: Continuity + switching condition
• 🎯 Challenge: Multi-phase structure with free transition time

Expected BK.jl capabilities:

• Handling piecewise-defined systems
• Interior point constraints (junction conditions)
• Multiple free parameters (τ, tf)
• Continuation through structural changes (switching structure)

---

Technical Requirements for BK.jl Extension

1. Input Interface

Instead of providing a pre-assembled shooting function, provide:

• Hamiltonian vector field: H_rhs(x, p, t, params)
• Boundary conditions: Initial and final (split conditions)
• Free parameters: (tf, and potentially others like junction times)

2. BVP Features to Support

Feature	Use Case 1	Use Case 2	Use Case 3
Split boundary conditions	✓	✓	✓
Free final time	✓	✓	✓
Parameter continuation	—	✓ (Tmax)	—
Multi-arc structure	—	—	✓
Interior point conditions	—	—	✓ (junction)
Symplectic integration	Optional	Preferred	Optional

3. Automatic Differentiation

• Support ForwardDiff for Jacobian computation
• Handle adaptive ODE solvers correctly (see ForwardDiff + ODE discussion)

4. Integration with OptimalControl.jl

• Use the Flow infrastructure from CTFlows.jl
• Leverage variable support in Hamiltonian flows
• Consistent interface with existing indirect method tools

---

References

• OptimalControl.jl Documentation
• Kepler.jl Examples
• BK Sprint Issue #654
• Orbit Transfer Bifurcations #224
• Caillau, Daoud, Gergaud (2021). "On Some Riemannian Aspects of Two and Three-Body Controlled Problems", ECC 2021

---

Status Summary:

• ✅ Use Case 1: Ready to test immediately
• 🔧 Use Case 2: Needs Hamiltonian generation, bifurcations expected
• 🏗️ Use Case 3: Most advanced, multi-arc structure

[jbcaillau]

@rveltz au boulot bordel !!! 🫢

[jbcaillau]

TODO: when bvp is ready

• test example no. 1 without using OC.jl flow
• idem for example no. 2 (continuation on Tmax)