Applications

[jbcaillau]

Bactery growth

@agustinyabo Try this. The length three bracket $H_{101}$ not being zero at this particular point, it is not identically zero (while it is zero on ${H_1=H_{01}=0}$ as your previous computation shows). So the singular is of local order two.

# batch.jl
using OptimalControl

const kᵣ = 1.1
const kₘ = 1.2
const Kᵣ = 1.3
const Kₘ = 1.4

wᵣ(p) = kᵣ * p / (Kᵣ + p)
wₘ(s) = kₘ * s / (Kₘ + s)

F0(x) = begin
    s, p, r, V = x
    return [ -wₘ(s) * (1 - r) * V
              wₘ(s) * (1 - r) - wᵣ(p) * r * (p + 1)
             -wᵣ(p) * r^2
              wᵣ(p) * r * V ]
end

F1(x) = begin
    s, p, r, V = x
    return [ 0, 0, wᵣ(p) * r, 0 ]
end

H0 = Lift(F0)
H1 = Lift(F1)

H01 = @Lie { H0, H1 }
H101 = @Lie { H1, H01 }

x = [ 1, 1, 1, 1 ]
λ = [ 1, 1, 1, 1 ]

H101(x, λ)

[jbcaillau]

Application ravageurs

@fkemayou @PierreMartinon

• installer la dernière release stable (1.10.2)
  https://julialang.org/downloads/

using OptimalControl

tf = 15 # and other params

@def ocp begin
       t ∈ [ 0, tf ], time
       x ∈ R³,  state
       u ∈ R, control
       S = x₁
       I = x₂
       N = x₃
       ẋ(t) == [ λ-β*(1-u(t))*S(t)*N(t)-μ*S(t), 0, 0 ] # please update!
       S(0) == 10
       I(0) == 2
       N(0) == 15
       0 ≤ u(t) ≤ 1
       ∫( u(t) - b*S(t) ) → min
end

sol = solve(ocp)

plot(sol)

• @fkemayou ajoute la description de ton problème (par ex. lien vers un papier le décrivant)
• documentation en ligne (par ex. sur fonction solve)
• possibilité de résoudre par une méthode indirecte (voir cet exemple)
• possibilité en Julia (pas avec OptimalControl.jl) de traiter aussi des pb de contrôle sur des EDP (avec JuMP.jl)
• @PierreMartinon ajouter un lien vers un exemple avec warm start (par ex. pour solve avec nouveau tf)

[fkemayou]

@jeanbaptiste et @PierreMartinon j'ai pourtant la version (1.10.2) de julia. Je parviens bien a installer ceratins packages (using plots, using DifferentialEquations) mais au moment d'installer "import Pkg; Pkg.add("OptimalControl")" j'obtiens le message d'erreur suivant depuis hier mardi

The following package names could not be resolved:

• OptimalControl (not found in project, manifest or registry)
  Suggestions: TwoStageOptimalControl RobustAndOptimalControl

Stacktrace:
[1] pkgerror(msg::String)
@ Pkg.Types C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\Types.jl:70
[2] ensure_resolved(ctx::Pkg.Types.Context, manifest::Pkg.Types.Manifest, pkgs::Vector{Pkg.Types.PackageSpec}; registry::Bool)
@ Pkg.Types C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\Types.jl:1030
[3] ensure_resolved
@ C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\Types.jl:981 [inlined]
[4] add(ctx::Pkg.Types.Context, pkgs::Vector{Pkg.Types.PackageSpec}; preserve::Pkg.Types.PreserveLevel, platform::Base.BinaryPlatforms.Platform, kwargs::@kwargs{io::IJulia.IJuliaStdio{Base.PipeEndpoint}})
@ Pkg.API C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\API.jl:267
[5] add(pkgs::Vector{Pkg.Types.PackageSpec}; io::IJulia.IJuliaStdio{Base.PipeEndpoint}, kwargs::@kwargs{})
@ Pkg.API C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\API.jl:159
[6] add(pkgs::Vector{Pkg.Types.PackageSpec})
@ Pkg.API C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\API.jl:148
[7] add
@ C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\API.jl:147 [inlined]
[8] add(pkg::String)
@ Pkg.API C:\Users\PC.julia\juliaup\julia-1.10.2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Pkg\src\API.jl:146
[9] top-level scope
@ In[5]:1
/

