control-toolbox: solving control problems within Julia

Joseph Gergaud, Jean-Baptiste Caillau, Olivier Cots, Pierre Martinon

► To cite this version:

Joseph Gergaud, Jean-Baptiste Caillau, Olivier Cots, Pierre Martinon. control-toolbox: solving control problems within Julia. FGS 2024 - 21st French-German-Spanish conferences on Optimization, Oviedo University, Jun 2024, Gijón, Asturias, Spain. à paraître. hal-04651475

HAL Id: hal-04651475

https://hal.science/hal-04651475v1

Submitted on 23 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

# control-toolbox: solving optimal control problems within Julia

Jean-Baptiste Caillau, Olivier Cots, Joseph Gergaud, Pierre Martinon, Sophia Sed
## What it's about

* Nonlinear optimal control of ODEs:

$$g(x(t_0), x(t_f)) + \int_{t_0}^{t_f} f^0(x(t), u(t)) dt \rightarrow \min$$

subject to

$$\dot{x}(t) = f(x(t), u(t)), \quad t \in [t_0, t_f]$$

plus boundary, control and state constraints

* Our core interests: numerical & geometrical methods in control, applications

## Where it comes from

* BOCOP: the optimal control solver

* HamPath: indirect and Hamiltonian pathfollowing
* Coupling direct and indirect solvers, examples

## OptimalControl.jl

* Double integrator

* Batch processing

* Goddard problem

## Wrap up

* High level modelling of optimal control problems

* Efficient numerical resolution coupling direct and indirect methods

* Collection of examples

## Future

* New applications (biology, space mechanics, quantum mechanics and more)

* Additional solvers: direct shooting, collocation for BVP, Hamiltonian pathfollowing...

* ... and open to contributions!

* CTProblems.jl
* Open toolbox

* Collection of Julia Packages rooted at OptimalControl.jl

## control-toolbox

<figure><img src="image_8.png" /></figure>

The control-toolbox ecosystem gathers Julia packages for mathematical control and applications. It is an outcome of a research initiative supported by the Centre Inria of Université Côte d'Azur and a sequel to previous developments, notably Bocop and Hampath. See also: ct gallery. The root package is OptimalControl.jl which aims to provide tools to solve optimal control problems by direct and indirect methods.

## Documentation

doc OptimalControl.jl

## Installation

See the installation page.

## Partners

UT
Université
de Toulouse

UNIVERSITÉ
CÔTE D'AZUR

<figure><img src="image_9.png" /></figure>

Inria

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

PROGRAMME
DE RECHERCHE
INTELLIGENCE
ARTIFICIELLE

## Credits (not exhaustive!)

* ADNLPModels.jl

anr
agence nationale
de la recherche
* DifferentiationInterface.jl

* DifferentialEquations.jl

* MLStyle.jl

« Catenoid solution

JuliaCon 2024 »

Powered by Documenter.jl and the Julia Programming Language.

© 2025 control-toolbox
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

# Double integrator: time minimisation

The problem consists in minimising the final time $t_f$ for the double integrator system

$$ \dot{x}_1(t) = x_2(t), \quad \dot{x}_2(t) = u(t), \quad u(t) \in [-1, 1], $$

and the limit conditions

$$ x(0) = (1, 2), \quad x(t_f) = (0, 0). $$

This problem can be interpreted as a simple model for a wagon with constant mass moving along a line without friction.

<figure><img src="image_17.png" /></figure>

First, we need to import the `OptimalControl.jl` package to define the optimal control problem and `NLPModelsIpopt.jl` to solve it.

We also need to import the `Plots.jl` package to plot the solution.

using OptimalControl

<figure><img src="image_18.png" /></figure>
using NLPModelsIpopt

using Plots

Then, we can define the problem

@def ocp begin

tf ∈ R, variable

t ∈ [ 0, tf ], time

x = (q, v) ∈ R², state

u ∈ R, control

tf ≥ 0

-1 ≤ u(t) ≤ 1

q(0) == 1

v(0) == 2

q(tf) == 0

v(tf) == 0

0 ≤ q(t) ≤ 5, (1)

-2 ≤ v(t) ≤ 3, (2)

dot(t) == [ v(t), u(t) ]

tf → min

end
Nota bene

In order to ensure convergence of the direct solver, we have added the state constraints labelled (1) and (2):

$$0 \leq q(t) \leq 5, \quad -2 \leq v(t) \leq 3, \quad t \in [0, t_f].$$

Solve it

sol = solve(ocp)

Method = (:direct, :adnlp, :ipopt)

This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...: 1004

