From optimal solution to real solution

The goal is now to present how we can approximate an optimal solution to a real control, in which the control can only take the values $-1$ or $1$.

To do this, let first define a model

using Filtration
using Plots

vlim = (0,200); mlim = (0,100)
ap = 0; bp = 1; cp = 0.16; dp = 3; Qp = 1
ar = 0.5; cr = 0.64; dr = 12; Qr = 2
p = [ap, bp, cp, dp, Qp, ar, cr, dr, Qr]

f₁(m,p) = -p[1]*m + p[2]
f₂(m,p) = -p[6]*m
g₁(m,p) = p[5]
g₂(m,p) = -p[9]
e₁(m,p) = p[3]*m + p[4]
e₂(m,p) = p[7]*m + p[8]

model = Model(e₁, e₂, f₁, f₂, g₁, g₂, p, :min; initial_guess = 15.)

synthesis = Synthesis(model, vlim, mlim)
plt_synthesis = plot(synthesis, xlabel = "Volume (v)", ylabel = "Mass (m)", legend = :topleft)

Let fix for instance $x_1(t_0) = m_0 = 40$ and $x_2(t_0) = v_0 = 5$, and compute the solution of the problem.

v0 = 5; m0 = 40.
sol = solution(synthesis, v0, vlim[2], m0)
plot(sol, label = "optimal ")

add_solution_to_synthesis!(plt_synthesis, sol, axis_arrow = :x, label = "optimal ")

As we can see, the control $u(\cdot)$ on the singular arc take a value between $-1$ and $1$, which is not really applicable on a real membrane. We propose three method to provide a real control solution. For all these method, we need to fix a number $N$ of swiching, which will be in this example fixed to $N = 5$

N = 5;
nothing; #hide

Method 1 : Control approximation

The first method, named :control, approximate the singular control $u(\cdot) = u_s$ by $-1$ and $1$ value on each subsegments.

Let assume that the optimal solution have a singular arc on $[T_1, T_2]$. This interval of time is divised by $N$ subintervals of time $\Delta_i = [t_i, t_{i+1}]$ where $i \in \{ 1, ... , N \}$ and where

\[t_i = T_1 + (i-1) * \Delta t, \quad \Delta t = \frac{T_2 - T_1}{N}\]

On each subintervals of time $\Delta_i$, the time spend with $u = +1$ with respect to the time spend with $u = -1$ is proportional to the value of the singular control $u_s$ in $[-1, 1]$. In other words, the (real) control $u(\cdot)$ is defined for all $i \in \{ 1, ... , N \}$ on the singular arc by

\[u(t) = \left \{ \begin{array}{ll} +1 & \quad \forall t \in [t_i, t_i + \delta t] \\ -1 & \quad \forall t \in [t_i + \delta t, t_{i+1}] \end{array} \right. \quad \text{ where } \quad \delta t = \Delta t \frac{u_s + 1}{2}. \]

This following code shows the results of this method.

real_sol₁ = real_solve(model, v0, vlim[2], m0, N, method = :control)
plt_solution = plot(sol, label = "optimal ")
plot!(plt_solution, real_sol₁, label = ":control ", color = :blue)

plt_synthesis = plot(synthesis, xlabel = "Volume (v)", ylabel = "Mass (m)", legend = :topleft)
add_solution_to_synthesis!(plt_synthesis, sol, axis_arrow = :x, label = "optimal ", linewidth = 2)
add_solution_to_synthesis!(plt_synthesis, real_sol₁, axis_arrow = :x, color = :blue, label = ":control ")

Method 2 : Singular state approximation

The second method, named :singular, is based on the singular arc approximation. The main idea is to subdivise the singular arc intro $N$ subsegments and to find the optimal control of structure $u_+ \circ u_- \circ u_+$ that minimise the energy and that satisfy that the final state (on $x_1$ and $x_2$) corresponds to the one of the end of the subsegments of the singular arc defined above.

Denoting $T_1$ the begining time of the singular arc, the singular control is defined for all $i \in \{ 1, ... N \}$ by

\[u(t) = \left \{ \begin{array}{ll} +1 & \quad \forall t \in [T_1 + (i-1) \delta t_3, T_1 + (i-1) \delta t_3 + \delta t_1] \\ -1 & \quad \forall t \in [T_1 + (i-1) \delta t_3 + \delta t_1, T_1 + (i-1) \delta t_3 + \delta t_2] \\ +1 & \quad \forall t \in [T_1 + (i-1) \delta t_3 + \delta t_2, T_1 + i \delta t_3] \end{array} \right.\]

where $\delta t_1$, $\delta t_2$ and $\delta t_3$ are solution of the following optimization problem

\[\left \{ \begin{array}{l} \min_{\delta t_1, \delta t_2, \delta t_3} \phi_0^{+-+}(\delta t_1, \delta t_2, \delta t_3) \\ \textrm{s.t. } \phi_1^{+-+}(\delta t_1, \delta t_2, \delta t_3) = x_s \\ \phantom{\textrm{s.t. }} \phi_2^{+-+}(\delta t_1, \delta t_2, \delta t_3) = \delta x_2 \end{array}\right.\]

where $\phi_k^{+-+}$ corresponds to the value of the final state $k$ ($0$ for the cost) by applying the control $+1$ on $[0, \delta t_1]$, the control $-1$ on $[\delta t_1, \delta t_2]$ and the control $+1$ on $[\delta t_2, \delta t_3]$. The value $\delta x_2$ is computed by the subdivision of the singular arc and $x_s$ corresponds to the value of the singular state.

Here is the result for this method

real_sol₂ = real_solve(model, v0, vlim[2], m0, N, method = :singular)

plt_solution = plot(sol, label = "optimal ")
plot!(plt_solution, real_sol₂, label = ":singular ", color = :red)

plt_synthesis = plot(synthesis, xlabel = "Volume (v)", ylabel = "Mass (m)", legend = :topleft)
add_solution_to_synthesis!(plt_synthesis, sol, axis_arrow = :x, label = "optimal ", linewidth = 2)
add_solution_to_synthesis!(plt_synthesis, real_sol₂, axis_arrow = :x, color = :red, label = ":singular ")

Method 3 : Optimal real solution

The third method computes the real optimal real solution of the problem. Here we define the large optimization problem which consists to find the values of the switching time which minimize the cost and satisfy the boundary conditions.

Here the code to obtain the optimal real solution

real_sol₃ = real_solve(model, v0, vlim[2], m0, N, method = :optimal)

plt_solution = plot(sol, label = "optimal ")
plot!(plt_solution, real_sol₃, label = ":optimal ", color = :green)

plt_synthesis = plot(synthesis, xlabel = "Volume (v)", ylabel = "Mass (m)", legend = :topleft)
add_solution_to_synthesis!(plt_synthesis, sol, axis_arrow = :x, label = "optimal ", linewidth = 2)
add_solution_to_synthesis!(plt_synthesis, real_sol₃, axis_arrow = :x, color = :green, label = ":optimal ")

Comparison

We compare these three methods in terms of optimality and computation time.