[PierreMartinon]

@fkemayou Bonjour Frank,
Oui, c'est normal, notre package OptimalControl n'est pas [encore] dans le registre officiel de Julia.
Une fois dans le pkg> (avec la touche ], et backspace pour revenir a julia>), avant d'ajouter OptimalControl, il faut lui donner notre registre

pkg> registry add https://github.com/control-toolbox/ct-registry.git
pkg> add OptimalControl

Normalement ca marche aussi sans entrer dans le pkg>

julia> import Pkg; Pkg.Registry.add(RegistrySpec(url = "https://github.com/control-toolbox/ct-registry.git")); Pkg.add("OptimalControl")

[jeanbaptiste]

Bonjour @PierreMartinon je pense que vous souhaitez ajouter @jbcaillau plutôt que moi sur ce thread. Bravo pour votre travail sur CT. Bien à vous, @jeanbaptiste (un ancien de l'inria)

[PierreMartinon]

Bonjour Jean-Baptiste, en effet, desole !

[jeanbaptiste]

Pas de souci! Je découvre BTW qu'on peut juste mentionner qui on veut sur un thread et on peut poster dessus sans aucune validation... curieux de voir ce que ça va donner niveau spam 😉

[fkemayou]

Bonjour Jean Baptiste Caillau et Pierre Martinon

Ci joint mon code julia qui marche parfaitement jusqu’à tf = 70 mais des lors que je passe a tf = 80 je rencontre des problèmes de non convergence pourtant je dois arriver jusqu’à tf = 300. Ma question est donc la suivante: comme sur le v2 comment redémarrer d'une solution initiale en utilisant par exemple #init=OptimalControlInit() dans la fonction solve.

using OptimalControl
using Plots
tf = 80 # and other params
λ = 300
β =1e-07
μ = 0.05
d = 0.0002
a = 60
r = 0.15
rho = 400
K = 1000
nu = 0.045
b1 = 2500
b2 = 2
b3 = 5 
@def ocp begin
       t ∈ [ 0, tf ], time
       x ∈ R³,  state 
       u ∈ R, control
       S = x₁
       I = x₂
       N = x₃
       ẋ(t) == [ λ-β*(1-u(t))*S(t)*N(t)-μ*S(t), 
                 β*(1-u(t))*S(t)*N(t)-μ*I(t)-d*N(t)*I(t)/(a+I(t)), 
                 (r+rho*d*I(t)/(a+I(t)))*N(t)*(1-N(t)/(K*I(t)))-nu*N(t) ] 
       S(0) == 100
       I(0) == 100
       N(0) == 1000
       0 ≤ u(t) ≤ 1
       b3*I(tf) + ∫(b1*u(t) - b2*S(t) ) → min 
end

sol = solve(ocp, grid_size = 100, max_iter = 2000, tol = 10e-20, print_level=5, mu_strategy="adaptive")
plot(sol)

[PierreMartinon]

Bonjour Frank,
Regarde le tutoriel sur la documentation 'dev': https://control-toolbox.org/CTDirect.jl/dev/continuation.html
On va essayer de faire une release bientot, mais quelque chose comme

init=OptimalControlInit(sol)

devrait marcher avec la version actuelle de CTDirect.
Pierre

edit: release 0.4.6 faite, tu peux mettre a jour le package avec la commande 'up' dans le package REPL (accessible avec ])

[jbcaillau]

@fkemayou @PierreMartinon bonsoir tous les deux, merci frank pour les infos / tests et pierre pour la nouvelle release.

