Skip to content

Commit

Permalink
foo
Browse files Browse the repository at this point in the history
  • Loading branch information
@ocots
ocots committed Oct 5, 2023
1 parent 53b2ac6 commit 2cf5dfa
Show file tree
Hide file tree
Showing 2 changed files with 125 additions and 56 deletions.
17 changes: 15 additions & 2 deletions docs/src/juliaopt2023.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,47 +3,60 @@
```

# Solving optimal control problems with Julia

### [Jean-Baptiste Caillau](http://caillau.perso.math.cnrs.fr), [Olivier Cots](https://ocots.github.io), [Joseph Gergaud](https://scholar.google.com/citations?user=pkH4An4AAAAJ&hl=fr), [Pierre Martinon](https://www.linkedin.com/in/pierre-martinon-b4603a17), [Sophia Sed](https://iww.inria.fr/sed-sophia)

```@raw html
<img width="800" alt="affiliations" src="./assets/affil.jpg">
```

# What it's about

- Nonlinear optimal control of ODEs:

```math
g(x(t_0),x(t_f)) + \int_{t_0}^{t_f} f^0(x(t), u(t)) \to \min
```

subject to

```math
\dot{x}(t) = f(x(t), u(t)),\quad t \in [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](https://www.bocop.org)
- [HamPath: indirect and Hamiltonian pathfollowing](http://www.hampath.org)
- [Coupling direct and indirect solvers, examples](https://ct.gitlabpages.inria.fr/gallery//notebooks.html)

# OptimalControl.jl

- [Basic example: double integrator (1/3)](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-basic-example-f.html)
- [Basic example: double integrator (2/3)](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-basic-example.html)
- [Basic example: double integrator (3/3)](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-double-integrator.html)
- [Indirect simple shooting](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-iss.html)
- [Advanced example: Goddard problem](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-goddard.html)

# Wrap up

- [X] High level modelling of optimal control problems
- [X] Efficient numerical resolution coupling direct and indirect methods
- [X] Collection of examples
- [X] Collection of examples

# Future

- [ct_repl](./assets/repl.mp4)
- Additional solvers: direct shooting, collocation for BVP, Hamiltonian pathfollowing...
- ... and open to contributions!
- [CTProblems.jl](https://control-toolbox.org/docs/ctproblems/stable/problems-list.html)

# control-toolbox.org

- Open toolbox
- Collection of Julia Packages rooted at [OptimalControl.jl](https://control-toolbox.org/docs/optimalcontrol)

Expand All @@ -52,6 +65,7 @@ plus boundary, control and state constraints
```

# Credits (not exhaustive!)

- [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl)
- [JuMP](https://jump.dev/JuMP.jl),
[InfiniteOpt.jl](https://docs.juliahub.com/InfiniteOpt/p3GvY/0.4.1),
Expand All @@ -62,4 +76,3 @@ plus boundary, control and state constraints
[Zygote.jl](https://fluxml.ai/Zygote.jl))
- [MLStyle.jl](https://thautwarm.github.io/MLStyle.jl)
- [REPLMaker.jl](https://docs.juliahub.com/ReplMaker)

164 changes: 110 additions & 54 deletions docs/src/tutorial-iss.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,20 @@ using OptimalControl
using MINPACK # NLE solver
```

Let us consider the following optimal control problem.
Let us consider the following optimal control problem:

```math
\left\{
\begin{array}{l}
\min \displaystyle \frac{1}{2} \int_{t_0}^{t_f} {u(t)}^2 \, \mathrm{d} t\\[1.0em]
\dot{x}(t) = \displaystyle -x(t)+\alpha x^2(t)+u(t), \quad u(t) \in \R,
\quad t \in [t_0, t_f] \text{ a.e.}, \\[0.5em]
x(t_0) = x_0, \quad x(t_f) = x_f,
\end{array}
\right.%
```

with $t_0 = 0$, $t_f = 1$, $x_0 = -1$, $x_f = 0$, $\alpha=1.5$ and $\forall\, t \in [t_0, t_f]$, $x(t) \in \R$.

```@example main
t0 = 0
Expand All @@ -33,13 +46,49 @@ end;
nothing # hide
```

From the Pontryagin Maximum Principle, we have that the maximising control is given by
The **pseudo-Hamiltonian** of this problem is

```math
H(x,p,u) = p \, (-x+u) + p^0 u^2 /2,
```

where $p^0 = -1$ since we are in the normal case. From the Pontryagin Maximum Principle, the maximising control is given by

```math
u(x, p) = p
```

and we can define the flow of the associated Hamiltonian vector field
since $\partial^2_{uu} H = p^0 = - 1 < 0$. Plugging this control in feedback form into the pseudo-Hamiltonian, and considering the limit conditions, we obtain the following two-points boundary value problem (BVP).

```math
\left\{
\begin{array}{l}
\dot{x}(t) = \phantom{-} \nabla_p H[t] = -x(t) + u(x(t), p(t)) = -x(t)+p(t), \\[0.5em]
\dot{p}(t) = - \nabla_x H[t] = p(t), \\[0.5em]
x(0) = x_0, \quad x(t_f) = x_f,
\end{array}
\right.
```

where $[t]~= (x(t),p(t),u(x(t), p(t)))$.

!!! note "Our goal"

Our goal is to solve this (BVP).

To achive our goal, let us first introduce the pseudo-Hamiltonian vector field

```math
\vec{H}(z,u) = \left( \nabla_p H(z,u), -\nabla_x H(z,u) \right), \quad z = (x,p),
```

and then denote by $z(\cdot,x_0,p_0)$ the solution of

```math
\dot{z}(t) = \vec{H}(z(t), u(z(t))), \quad z(0) = (x_0,p_0).
```

To compute $z$ with the `OptimalControl` package, we define the flow of the associated Hamiltonian vector field:

```@example main
u(x, p) = p
Expand All @@ -48,14 +97,33 @@ f = Flow(ocp, u)
nothing # hide
```

We define an auxiliary exponential map.
!!! note "Nota bene"

Actually, $z(\cdot, x_0, p_0)$ is also solution of

```math
\dot{z}(t) = \vec{\mathbf{H}}(z(t)), \quad z(0) = (x_0, p_0),
```
where $\mathbf{H}(z) = H(z, u(z))$ and $\vec{\mathbf{H}} = (\nabla_p \mathbf{H}, -\nabla_x \mathbf{H})$. This is what is actually computed by `Flow`.

We define also an auxiliary exponential map for clarity.

```@example main
exp(p0; saveat=[]) = f((t0, tf), x0, p0, saveat=saveat).ode_sol
nothing # hide
```

We are now in position to define the shooting function and solve the shooting equations.
Now, to solve the (BVP) we introduce the **shooting function**.

```math
\begin{array}{rlll}
S \colon & \R & \longrightarrow & \R \\
& p_0 & \longmapsto & S(p_0) = \pi(z(t_f,x_0,p_0)) - x_f,
\end{array}
```

where $\pi(x,p) = x$. At the end, solving (BVP) leads to solve $S(p_0) = 0$.
This is what we call the **indirect simple shooting method**.

```@example main
S(p0) = exp(p0)(tf)[1] - xf; # shooting function
Expand All @@ -80,53 +148,41 @@ nothing # hide
We get:

```@example main
# parameters
times = range(t0, tf, length=2)
p0min = -0.5
p0max = 2
# plot of the flow
plt_flow = plot()
p0s = range(p0min, p0max, length=20)
for i ∈ 1:length(p0s) # plot some trajectories
sol = exp(p0s[i])
x = [sol(t)[1] for t ∈ sol.t]
p = [sol(t)[2] for t ∈ sol.t]
label = i==1 ? "extremals" : false
plot!(plt_flow, x, p, color=:blue, label=label)
end
p0s = range(p0min, p0max, length=200)
xs = zeros(length(p0s), length(times))
ps = zeros(length(p0s), length(times))
for i ∈ 1:length(p0s)
sol = exp(p0s[i], saveat=times)
xs[i, :] = [z[1] for z ∈ sol.(times)]
ps[i, :] = [z[2] for z ∈ sol.(times)]
end
for j ∈ 1:length(times) # plot intermediate points
label = j==1 ? "flow at times" : false
plot!(plt_flow, xs[:, j], ps[:, j], color=:green, linewidth=2, label=label)
end
plot!(plt_flow, xlims = (-1.1, 1), ylims = (p0min, p0max))
plot!(plt_flow, [0, 0], [p0min, p0max], color=:black, xlabel="x", ylabel="p", label="x=xf")
# solution
sol = exp(p0_sol)
x = [sol(t)[1] for t ∈ sol.t]
p = [sol(t)[2] for t ∈ sol.t]
plot!(plt_flow, x, p, color=:red, linewidth=2, label="extremal solution")
plot!(plt_flow, [x[end]], [p[end]], seriestype=:scatter, color=:green, label=false) # solution
# plot of the shooting function
plt_shoot = plot(xlims=(p0min, p0max), ylims=(-2, 4), xlabel="p₀", ylabel="y")
plot!(plt_shoot, p0s, S, linewidth=2, label="S(p₀)", color=:green)
plot!(plt_shoot, [p0min, p0max], [0, 0], color=:black, label="y=0")
plot!(plt_shoot, [p0_sol, p0_sol], [-2, 0], color=:black, label="p₀ solution", linestyle=:dash)
plot!(plt_shoot, [p0_sol], [0], seriestype=:scatter, color=:green, label=false) # solution
# plots
plot(plt_flow, plt_shoot, layout=(1,2), size=(800, 450))
times = range(t0, tf, length=2) # hide
p0min = -0.5 # hide
p0max = 2 # hide
plt_flow = plot() # hide
p0s = range(p0min, p0max, length=20) # hide
for i ∈ 1:length(p0s) # hide
sol = exp(p0s[i]) # hide
x = [sol(t)[1] for t ∈ sol.t] # hide
p = [sol(t)[2] for t ∈ sol.t] # hide
label = i==1 ? "extremals" : false # hide
plot!(plt_flow, x, p, color=:blue, label=label) # hide
end # hide
p0s = range(p0min, p0max, length=200) # hide
xs = zeros(length(p0s), length(times)) # hide
ps = zeros(length(p0s), length(times)) # hide
for i ∈ 1:length(p0s) # hide
sol = exp(p0s[i], saveat=times) # hide
xs[i, :] = [z[1] for z ∈ sol.(times)] # hide
ps[i, :] = [z[2] for z ∈ sol.(times)] # hide
end # hide
for j ∈ 1:length(times) # hide
label = j==1 ? "flow at times" : false # hide
plot!(plt_flow, xs[:, j], ps[:, j], color=:green, linewidth=2, label=label) # hide
end # hide
plot!(plt_flow, xlims = (-1.1, 1), ylims = (p0min, p0max)) # hide
plot!(plt_flow, [0, 0], [p0min, p0max], color=:black, xlabel="x", ylabel="p", label="x=xf") # hide
sol = exp(p0_sol) # hide
x = [sol(t)[1] for t ∈ sol.t] # hide
p = [sol(t)[2] for t ∈ sol.t] # hide
plot!(plt_flow, x, p, color=:red, linewidth=2, label="extremal solution") # hide
plot!(plt_flow, [x[end]], [p[end]], seriestype=:scatter, color=:green, label=false) # hide
plt_shoot = plot(xlims=(p0min, p0max), ylims=(-2, 4), xlabel="p₀", ylabel="y") # hide
plot!(plt_shoot, p0s, S, linewidth=2, label="S(p₀)", color=:green) # hide
plot!(plt_shoot, [p0min, p0max], [0, 0], color=:black, label="y=0") # hide
plot!(plt_shoot, [p0_sol, p0_sol], [-2, 0], color=:black, label="p₀ solution", linestyle=:dash) # hide
plot!(plt_shoot, [p0_sol], [0], seriestype=:scatter, color=:green, label=false) # hide
plot(plt_flow, plt_shoot, layout=(1,2), size=(800, 450)) # hide
```

0 comments on commit 2cf5dfa

Please sign in to comment.