Number of nonzeros in inequality constraint Jacobian.: 0

Number of nonzeros in Lagrangian Hessian............: 202

Total number of variables....................: 304

variables with only lower bounds: 1

variables with lower and upper bounds: 303

variables with only upper bounds: 0

Total number of equality constraints............: 204

Total number of inequality constraints............: 0

inequality constraints with only lower bounds: 0

inequality constraints with lower and upper bounds: 0

inequality constraints with only upper bounds: 0
<table><thead><tr><th>iter</th><th>objective</th><th>inf_pr</th><th>inf_du</th><th>lg(mu)</th><th>||d||</th><th>lg(rg)</th><th>alpha_du</th><th>alpha_pr</th><th>ls</th></tr></thead><tbody><tr><td>0</td><td>1.0000000e-01</td><td>1.90e+00</td><td>0.00e+00</td><td>0.0</td><td>0.00e+00</td><td>-</td><td>0.00e+00</td><td>0.00e+00</td><td>0</td></tr><tr><td>1</td><td>9.9999010e-04</td><td>1.88e+00</td><td>2.32e+01</td><td>-5.7</td><td>1.10e+01</td><td>-</td><td>8.53e-02</td><td>8.99e-03h</td><td>1</td></tr><tr><td>2</td><td>9.9900010e-06</td><td>1.88e+00</td><td>6.44e+03</td><td>-5.7</td><td>1.84e+03</td><td>-</td><td>5.27e-04</td><td>1.07e-04h</td><td>1</td></tr><tr><td>3r</td><td>9.9900010e-06</td><td>1.88e+00</td><td>9.99e+02</td><td>0.3</td><td>0.00e+00</td><td>-</td><td>0.00e+00</td><td>2.68e-07R</td><td>3</td></tr><tr><td>4r</td><td>1.8899160e-03</td><td>1.74e+00</td><td>4.65e+02</td><td>-0.7</td><td>4.77e-01</td><td>-</td><td>4.73e-01</td><td>5.34e-01f</td><td>1</td></tr><tr><td>5r</td><td>1.4658984e-01</td><td>1.73e+00</td><td>2.09e+01</td><td>-1.1</td><td>1.45e-01</td><td>2.0</td><td>9.79e-01</td><td>1.00e+00f</td><td>1</td></tr><tr><td>6r</td><td>2.9193546e-01</td><td>1.70e+00</td><td>4.08e+01</td><td>-0.8</td><td>1.70e-01</td><td>2.4</td><td>1.00e+00</td><td>8.53e-01f</td><td>1</td></tr><tr><td>7r</td><td>3.5836667e-01</td><td>1.70e+00</td><td>4.73e+01</td><td>-0.9</td><td>6.64e-02</td><td>2.9</td><td>1.00e+00</td><td>1.00e+00f</td><td>1</td></tr><tr><td>8r</td><td>5.9834814e-01</td><td>1.70e+00</td><td>7.36e+01</td><td>-0.9</td><td>3.75e-01</td><td>2.4</td><td>9.65e-01</td><td>6.40e-01f</td><td>1</td></tr><tr><td>9r</td><td>7.2398061e-01</td><td>1.70e+00</td><td>9.10e+01</td><td>-1.7</td><td>1.48e-01</td><td>2.8</td><td>1.00e+00</td><td>8.51e-01f</td><td>1</td></tr></tbody></table>