• @fkemayou mets stp un lien vers un papier décrivant ton problème et ajoute quelques plots de tes simus qu'on voie à quoi ça ressemble. pierre t'a indiqué comment faire un warm start (résoudre en reprenant la solution précédente), tiens nous au courant 🤞🏾

[fkemayou]

Bonjour @Jbcaillau et @PierreMartinon. Ci-joint un le lien vers un papier décrivant mon problème

papierfrank

Aussi, je n'arrive pas toujours a faire le warm start. J'obtiens toujours cette erreur

UndefVarError: OptimalControlInit not defined

[fkemayou]

Dans ce papier, j'ai ajoute un exemple de simul que j'obtiens. Bref je cherche toujours a faire le warm start afin de voir comment ca marche. Si vous pouver m'aider avec un peu plus d'indications ou un exemple concret, cela me conviendrait.

Merci

[jbcaillau]

Bonjour @fkemayou

• pour l'init, je laisse @PierreMartinon répondre (mise à jour, de fait, cf. cette issue)
• pour le papier, où est le lien vers la référence dont tu parles ?

[PierreMartinon]

Bonjour @fkemayou ,
Le premier exemple de la page 'Continuation' montre justement une continuation sur le temps final
https://control-toolbox.org/CTDirect.jl/stable/continuation.html

[fkemayou]

C'est ca le lien
papierfrank

[jbcaillau]

merci @fkemayou ; par contre, l'accès à ce doc sur drive google demande une permission d'accès. dépose plutôt le pdf dans cette discussion, stp.

[fkemayou]

Merci @jbcaillau
Le pdf en fichier joint

OCP_Frank.pdf

[jbcaillau]

@PierreMartinon peux-tu stp renvoyer à @fkemayou un exemple simple de warm start ?

[fkemayou]

@fkemayou Pour ton export de solution, c'est pour faire quoi exactement ? La sauvegarde sous julia fonctionne mais c'est un format qu'on ne peut pas lire en dehors. Le format texte (JSON) est en cours.

Aussi vu que je ne trouve pas d'arc singulier, comment montrer alors que c'est la meilleure solution optimale?

Sans etre une preuve, pour un probleme de petite dimension sur l'etat on peut essyer une methode globale comme bocopHJB et voir si la solution ressemble. Ici, n=3, avec une discretisation grossiere ca pourrait marcher. Le code doit etre sur http://bocop.org, il doit etre plus facile a installer que bocop2 car il n'y a ni solver NLP ni differentiation automatique.

@PierreMartinon , Mon export était exactement pour les points suivants:

