Towards a package

[ocots]

We could imagine to make a package to solve optimal control problems with regions of loss control.

First, we need to automatise the direct part. This may be done in several steps. First, we need a modelisation of the optimal control problem with loss control region. Then, we need an algorithm with transform the problem into a classical optimal control problem that can be solved by OptimalControl.jl. Note that the modelisation may be helped also by OptimalControl.jl. At the end the resolution is done by OptimalControl.jl but there will be maybe some parameters to tune or homotopic parameters for instance for the regularisation. The $\varepsilon$ parameter may be fixed in the algorithm or updated.

We can use this issue to share ideas.

[ocots]

@AnasXbouali Can you write here the general definition of an optimal control problem with loss control regions?

[AnasXbouali]

Given a state space partition:

$$\mathbb{R}^n = \bigcup_{j=1}^N \overline{X}_j,$$

to which we associate the following indexation

$$q_j := \left\{\begin{array}{l} 1 \text{ if } X_j \text{ is a control region,} \\\ 0 \text{ if } X_j \text{ is a loss control region}. \end{array} \right.$$

Now, we define the following OCP with loss control regions

$$\begin{array}{lcl} \text{minimize}& & \phi(x(0),x(T)),\\[4pt] \text{subject to}& & (x,u) \in \mathrm{AC0}([0,T],\mathbb{R}^n) \times \mathrm{L}^\infty([0,T],\mathbb{R}^m), \\[2pt] & & \dot{x}(t) = f(x(t),u(t)), \quad \text{a.e.\ } t\in [0,T], \\[2pt] & & g(x(0),x(T)) \in \mathrm{S}, \\[2pt] & & u(t) \in \mathrm{U} , \quad \text{a.e.\ } t\in [0,T], \\[2pt] & & u \text{ is constant when } x \text{ is in a loss control region}. \end{array}$$

[AnasXbouali]

@ocots As we did in JDE, we need to:

• Consider a dynamics regularization for a given partition using smooth.jl.
• Add a regionally switching parameter and an auxiliary control v to model the potential changes in constant values within loss control regions. This approach also effectively handles situations where a trajectory visits a loss control region multiple times.
• Include a penalization of both the auxiliary control v and the control u, which may not be active at all times.

A first step would be to incorporate these three aspects?

[ocots]

I think that we should write in this issue how a user would interact with the package to solve this example. The transformation to a classical problem is internal and the user won't see it.

The key point is how to define the loss control regions.

First attempt:

tf  = 8
ocp = @def begin

    t ∈ [ 0, tf ], time
    x ∈ R^2, state
    u ∈ R, control

    0.5 ≤  x₂(t)  ≤ 3.5, loss_region

    x₁(0) == 0
    x₂(0) == 0
    x₂(tf) == 4

    -π/2  ≤  u(t)  ≤ π/2

    ẋ(t) == [ x₂(t) + cos(u(t)), sin(u(t)) ]

    -x₁(tf) → min

end

sol = solve(ocp; ε₀ = 1e-3, ε = 1e-3)
plot(sol)

The names of the parameters ε₀ and ε are not good. We should find explicit ones.

[ocots]

We should write some examples of definitions of loss regions. Below, each line is a region:

2 ≤ x₂(t) ≤ 3, loss_region
0 ≤ x₁(t)^2 + x₂(t)^2 ≤ 1, loss_region
0 ≤ x₁(t)^2 ≤ 1 and 2 ≤ x₂(t) ≤ 3, loss_region

# or something like that which is less usual with current version of OptimalControl
( x₁(t), x₂(t) ) ∈ [0, 1] × [2, 3], loss_region

[ocots]

Important remark

Since it is not so simple to parse the code which defines the loss control regions, we can in a first attempt make something mixing abstract definition and functional definition. Note that the abstract code call the functional one so the functional code has to be written first anyway.

tf  = 8
ocp = @def begin

    t ∈ [ 0, tf ], time
    x ∈ R^2, state
    u ∈ R, control

    x₁(0) == 0
    x₂(0) == 0
    x₂(tf) == 4

    -π/2  ≤  u(t)  ≤ π/2

    ẋ(t) == [ x₂(t) + cos(u(t)), sin(u(t)) ]

    -x₁(tf) → min

end

loss_regions!(ocp; rg = 2, lb = 2, ub = 3)

sol = solve(ocp; ε₀ = 1e-3, ε = 1e-3)

plot(sol)

Voir constraint!.

[ocots]

Let split the methods:

# definition of a classical ocp
tf  = 8
ocp = @def begin

    t ∈ [ 0, tf ], time
    x ∈ R^2, state
    u ∈ R, control

    x₁(0) == 0
    x₂(0) == 0
    x₂(tf) == 4

    -π/2  ≤  u(t)  ≤ π/2

    ẋ(t) == [ x₂(t) + cos(u(t)), sin(u(t)) ]

    -x₁(tf) → min

end

# switch to problem with loss control
loss_ocp = LossControlModel(ocp)

# add loss control regions
loss_region!(loss_ocp; rg = 2, lb = 2, ub = 3)

# convert loss control problem to a classical one
ocp = OptimalControlModel(loss_ocp; ε₀ = 1e-3, ε = 1e-3))