<table><thead><tr><th>iter</th><th>objective</th><th>inf_pr</th><th>inf_du</th><th>lg(mu)</th><th>||d||</th><th>lg(rg)</th><th>alpha_du</th><th>alpha_pr</th><th>ls</th></tr></thead><tbody><tr><td>10r</td><td>8.4184963e-01</td><td>1.71e+00</td><td>4.59e+02</td><td>-1.7</td><td>6.10e-01</td><td>2.3</td><td>1.00e+00</td><td>1.93e-01f</td><td>1</td></tr><tr><td>11r</td><td>9.4835839e-01</td><td>1.72e+00</td><td>7.81e+02</td><td>-1.3</td><td>3.80e-01</td><td>1.8</td><td>1.00e+00</td><td>2.80e-01f</td><td>1</td></tr><tr><td>12r</td><td>1.0283813e+00</td><td>1.71e+00</td><td>8.78e+02</td><td>-1.6</td><td>5.91e-01</td><td>1.4</td><td>1.00e+00</td><td>1.67e-01f</td><td>1</td></tr><tr><td>13r</td><td>1.0666091e+00</td><td>1.66e+00</td><td>2.42e+02</td><td>-2.8</td><td>1.71e-01</td><td>1.8</td><td>1.00e+00</td><td>5.19e-01f</td><td>1</td></tr><tr><td>14</td><td>5.7360280e-01</td><td>1.56e+00</td><td>5.29e+01</td><td>0.3</td><td>1.05e+01</td><td>-</td><td>2.58e-01</td><td>6.17e-02h</td><td>1</td></tr><tr><td>15</td><td>6.1993703e-01</td><td>1.56e+00</td><td>7.19e+01</td><td>-5.7</td><td>4.08e+02</td><td>-</td><td>2.38e-03</td><td>2.47e-03h</td><td>1</td></tr><tr><td>16</td><td>7.2554684e-01</td><td>1.55e+00</td><td>7.04e+01</td><td>-5.7</td><td>1.98e+02</td><td>-</td><td>6.46e-03</td><td>2.79e-03h</td><td>1</td></tr><tr><td>17</td><td>9.7142594e-01</td><td>1.53e+00</td><td>5.07e+01</td><td>-5.7</td><td>6.37e+01</td><td>-</td><td>1.85e-02</td><td>1.07e-02h</td><td>1</td></tr><tr><td>18</td><td>1.6662909e+00</td><td>1.45e+00</td><td>1.50e+01</td><td>0.5</td><td>1.28e+01</td><td>-</td><td>1.10e-01</td><td>5.58e-02f</td><td>1</td></tr><tr><td>19</td><td>7.3295731e+00</td><td>1.70e-01</td><td>4.17e+02</td><td>0.3</td><td>6.42e+00</td><td>-</td><td>6.24e-01</td><td>8.83e-01f</td><td>1</td></tr></tbody></table>

<table><thead><tr><th>iter</th><th>objective</th><th>inf_pr</th><th>inf_du</th><th>lg(mu)</th><th>||d||</th><th>lg(rg)</th><th>alpha_du</th><th>alpha_pr</th><th>ls</th></tr></thead><tbody><tr><td>20</td><td>6.1520840e+00</td><td>1.78e-02</td><td>1.48e+02</td><td>0.1</td><td>1.41e+00</td><td>-</td><td>1.00e+00</td><td>9.21e-01h</td><td>1</td></tr><tr><td>21</td><td>6.2752208e+00</td><td>2.09e-04</td><td>1.72e+01</td><td>-0.5</td><td>2.95e-01</td><td>-</td><td>9.95e-01</td><td>9.90e-01f</td><td>1</td></tr><tr><td>22</td><td>6.2660417e+00</td><td>1.86e-06</td><td>6.24e+00</td><td>-2.2</td><td>1.23e-02</td><td>-</td><td>1.00e+00</td><td>9.91e-01h</td><td>1</td></tr><tr><td>23</td><td>5.6695137e+00</td><td>1.84e-03</td><td>2.90e-02</td><td>-2.8</td><td>5.97e-01</td><td>-</td><td>1.00e+00</td><td>1.00e+00f</td><td>1</td></tr><tr><td>24</td><td>5.5118889e+00</td><td>7.92e-04</td><td>6.45e+02</td><td>-3.3</td><td>5.19e-01</td><td>-</td><td>9.78e-01</td><td>1.00e+00h</td><td>1</td></tr><tr><td>25</td><td>5.4861365e+00</td><td>1.97e-04</td><td>6.37e+03</td><td>-3.8</td><td>5.63e-01</td><td>-</td><td>1.00e+00</td><td>8.15e-01h</td><td>1</td></tr></tbody></table>
<table><tr><td>26</td><td>5.4663359e+00</td><td>4.94e-05</td><td>6.71e+01</td><td>-4.8</td><td>5.00e-01</td><td>-</td><td>1.00e+00</td><td>9.97e-01h</td><td>1</td></tr><tr><td>27</td><td>5.4649557e+00</td><td>7.68e-08</td><td>1.39e-07</td><td>-5.8</td><td>6.24e-03</td><td>-</td><td>1.00e+00</td><td>1.00e+00h</td><td>1</td></tr><tr><td>28</td><td>5.4648096e+00</td><td>1.42e-09</td><td>3.60e-01</td><td>-11.0</td><td>1.33e-03</td><td>-</td><td>9.97e-01</td><td>1.00e+00h</td><td>1</td></tr><tr><td>29</td><td>5.4648096e+00</td><td>8.88e-16</td><td>4.76e-16</td><td>-10.3</td><td>1.55e-06</td><td>-</td><td>1.00e+00</td><td>1.00e+00h</td><td>1</td></tr></table>

Number of Iterations....: 29

