Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
name = "ControlLoss"
uuid = "6a3a0e7e-9c90-4c26-bc64-8dd630726a7b"
authors = ["Olivier Cots <olivier.cots@irit.fr>"]
version = "0.1.1"
version = "0.1.3"
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -17,4 +17,4 @@ The control-toolbox ecosystem gathers Julia packages for mathematical control an

| **Documentation** | **Code Status** |
|:-------------------|:-----------------|
| [![Documentation][doc-stable-img]][doc-stable-url] [![Documentation][doc-dev-img]][doc-dev-url] | [![Build Status][ci-img]][ci-url] [![Covering Status][co-img]][co-url] |
| [![Documentation][doc-stable-img]][doc-stable-url] [![Documentation][doc-dev-img]][doc-dev-url] | [![Build Status][ci-img]][ci-url] [![Covering Status][co-img]][co-url] |
4 changes: 3 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,16 @@ makedocs(
sitename = "Control loss",
format = Documenter.HTML(
prettyurls = false,
size_threshold_ignore = ["zermelo.md", "ho.md"],
size_threshold_ignore = ["statement.md", "numerical.md","zermelo.md", "ho.md"],
assets=[
asset("https://control-toolbox.org/assets/css/documentation.css"),
asset("https://control-toolbox.org/assets/js/documentation.js"),
],
),
pages = [
"Introduction" => "index.md",
"Statement of the problem" => "statement.md",
"Numerical approach" => "numerical.md",
"Zermelo navigation problem" => "zermelo.md",
"Harmonic oscillator problem" => "ho.md",

Expand Down
74 changes: 65 additions & 9 deletions docs/src/ho.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,3 @@



## Harmonic oscillator problem

```math
Expand Down Expand Up @@ -41,6 +38,7 @@
using OptimalControl
using NLPModelsIpopt
include("smooth.jl");
nothing # hide
```


Expand All @@ -53,7 +51,7 @@ plot(fNC,-1., 1, label="fNC")


```@example main
@def ocp begin
@def ocp begin
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the end of the code, you can add nothing # hide if you want.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, I used nothing # hide for most of the code and it worked well except in

N = 500
sol = solve(ocp; grid_size=N); 
nothing # hide

where it (still) shows the output.

Copy link
Copy Markdown
Member

@ocots ocots Aug 20, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@AnasXbouali This is not really an output but a trace during the optimization process. To turn it off you can use the option display:

N = 500
sol = solve(ocp; grid_size=N, display=false)

and if you don't want to print the output, here print sol, then you can add:

N = 500
sol = solve(ocp; grid_size=N, display=false)
nothing # hide

Note. You can remove ";" when using nothing # hide.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Noted, Thanks!

Copy link
Copy Markdown
Member

@ocots ocots Aug 20, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@AnasXbouali If you want the trace but there are many iterations and you want a shorter display, you can use the print_level option from Ipopt.

sol = solve(ocp; grid_size=N, print_level=4)

See https://control-toolbox.org/OptimalControl.jl/stable/tutorial-double-integrator.html#Solve-and-plot


ε = 1e-3

Expand All @@ -74,8 +72,8 @@ plot(fNC,-1., 1, label="fNC")
xv(0) == 0

#final condition
x1(tf) == 1e-6
x2(tf) == -1e-6
x1(tf) == 0
x2(tf) == 0

#control constraint
-1. ≤ u(t) ≤ 1.
Expand All @@ -90,15 +88,73 @@ plot(fNC,-1., 1, label="fNC")

#cost function
tf + ε*xv(tf) + xu(tf) → min
end
end;
nothing # hide
```


```@example main
N = 400
sol = solve(ocp; init = (state = t -> [0.1, 0.1, 1, 0, 0], control =[-1, 0], variable =15), grid_size=N);
nothing # hide
```

```@example main
plot(sol; layout=:group, size=(800, 300))
```

```@example main
tt = sol.times
tf = tt[end]
x1(t) = sol.state(t)[1]
x2(t) = sol.state(t)[2]
λ(t) = sol.state(t)[3]
u(t) = sol.control(t)[1];
nothing # hide
```

```@example main
sol = solve(ocp; init = (state = t -> [0.1, 0.1, 1, 0, 0], control =[-1, 0], variable =15), max_iter=650)
# Plot the optimal trajectory
plot(x1, x2, 0, tf, label="optimal trajectory", color="blue", linewidth=2)
plot!([-4, 4], [0, 0], color=:black, label=false, linewidth=2)
```

```@example main
plot(tt, u, label="optimal control", color="red", linewidth=2)
plot!(tt, λ, label="state λ", color="green", linewidth=2)
```

```@example main
plot(sol; layout=:group, size=(800, 300))
# Find the crossing times based on conditions for x1
t1_index = findfirst(t -> x2(t) ≤ 0, tt)
t2_index = nothing
t3_index = nothing
t4_index = nothing

# If t1 is found, find the next crossing times
if t1_index !== nothing
t2_index = findfirst(t -> x2(t) ≥ 0, tt[t1_index+1:end])
t2_index = t2_index !== nothing ? t2_index + t1_index : nothing
end

if t2_index !== nothing
t3_index = findfirst(t -> x2(t) ≤ 0, tt[t2_index+1:end])
t3_index = t3_index !== nothing ? t3_index + t2_index : nothing
end

if t3_index !== nothing
t4_index = findfirst(t -> x2(t) ≥ 0, tt[t3_index+1:end])
t4_index = t4_index !== nothing ? t4_index + t3_index : nothing
end

# Convert indices to times
t11 = t1_index !== nothing ? tt[t1_index] : "No such t1 found"
t22 = t2_index !== nothing ? tt[t2_index] : "No such t2 found"
t33 = t3_index !== nothing ? tt[t3_index] : "No such t3 found"
t44 = t4_index !== nothing ? tt[t4_index] : "No such t4 found"

println("First crossing time: ", t11)
println("Second crossing time: ", t22)
println("Third crossing time: ", t33)
println("Fourth crossing time: ", t44)
```
42 changes: 19 additions & 23 deletions docs/src/index.md
Original file line number Diff line number Diff line change
@@ -1,39 +1,35 @@
# Loss control regions in optimal control problems

## Introduction

**General context.** *Optimal control theory* studies controlled systems to achieve desired targets with minimal cost. The *Pontryagin maximum principle* (PMP) provides necessary conditions for optimality, ensuring an *adjoint vector* (or *costate*) meets the *Hamiltonian maximization condition*.
**General context.** *Optimal control theory* studies controlled systems to achieve desired targets with minimal cost. The *Pontryagin maximum principle* (PMP, in short) provides necessary conditions for optimality, ensuring an *adjoint vector* (or *costate*) meets the *Hamiltonian maximization condition*[^1].

Typically, optimal control involves *permanent control*, allowing modification of the control function at each time instant. However, practical constraints can lead to *nonpermanent control*. For instance, digital controls result in *sampled-data control* with discrete changes. In aerospace, *eclipse constraints* limit control for solar-powered satellites in a shadow region where the control is reduced to zero. Hence, it is desirable to keep the system outside these regions.
Typically, optimal control involves *permanent control*, allowing modification of the control function at each time instant. However, practical constraints can lead to *nonpermanent control*. For instance, digital controls result in *sampled-data control* with discrete changes[^3],[^4]. In aerospace, *eclipse constraints* limit control for solar-powered satellites in a shadow region where the control is reduced to zero[^5]. Hence, it is desirable to keep the system outside these regions.

![aerospace](resources/aerospace.jpg)
![aerospace](resources/aerospace.jpg)

```math
\begin{equation}
\begin{aligned}
&\text{minimize} && -x_1(8), \\[10pt]
&\text{subject to} && (x,u) \in \mathrm{AC}([0,8],\mathbb{R}^2) \times \mathrm{L}^\infty([0,8],\mathbb{R}), \\[6pt]
& && \dot{x}_1(t) = x_2(t) + \cos(u(t)), \quad \text{a.e. } t \in [0,8], \\[6pt]
& && \dot{x}_2(t) = \sin(u(t)), \quad \text{a.e. } t \in [0,8], \\[6pt]
& && x(0)=0_{\mathbb{R}^2}, \quad x_2(8)=4, \\[6pt]
& && u(t) \in \left[-\frac{\pi}{2},\frac{\pi}{2}\right], \quad \text{a.e. } t \in [0,8], \\[6pt]
& && u \text{ is constant when } x \text{ is in the loss control region } \{x \in \mathbb{R}^2 \mid 0.5 < x_2 < 3.5 \}.
\end{aligned}
\end{equation}
```

**Objective and approach.** Here, we address optimal control problems with *loss control regions*[^1], where the state space is divided into *control regions* and *loss control regions*. In control regions, control can change at any time, while in loss control regions, control must remain constant, though its value is to be optimized and can vary with each visit.
**Objective and approach.** Here, we address optimal control problems with *loss control regions* (we refer to [^2] for further details), where the state space is divided into *control regions* and *loss control regions*. In control regions, control can change at any time, while in loss control regions, control must remain constant, though its value is to be optimized and can vary with each visit.

We extend our previous work by using a permanent control for control regions and a *regionally switching parameter* for loss control regions. This leads to a discontinuous dynamics framework, fitting into *spatially hybrid optimal control theory*. The *hybrid maximum principle* (HMP) extends the PMP to hybrid settings, with a piecewise absolutely continuous adjoint vector.
We extend our previous work by using a permanent control for control regions and a *regionally switching parameter* for loss control regions. This leads to a discontinuous dynamics framework, fitting into *optimal control problems involving spatially heterogenous dynamics*[^2],[^5]. The *hybrid maximum principle* (HMP, in short) extends the PMP to hybrid settings[^2],[^5], with a piecewise absolutely continuous adjoint vector.

**Numerical contribution.** In this note we illustrate a two-step numerical method for optimal control problems with loss control regions. First, a direct numerical approach is applied to a *regularized* problem to manage discontinuities and outline the optimal trajectory's structure. Second, this helps initialize an indirect numerical method for the original problem, using the PMP from Theorem~1. The method incorporates the averaged Hamiltonian gradient condition and adjoint vector discontinuities to define an appropriate shooting function, adding to classical terms for non-hybrid optimal control problems (see ref).
**Numerical contribution.** In this note we illustrate a two-step numerical method for optimal control problems with loss control regions. First, a direct numerical approach is applied to a *regularized* problem to manage discontinuities and outline the optimal trajectory's structure. Second, this helps initialize an indirect numerical method for the original problem, using the PMP with loss control regions. The method incorporates the *averaged Hamiltonian gradient condition*[^3],[^4] and adjoint vector discontinuities to define an appropriate shooting function, adding to classical terms for non-hybrid optimal control problems.

!!! note "Contents"
- Provide a statement of the PMP with loss control regions.
- Provide a direct method for solving optimal control problems with loss control regions (based on a regularization technique).
- Provide an indirect method (shooting method) for solving optimal control problems with loss control regions using the PMP with loss control regions.
- Provide a statement of the **PMP with loss control regions**.
- Provide a **direct method** for solving optimal control problems with loss control regions (based on a **regularization** technique).
- Provide an **indirect method** (shooting method) for solving optimal control problems with loss control regions using the PMP with loss control regions.


[^1]: L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, E. F. Mishchenko, *The Mathematical Theory of Optimal Processes*, A Pergamon Press Book. The Macmillan Co., New York, 1964.

[^2]: T. Bayen, A. Bouali, L. Bourdin & O. Cots, *Loss control regions in optimal control problems*, Journal of Differential Equations, **405** (2024), 359-397.

[^3]: P. Bettiol, L. Bourdin, *Pontryagin Maximum Principle for State Constrained Optimal Sampled-Data Control Problems on Time Scales*, ESAIM Control Optim. Calc. Var., **27** (2021) 51.

[^4]: L. Bourdin, G. Dhar, *Optimal Sampled-Data Controls with Running Inequality State Constraints: Pontryagin Maximum Principle and Bouncing Trajectory Phenomenon*, Mathematical Programming, **191** (2022) 907–951.

[^1]: T. Bayen, A. Bouali, L. Bourdin & O. Cots, *Loss control regions in optimal control problems*, Journal of Differential Equations, **12** (2024) 405, 359-397.
[^5]: T. Haberkorn, E. Trélat, *Convergence Results for Smooth Regularizations of Hybrid Nonlinear Optimal Control Problems*, SIAM Journal on Control and Optimization, **49** (2011) 1498–1522.

## Dependencies

Expand Down
91 changes: 91 additions & 0 deletions docs/src/numerical.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
# A numerical approach for optimal control problems with loss control regions

In optimal control theory, there are several ways for solving numerically an optimal control problem. Direct and indirect methods represent an important class of methods that we will use hereafter.

Direct methods involve discretizing the state and control variables, simplifying the problem into a finite-dimensional nonlinear optimization problem. On the other hand, indirect methods tackle the problem by solving a boundary value problem, based on the PMP, through the use of a shooting method (see, e.g.[^1],[^2]).

It is important to note that neither of these methods is fundamentally better than the other. Indeed, each of these methods has its pros and cons. For instance, although the direct method is simple to implement, more robust and less sensitive to the choice of the initial condition, it should be noted that it yields less precise results and can converge to local minima that significantly deviate from the optimal solution. Additionally, this method requires a large amount of memory. On the other hand, the indirect method is known for its extreme precision. However, it is based only on necessary optimality conditions from the PMP and often requires knowledge of the structure of the optimal solution. Moreover, it is quite sensitive to the choice of the initial condition, which must be chosen carefully to ensure convergence.

Often in the literature, one proceeds in two steps. The first step is to implement a direct method to determine the optimal solution's structure and extract the associated adjoint vector. The second step involves constructing an indirect shooting method, where the initial condition is based on the numerical results obtained from the direct method.


## Description of the direct method

For some $\omega_0 \in \mathrm{U}$, some $\varepsilon_0>0$ and $\varepsilon>0$ small enough, we introduce the *regularized problem* given by
```math
\begin{equation}
\begin{array}{lcl}
\text{minimize}& & \phi(x(0),x(T)) + \varepsilon_0 \displaystyle\int_0^T v^2(t) \, dt + \displaystyle\int_0^T (1- \Psi_{\varepsilon}(x(t)))\|u(t)-\omega_0\|_{\R^m}^2 \, dt,\\[10pt]
\text{subject to}& & (x,\lambda,u,v) \in \mathrm{AC}([0,T],\R^n) \times \mathrm{AC}([0,T],\R^m) \times \mathrm{L}^\infty([0,T],\R^m)\times \mathrm{L}^\infty([0,T],\R), \\[2pt]
& & \dot{x}(t) = \Psi_{\varepsilon}(x(t)) f(x(t),u(t)) + (1-\Psi_{\varepsilon}(x(t))) f(x(t),\lambda(t)), \quad \text{a.e.\ } t\in [0,T], \\[2pt]
& & \dot{\lambda}(t) =\Psi_{\varepsilon}(x(t))v(t), \quad \text{a.e.\ } t\in [0,T], \\[2pt]
& & g(x(0),x(T)) \in \mathrm{S}, \\[2pt]
& & \lambda(t) \in \mathrm{U}, \quad \text{a.e.\ } t\in [0,T],\\[2pt]
& & (u(t),v(t)) \in \mathrm{U} \times \R , \quad \text{a.e.\ } t\in [0,T],
\end{array}
\end{equation}
```
where $\Psi_{\epsilon} : \R^n \to \R$ is the regularization of the characteristic function of $\cup_{q_j = 1} \overline{X_j}$ given by
```math
\Psi_{\varepsilon}(x) := \sum_{q_j = 1} e^{-\frac{1}{2 \varepsilon} \mathrm{d}^2_{j}(x)},
```
for all $x\in \R^n$, where $\mathrm{d}_{j}: \R^n \to \R$ stands for the distance function to the set $\overline{X_j}$ defined by $\mathrm{d}_{j}(x) := \inf_{y \in \overline{X_j}} \|x-y\|_{\R^n}$ for all $x\in \R^n$ and every $j \in \mathcal{J}$.

## Description of the indirect method

Recall that the direct method has captured the structure of the optimal pair $(x^*,u^*)$. In the indirect method, we address each arc separately.
We begin by defining the **flow** of the Hamiltonian associated with each arc. To accomplish this, we use the function `Flow` that can be found in the [`CTFlows.jl`](https://github.com/control-toolbox/CTFlows.jl) package. This latter allows to solve the Hamiltonian system over a given time interval from given initial values of the state and the adjoint vector. This function requires necessary libraries such as \texttt{ForwardDiff} for calculating gradients and Jacobians and `OrdinaryDiffEq` for solving ordinary differential equations.

In the setting of the present paper (including control regions and loss control regions), we distinguish between two types of Hamiltonian flows:

- **Hamiltonian flows in control regions** We recall that the Hamiltonian $H$ associated with the optimal control problem given above is defined by
```math
H(x,u,p) := \langle p, f(x,u)\rangle_{\R^n},
```
for all $(x,u,p)\in \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R}^n$. Using the theorem given above and more specifically the Hamiltonian maximization condition, we obtain an expression of the control $u^*$ which can generate a sequence of arcs. Then it remains to define a \textit{pseudo-Hamiltonian}\footnote{the pseudo-Hamiltonian stands for the Hamiltonian flow associated with each arc.} associated with each arc. Finally we define the flow associated with each arc, which allows the resolution of the boundary value problem on a time interval satisfied by the pair $(x^*,p)$ with an initial condition.

- **Hamiltonian flows in loss control regions** Recall that, in loss control regions, $u^*$ satisfies an averaged Hamiltonian gradient condition (see Theorem above). Here, the difficulty lies in the fact that this condition is given in an integral and implicit form. Therefore, to overcome this difficulty, we first introduce new states $\lambda$ and $y$. First, the state $\lambda$ comes from the augmentation technique (we refer to [^3] for more details) to handle the constant value $u^*_k$, so it satisfies the dynamics $\dot{\lambda}(t)=0$ and the initial condition $\lambda(\tau^*_{k-1}) = u^*_k$. Second, the state $y$ satisfies $\dot{y}(t) = 0$ and $y(\tau^*_{k-1})=0$ (and thus $y=0$). Now we define the new Hamiltonian $\tilde{H}$ as follows:
```math
\tilde{H}(x,u,y,p)
:= H(x,u,p) - y \nabla_{u} H(x,u,p) = \langle p, f(x,u)\rangle_{\R^n} - y \nabla_{u} H(x,u,p),
```
for all $(x,u,y,p)\in \mathbb{R}^n \times \mathbb{R}^m \times \R \times \mathbb{R}^n $.
It is important to note that, since $y=0$, we recover the same Hamiltonian $H$. But, the actual utility of introducing the state $y$ is that it allows us to rewrite the integral expressed in the averaged Hamiltonian gradient condition as a terminal value of an adjoint vector. This makes it easier to take into account in the shooting function. Here is the justification of this point. First, we define $p_{y}$ as the solution to the system
```math
\left\{\begin{array}{l}
\dot{p}_{y}(t) = - \nabla_y \tilde{H}(x^*(t),u^*_k,y(t),p(t)), \quad
\text{for a.e. } t\in [\tau^*_{k-1}, \tau^*_k], \\
{p}_{y}(\tau^*_{k-1}) = 0_{\R^n}.
\end{array}
\right.
```
Second, since $u_k^*$ is assumed to be an interior value to $\mathrm{U}$, we get that
```math
\begin{equation*}
\int_{\tau^*_{k-1}}^{\tau^*_k} \nabla_{u} H(x^*(t),u^*_k,p(t)) \, dt = {p}_{y}(\tau^*_k)=0.
\end{equation*}
```
Hence there is no need to compute an integral in order to take into account the averaged Hamiltonian gradient condition.




















[^1]: S. Aronna, F. Bonnans, P. Martinon, *A shooting algorithm for optimal control problems with singular arcs*, Journal of Optimization Theory and Applications, **158**, pp. 419–459, Springer, 2013.
[^2]: J.F. Bonnans, *The shooting approach to optimal control problems*, IFAC Proceedings Volumes, **46**, no. 11, pp. 281–292, Elsevier, 2013.
[^3]: T. Bayen, A. Bouali, L. Bourdin & O. Cots, *Loss control regions in optimal control problems*, Journal of Differential Equations, **405** (2024), 359-397.
Loading