1. Comment augmenter la taille des figures afin que les valeurs de l'axe des ordonnées ne soient bien visibles vu qu'elles sont très serrées.
2. Est - il aussi possible d'ajouter l'épaisseur des courbes. ( j'ai essaye un truc comme ca " inewidth " mais ca ne marche pas)
3. Généralement en contrôle optimal, nous aimons bien souvent avoir sur la même figure les courbes lorsque le contrôle est nul et lorsqu'on applique le contrôle optimal. Est il alors possible de le faire directement dans ce cas?
   Ce sont les raisons de mon export. Avec la v2, en faisais le export, je pouvais facilement faire tout ceci afin d'avoir des courbes et peu plus propre.

[ocots]

Pour les paramètres des plots, il y a un tutoriel ici : https://control-toolbox.org/docs/optimalcontrol/stable/tutorial-plot.html

Ça m'étonne que l'on ne puisse pas changer l'épaisseur des courbes directement. Je vais checker ça.

[fkemayou]

@fkemayou as-tu mis ton code / notebook visible qqpart sur un repo git ?

Bonjour @jbcaillau pour le code, il n'est visible que sur un report git !!!!

[jbcaillau]

@fkemayou OK, sur quel repo (public) est-il ?

[fkemayou]

@fkemayou OK, sur quel repo (public) est-il ?

Bonjour @jbcaillau, désole du retard j'avais des problèmes avec mon PC. Le code est en effet sur ce report public fkemayou/example

[jbcaillau]

@ocots check this example

[jbcaillau]

@rand-asswad n'hésite pas à utiliser cette discussion pour échanger avec @PierreMartinon et les autres personnes impliquées dans control-toolbox. On est bien sûr intéressés par ton retour d'expérience, si tu as un exemple pour alimenter les exemples / tutos déjà présents ici, parfait 🤞🏾

[jbcaillau]

control-toolbox/CTDirect.jl#111 (comment)

@rand-asswad n'hésite pas à utiliser cette discussion pour échanger avec @PierreMartinon et les autres personnes impliquées dans control-toolbox. On est bien sûr intéressés par ton retour d'expérience, si tu as un exemple pour alimenter les exemples / tutos déjà présents ici, parfait 🤞🏾

[jbcaillau]

@fkemayou tu peux merger cette PR fkemayou/example#1 (cliquer sur le bouton merge)

[ocots]

@fkemayou @agustinyabo Bonjour ! Si vous souhaitez ajouter une application à la liste des applications que l'on trouve sur la documentation de OptimalControl.jl, vous êtes les bienvenus.

Pour cela, il vous suffit de suivre le tutoriel ici.

[agustinyabo]

C'est beau, je regarde ça dès que possible !

[jbcaillau]

@agustinyabo please restart from scratch and create a new repo from @ocots tutorial (this repo will replace the current one) 🤞🏾

[fkemayou]

Très bien @ocots, je le regarde ça

[jbcaillau]

@ocots quantum stuff
https://arxiv.org/abs/2205.15044

[amontoison]

@jbcaillau
An opportunity to sell your packages:
https://discourse.julialang.org/t/optimizing-a-function-in-optimal-control-problem/125483/14

[jbcaillau]

thanks @amontoison @rveltz for the heads-up 👍🏽
https://discourse.julialang.org/t/optimizing-a-function-in-optimal-control-problem/125483/15

[jbcaillau]

@rveltz @amontoison https://discourse.julialang.org/t/optimizing-a-function-in-optimal-control-problem/125483/16?u=jbcaillau

[rveltz]

My plaisir non dissimulated :D

[jbcaillau]

@ocots rings a bell 🙂?
https://discourse.julialang.org/t/how-to-make-a-locally-defined-function-type-stable/125639

[jbcaillau]

@amontoison again, yes; exactly that casadi/casadi#386 (comment)

as mentioned by @PierreMartinon, any variant of the Goddard problem with an integral cost and free final time will do. minimising the $L^1$ norm of $u$ is indeed a good example as it is what people want to do in practice; you can also think of a (squared) $L^2$ norm t regularise the cost, bref.

note that:

• keeping the integral cost instead of adding an additional state will reduce by $N$ (= grid size) the number of unknowns of the problem, up to adding a full row
• similarly, one could conversely remove the dense column by changing the final time into an additional state, up to the price of adding $N$ variables

[jbcaillau]

@amontoison some details below. @PierreMartinon and @ocots feel free to complete / correct.

Consider an optimal control with state $x(t) \in \mathbf{R}^n$, control $u(t) \in \mathbf{R}^m$, and an extra scalar parameter $y \in \mathbf{R}$ (the same reasonning readily extends to a vector); this parameter, that is also to be optimised, could for instance be the final time, if it is free, of the problem. Assume that the dynamics is modelled by the ODE

$$ \dot{x}(t) = f(x(t), u(t), y), $$

that also depends on the parameter, with $f : \mathbf{R}^n \times \mathbf{R}^m \times \mathbf{R} \to \mathbf{R}^n$ is a smooth function. Using a one-step numerical scheme to discretise this ODE on a time grid $t_0, t_1, \dots, t_N$ of size $N + 1$ results in a set of equality constraints; for instance, with a Crank-Nicolson (trapezoidal) scheme, denoting $h_i = t_{i+1} - t_i$, one has ($x_i \simeq x(t_i)$, $u_i \simeq u(t_i)$)

$$ x_{i+1} - x_i - (h_i/2) ( f(x_i, u_i, y) + f(x_{i+1}, u_{i+1}, y) ) = 0,\quad i = 0, \dots, N-1. $$

Let $F(X, U, y)$ be the associated constraint function on $X = (x_0, \dots, x_N) \in \mathbf{R}^{(N+1)n}$, $U = (u_0, \dots, u_N) \in \mathbf{R}^{(N+1)m}$ and $y$. Clearly, $\partial F/\partial y$ is a dense column in the $Nn \times (N + 1)(n + m + 1)$ Jacobian of $F$. This dense column can be eliminated by replacing the parameter $y$ with a constant additional state, $y(t) \in \mathbf{R}$. The augmented ODE on the extended state $(x(t),y(t)) \in \mathbf{R}^{n+1}$ is

$$ \dot{x}(t) = f(x(t), u(t), y(t)),\quad \dot{y}(t) = 0, $$

and one can use the same scheme as before up to the price of adding the $N+1$ new unknowns $Y = (y_0, \dots, y_N)$ to the discretised problem.

Similarly, assume that a (scalar) integral constraint such ($t_f$ begin the final time)

$$ \int_0^{t_f} g(x(t), u(t)),\mathrm{d}t = 0 $$

is to be discretised. This can be accomodated using any quadrature, for instance a trapezoidal one to be consistant with the discussion on the ODE above:

$$ \sum_{i=0}^{N-1} (h_i/2) ( g(x_i, u_i) + g(x_{i+1}, u_{i+1}) ) = 0. $$

The associated constraint $G(X, U)$ has a Jacobian that is a dense row vector. As before, it is possible to eliminate this dense row by extending the state with an extra coordinate, $\xi(t)$, and augmenting the dynamics according to

$$ \dot{\xi}(t) = g(x(t), u(t)), $$

and the two boundary constraints $\xi(0) = \xi(t_f) = 0$. Again, the price to pay is adding $N+1$ variables, $\Xi = (\xi_0, \dots, \xi_N)$ to the discretised problem, associated with the linear constraints (with sparse Jacobians!) $\xi_0 = \xi_N = 0$.

[amontoison]

@ocots @jbcaillau @PierreMartinon
What do you think of the explanations?

[jbcaillau]

@amontoison it's OK. minor remarks:

• for the first point, I think it's clearer to illustrate the first point with another parameter than $t_f$ (rather banal situation); when it is indeed $t_f$, it is because once has rescaled on $[0,1]$ setting $s = t/t_f$ so that the dynamics depends on $t_f$ in a very particular way

$$ dx/ds = t_f f(x, u). $$

• I suggest to keep $x_i$ (and associated vector $X := (x_0, \dots, x_N)$, etc.) instead of $x(t_i)$; it is easier to understand when the actual finite dim variables are explicited

[amontoison]

@jbcaillau I have this version now.

[jbcaillau]

@Weng-W0 code below for your example (this one converges)

using OptimalControl
using NLPModelsIpopt
using Plots

const T = 10.0
const x0 = [1.0, 0.0] 
const α = [0.5, 0.35]

ocp = @def begin
    t ∈ [0, T], time
    x ∈ R², state
    u ∈ R, control
    x(0) == x0
    ẋ(t) == [x₂(t), u(t)]
    ∫( 0.5u(t)^2 + F(x(t), α) ) → min
end

function F(x, α)
    r = 0.
    for i ∈ 1:length(α)
        r = r + abs(x[i])^(1 + α[i]) / (1 + α[i]) # non-smooth!
    end
    return r
end

sol = solve(ocp)

plot(sol)

plus a link towards a shooting example:
https://control-toolbox.org/OptimalControl.jl/stable/tutorial-iss.html

[ocots]

@jbcaillau Do we close? Where are we with this discussion?

[jbcaillau]

i'd say we leave open for general discussion on further discussions (this is not an issue) on applications