<table><thead><tr><th></th><th>(scaled)</th><th>(unscaled)</th></tr></thead><tbody><tr><td>Objective............:</td><td>5.4648095606409361e+00</td><td>5.4648095606409361e+00</td></tr><tr><td>Dual infeasibility.....:</td><td>4.7600990651224600e-16</td><td>4.7600990651224600e-16</td></tr><tr><td>Constraint violation....:</td><td>8.8817841970012523e-16</td><td>8.8817841970012523e-16</td></tr><tr><td>Variable bound violation:</td><td>9.5654382192833509e-09</td><td>9.5654382192833509e-09</td></tr><tr><td>Complementarity........:</td><td>5.0382564513062877e-11</td><td>5.0382564513062877e-11</td></tr><tr><td>Overall NLP error......:</td><td>5.0382564513062877e-11</td><td>5.0382564513062877e-11</td></tr></tbody></table>

<table><tr><td>Number of objective function evaluations</td><td>= 34</td></tr><tr><td>Number of objective gradient evaluations</td><td>= 21</td></tr><tr><td>Number of equality constraint evaluations</td><td>= 34</td></tr><tr><td>Number of inequality constraint evaluations</td><td>= 0</td></tr><tr><td>Number of equality constraint Jacobian evaluations</td><td>= 31</td></tr><tr><td>Number of inequality constraint Jacobian evaluations</td><td>= 0</td></tr><tr><td>Number of Lagrangian Hessian evaluations</td><td>= 29</td></tr><tr><td>Total seconds in IPOPT</td><td>= 0.448</td></tr></table>

EXIT: Optimal Solution Found.

and plot the solution
plot.sol)

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>« Basic example (functional version)</figcaption></figure>

Powered by Documenter.jl and the Julia Programming Language.

<caption>Initial guess options »</caption>
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

# Batch processing

## Introduction

Let us consider a self-replicator model [1] describing the dynamics of a microbial population growing inside a closed bioreactor. The bacterial culture has a constant volume. At the beginning of the experience, there is an initial mass of substrate *S* inside the bioreactor, that is gradually consumed by the bacterial population, and transformed into precursor metabolites *P*. These precursors are intermediate metabolites used to produce proteins—such as ribosomes and enzymes—responsible for specific cellular functions; and metabolites of interest *X* which are excreted from the cell. The proteins forming bacterial cells are divided into three classes *M*, *R* and *Q*, associated with the following cellular functions:

* class *M* proteins of the metabolic machinery, responsible for the uptake of nutrients,

* class *S* from the medium (substrate),

* the production of precursor metabolites *P*,

* and the synthesis of metabolites of interest *X*.

Class *R* proteins of the gene expression machinery (such as ribosomes) actively involved in protein biosynthesis (i.e. in the production of proteins of classes *M*, *Q* and *R*). Class *Q* Growth are rate-independent proteins, such as housekeeping proteins responsible for cell maintenance, and ribosomes not involved in protein synthesis. There is an internal control, *α*, accounting for the behaviour of each
individual cell that has to decide between two pathways (transforming precursors *P* into metabolites *M* or into ribosomes *R*), and one external control leveraging the production of the product *X* (for instance using a light induced control of bioengineered cells).

<figure><img src="image_29.png" /></figure>

After some normalisations, a simplified version of the system (not describing the production of metabolite as we will focus on volume maximisation on this example) can be written in terms of the concentrations of the substrate *s*, precursors *p*, ribosomes *r*, and the of the volume *V* of the bacterial population.[2] Accordingly, the state $\varphi = (s, p, r, V)$ is four-dimensional, and there is only one control, $\alpha$:

$$ \dot{s} = -w_M(s)(1-r)V, $$

$$ \dot{p} = w_M(s)(1-r) - w_R(p)r(p+1), $$

$$ \dot{r} = (\alpha - r)w_R(p)r, $$

$$ \dot{V} = w_R(p)rV, $$
where velocities are taken linear in the concentrations, and where Michaelis-Menten kinetics are assumed:

$$v_R := V_R/V = w_R(p)r, \quad v_M := V_M/V = w_M(s)m,$$

with

$$w_R(p) = \frac{k_R p}{K_r + p}, \quad w_M(s) = \frac{k_m s}{K_m + s}.$$

## Biomass maximisation

We first import the needed packages.

```python
using OptimalControl
using NLPModelsIpopt
using Plots
```

<figure><img src="image_30.png" /></figure>

We are interested in maximising the biomass production [3] (final volume of the bacterial population) over a finite time horizon $[0, t_f]$.

To solve the problem, we first set up the boundary values,

```matlab
t0 = 0
```