# solve the classical problem
sol = solve(ocp)

# plot the classical solution
plot(sol)

# or plot a solution from a loss control problem
sol_loss = LossControlSolution(sol, ocp_loss)
plot(sol_loss)

[ocots]

In the previous comment, LossControlModel is used to get a struct that can save loss control regions. The function loss_regions permits to add loss control regions. The function OptimalControlModel is the method which regularize the problem, add the variable, etc.

[ocots]

En fait, ce serait mieux de rester au niveau abstrait. Il faudrait pouvoir passer de

tf  = 8
ocp = @def begin

    t ∈ [ 0, tf ], time
    x ∈ R^2, state
    u ∈ R, control

    0.5 ≤  x₂(t)  ≤ 3.5, loss_region

    x₁(0) == 0
    x₂(0) == 0
    x₂(tf) == 4

    -π/2  ≤  u(t)  ≤ π/2

    ∂x₁ = x₂(t) + cos(u(t))
    ∂x₂ = sin(u(t))

    ẋ(t) == [ ∂x₁, ∂x₂ ]

    -x₁(tf) → min

end

à la version classique

ε₀ = 1e-3
ε   = 1e-3
tf  = 8
ocp = @def begin

    t ∈ [ 0, tf ], time
    q = [ x₁, x₂, λ ] ∈ R^3, state
    ω = [u, v] ∈ R^2, control

    x₁(0) == 0
    x₂(0) == 0
    x₂(tf) == 4

    -π/2  ≤  u(t)  ≤ π/2
    -π/2  ≤  λ(t)  ≤ π/2

    ∂x₁ = x₂(t) + cos(u(t))
    ∂x₂ = sin(u(t))

    ∂x₁_ = x₂(t) + cos(λ(t))
    ∂x₂_ = sin(λ(t))

    I = Ind(t, q(t), ε₀)

    q̇(t) == [I*∂x₁_ + (1-I)*∂x₁, I*∂x₂_ + (1-I)*∂x₂, (1-I)*v(t)]

    -x1(tf) + ∫( ε*v(t)^2 + I*u(t)^2 )→ min
end

[AnasXbouali]

D'accord, je vais essayer de travailler sur cette transformation pour l'exemple de Zermelo 1. Ensuite, nous verrons les aspects qui peuvent être automatisés.

[ocots]

Ce serait en effet super de passer directement par l'abstrait mais cela me semble assez difficile. Je reviens donc sur ce que j'ai dit et il me semble que suivre dans un premier temps ceci #14 (comment) est plus simple. Un compromis serait la chose suivante.

On part de :

tf  = 8
ocp = @def begin

    t ∈ [ 0, tf ], time
    x ∈ R^2, state
    u ∈ R, control

    0.5 ≤  x₂(t)  ≤ 3.5, loss_region

    x₁(0) == 0
    x₂(0) == 0
    x₂(tf) == 4

    -π/2  ≤  u(t)  ≤ π/2

    ∂x₁ = x₂(t) + cos(u(t))
    ∂x₂ = sin(u(t))

    ẋ(t) == [ ∂x₁, ∂x₂ ]

    -x₁(tf) → min

end

• On parse tout d'abord avec @def de CTBase.jl en enlevant au préalable toutes les lignes contenant l'étiquette loss_region. Voir ici pour un exemple de parsing. Détail dans l'issue [Dev] First parsing #36.
• On ajoute autant d'états qu'il y a de contrôle et on ajoute les contraintes associées. Voir les fichiers OptimalControlModel, les getters et les setters.
• On ajoute le contrôle permettant de changer de valeur entre les régions de perte de contrôle.
• On parse et parcours les régions de perte de contrôle indiquées par l'étiquette loss_region et on met à jour la fonction indicatrice régularisée. On met aussi à jour la dynamique à chaque fois que l'on ajoute une région de perte de contrôle.
• On met à jour le critère.

Note

Il faut une structure de données permettant de stocker les régions de perte de contrôle, la fonction indicatrice, l'ocp...