<figure><img src="image_31.png" /></figure>

This documentation is not for the latest stable release, but for either the development version or an older release.

Click here to go to the documentation for the latest stable release.
V0 = 0.003

together with parameters and auxiliary functions defining the synthesis rates:

k_r = 1.1

k_m = 1.2

K_r = 1.3

K_m = 1.4

w_r(p) = k_r * p / (K_r + p)

w_m(s) = k_m * s / (K_m + s)

Then we define the optimal control problem setting time, state, control, boundary conditions, state and control constraints, dynamics and Mayer cost:

@def batch begin

t ∈ [ t0, tf ], time

φ = (s, p, r, V) ∈ R⁴, state

α ∈ R, control

s(t0) == s0

p(t0) == p0

r(t0) == r0

V(t0) == V0

s(t) ≥ 0

p(t) ≥ 0
$$0 \le r(t) \le 1$$

$$V(t) \ge 0$$

$$0 \le \alpha(t) \le 1$$

$$\dot{\varphi}(t) = F0(s(t), p(t), r(t), V(t)) + \alpha(t) * F1(s(t), p(t), r(t), V(t))$$

$$V(tf) \rightarrow \max$$

end

The dynamics is indeed affine in the control, $\dot{\varphi} = F_0(\varphi) + \alpha F_1(\varphi)$, with vector fields

$$F0(s, p, r, V) = \begin{bmatrix} -w_m(s) * (1 - r) * V \\ w_m(s) * (1 - r) - w_r(p) * r * (p + 1) \\ -w_r(p) * r^2 \\ w_r(p) * r * V \end{bmatrix}$$

$$F1(s, p, r, V) = \begin{bmatrix} 0 & 0 & w_r(p) * r & 0 \end{bmatrix}$$

<figure><img src="image_32.png" /></figure>

## Direct solve

Since the following result holds,[4]

**Proposition.** The Lie bracket $F_{101}$ belongs to the span of $F_1$ and $F_{01}$, so singular controls are at least of local order two.

one expects singular arcs connected with bang arcs through Fuller phenomenon (accumulation of switching times).
We first solve the problem using a uniform discretisation:

```pascal
sol0 = solve(batch; grid_size=1000, print_level=0)
println("Objective ", sol0/objective, " after ", sol0.iterations, " iterations")
```

```pascal
Method = (:direct, :adnlp, :ipopt)
Objective 0.07683200113875356 after 111 iterations
```

Although convergence is obtained, it is actually more efficient to first solve on a raw grid, then use a *warm start* to solve again on a finer (still uniform) grid:

```pascal
sol1 = solve(batch; grid_size=20, print_level=0)
println("Objective ", sol1/objective, " after ", sol1.iterations, " iterations")
```

```pascal
Method = (:direct, :adnlp, :ipopt)
Objective 0.0769994219213667 after 31 iterations
```

```pascal
sol2 = solve(batch; grid_size=1000, print_level=0, init=sol1)
println("Objective ", sol2/objective, " after ", sol2.iterations, " iterations")
```

```pascal
Method = (:direct, :adnlp, :ipopt)
Objective 0.07683200113876136 after 22 iterations
```

# Plotting
We eventually plot the solutions (raw grid + finer grid) and observe that the control exhibits the expected structure with a Fuller-in arc followed by a singular one, then a Fuller-out arc:

```matlab
plot.sol1; solution_label="(N=20)", size=(800, 1000))
plot.sol2; solution_label="(N=1000)"
```

# N is the grid size

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>
<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

## References
* 1 Giordano, N.; Mairet, F.; Gouzé, J.-L.; Geiselmann, J.; De Jong, H. Dynamical allocation of cellular resources as an optimal control problem: novel insights into microbial growth strategies. *PLoS comp. biol.* **12** (2016), e1004802.

* 2 Yabo, A. G.; Caillau, J.-B.; Gouzé, J.-L.; de Jong, H.; Mairet, F. Dynamical analysis and optimization of a generalized resource allocation model of microbial growth. *SIAM J. Appl. Dyn. Syst.* **21** (2022), no. 1, 137-165.

* 3 Yabo, A. G.; Caillau, J.-B.; Gouzé, J.-L. Optimal bacterial resource allocation strategies in batch processing. *SIAM J. Appl. Dyn. Syst.*, to appear.

* 4 Astruc, L.; Edery, N. Optimal allocation of bacterial resources in a bioreactor. Project report, Polytech Nice Sophia, Université Côte d'Azur (2023).

« Goddard problem

Solar sail »

Powered by Documenter.jl and the Julia Programming Language.

© 2025 control-toolbox