# Turnpike property and the singular flow¶

Author: O. Cots

Date: 10/07/2021

We consider a two dimensional optimal control problem in Lagrange form with an affine control system with respect to the (scalar) control. The control problem consists in minimizing the cost functional

$\frac{1}{2} \int_0^{T} \left( x^2(t)+y^2(t) \right) \,\mathrm{d}t \to \min,$

considering that the evolution of the state $$(x,y) \in \mathbb{R}^2$$ is governed by the control system

$\dot{x}(t) = u(t), \quad \dot{y}(t) = - u(t) - \frac{y(t)}{2}, \quad |u(t)| \leq M,$

fixing the initial condition to

$(x(0), y(0)) = (1, -1/2)$

and considering that the system has to reach the final condition given by

$(x(T), y(T)) = (1/2, 1/2).$

## The turnpike property in the singular flow: numerical evidence¶

In this part, we fix the control bound

$M = 1$
and solve the control problem for different values of the final time $$T$$. We use the modeling language Jump in Julia with the Ipopt solver to formulate a control problem as a finite dimensional optimization problem and solve it. We use the
Gadfly package to plot the solutions.

Remark. The optimal control law will contain bang arcs, where the control is given by $$u=\pm M$$ and singular arcs where the control satisfies $$u \in (-M, M)$$. To eliminate numerical chattering due to the discretization, we add a L2-norm regularization:

$\frac{1}{2} \int_0^{T} \left( x^2(t)+y^2(t)+ \alpha u^2(t) \right) \,\mathrm{d}t \to \min,$

where $$\alpha$$ is the regularization parameter. We set $$\alpha = 10^{-3}$$.

:

using JuMP   # NLP modelling
using Ipopt  # NLP solver

:

function ocpsolve(T, N, M)
#
# This function solves the OCP:
#
# min 1/2 int_0^T (x^2(t) + y^2(t)) dt
#
#     dx = u,         -M <= u <= M
#     dy = -u-y/2
#
#     x(0) = 1, x(T) = 1/2
#     y(0) = -1/2, y(T) = 1/2
#
# Remark: a L2-norm regularization with a factor of 1e-3
# is used to eliminate chattering along the singular arc.
#
# Inputs:
#
#    - T : final time
#    - N : discretization grid size
#    - M : control bound
#
# Outputs:
#
#    - t : times t = [t_1, ..., t_N]
#    - c : costs c = [c_1, ..., c_N]
#    - x : state x = [x_1, ..., x_N]
#    - y : state y = [y_1, ..., y_N]
#    - u : control u = [u_1, ..., u_N]
#    - pc : dual of the costs pc = [pc_1, ..., pc_{N-1}]
#    - px : costate px = [px_1, ..., px_{N-1}]
#    - py : costate py = [py_1, ..., py_{N-1}]
#
# Remark: the costate px, resp. py, is the Lagrange multiplier
# associated to the differential constraint dx=u, resp. dy=-u-y/2.
#

# Create JuMP model, using Ipopt as the solver
turnpike = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)) # or 5

# Parameters
tf = T    # final time
c0 = 0    # Initial cost
x0 = 1    # Initial position
xf = 0.5  # Final position
y0 = 0    # Initial position
yf = 0    # Final position
L  = 0.5  # Parameter in the dynamics
α  = 1e-3 # L2-norm regularization factor

Δt = tf/(N-1)  # Time step

# State and control variables
@variables(turnpike, begin
c[1:N]                # Cost
x[1:N]                # State 1
y[1:N]                # State 2
-M ≤ u[1:N] ≤ M       # Control
end)

# Objective
@objective(turnpike, Min, c[N])

# Initial and final conditions
@constraint(turnpike, con_c0, c - c0 == 0)
@constraint(turnpike, con_x0, x - x0 == 0)
@constraint(turnpike, con_y0, y - y0 == 0)
@constraint(turnpike, con_xf, x[N] - xf == 0)
@constraint(turnpike, con_yf, y[N] - yf == 0)

# Dynamics
@NLexpression(turnpike, dc[j = 1:N],  0.5*(y[j]^2+x[j]^2+α*u[j]^2))
@NLexpression(turnpike, dx[j = 1:N],  u[j])
@NLexpression(turnpike, dy[j = 1:N], -u[j]-L*y[j])

# Discretization of the dynamics
# See https://jump.dev/JuMP.jl/stable/manual/constraints/#constraint_arrays
# to define an array of constraints and name them
@NLconstraint(turnpike, con_c[j = 1:N-1],
c[j+1] == c[j] + 0.5 * Δt * (dc[j+1] + dc[j]))

@NLconstraint(turnpike, con_x[j = 1:N-1],
x[j+1] == x[j] + 0.5 * Δt * (dx[j+1] + dx[j]))

@NLconstraint(turnpike, con_y[j = 1:N-1],
y[j+1] == y[j] + 0.5 * Δt * (dy[j+1] + dy[j]))

# Solve
println("Solving for N = ", N, ", T = ", T, " and M = ", M, ".")
println()

#
status = optimize!(turnpike)

if termination_status(turnpike) == MOI.OPTIMAL
println("  Solution is optimal")
elseif  termination_status(turnpike) == MOI.LOCALLY_SOLVED
println("  (Local) solution found")
elseif termination_status(turnpike) == MOI.TIME_LIMIT && has_values(turnpike)
println("  Solution is suboptimal due to a time limit, but a primal solution is available")
else
error("  The model was not solved correctly.")
end
println("  objective value = ", objective_value(turnpike))
println()

# Outputs
c  = value.(c)[:]
x  = value.(x)[:]
y  = value.(y)[:]
u  = value.(u)[:]
t  = (1:N) * value.(Δt)

p_c = zeros(N)
p_x = zeros(N)
p_y = zeros(N)
for j in 1:N-1
p_c[j] = dual(con_c[j])
p_x[j] = dual(con_x[j])
p_y[j] = dual(con_y[j])
end

if (p_c > 0) # We prefer pc < 0 since we use later the Pontryagin Maximum Principle
p_c = -p_c; p_x = -p_x; p_y = -p_y
end

pc0 = dual(con_c0)
px0 = dual(con_x0)
pxf = dual(con_xf)
py0 = dual(con_y0)
pyf = dual(con_yf)

if(pc0*p_c[  1]<0); pc0 = -pc0; end
if(px0*p_x[  1]<0); px0 = -px0; end
if(pxf*p_x[N-1]<0); pxf = -pxf; end
if(py0*p_y[  1]<0); py0 = -py0; end
if(pyf*p_y[N-1]<0); pyf = -pyf; end

p_c[1:N] = [pc0; p_c[1:N-1]]
p_x[1:N] = [px0; (p_x[1:N-2]+p_x[2:N-1])/2; pxf] # We add the multiplier from the limit conditions
p_y[1:N] = [py0; (p_y[1:N-2]+p_y[2:N-1])/2; pyf]

return (t, c, x, y, u, p_c, p_x, p_y)
end;

:

mutable struct sollayers
x; y; c; u; px; py; rs; lp; lx;
end

function getlayers(t, c, x, y, u, p_c, p_x, p_y, col)

N = size(t, 1)
x_layer  = layer(x = t/t[N], y = x,   Geom.path, Theme(default_color=color(col)))
y_layer  = layer(x = t/t[N], y = y,   Geom.path, Theme(default_color=color(col)))
c_layer  = layer(x = t/t[N], y = c,   Geom.path, Theme(default_color=color(col)))
u_layer  = layer(x = t/t[N], y = u,   Geom.path, Theme(default_color=color(col)))
px_layer = layer(x = t/t[N], y = p_x, Geom.path, Theme(default_color=color(col)))
py_layer = layer(x = t/t[N], y = p_y, Geom.path, Theme(default_color=color(col)))

rs_layer = layer(x = x,   y = x-y,   Geom.path, Theme(default_color=color(col)))
lp_layer = layer(x = x+y, y = p_y,   Geom.path, Theme(default_color=color(col)))
lx_layer = layer(x = x+y, y = x,     Geom.path, Theme(default_color=color(col)))

layers = sollayers( x_layer, y_layer, c_layer, u_layer, px_layer, py_layer,  rs_layer, lp_layer, lx_layer)

return layers
end

function plotsolution(args...)
xp = (); yp = (); cp = (); up = (); pxp = (); pyp = (); rsp = (); lpp = (); lxp = ()
for layer in args
xp = (xp..., layer.x); yp = (yp..., layer.y); cp = (cp..., layer.c)
up = (up..., layer.u); pxp = (pxp..., layer.px); pyp = (pyp..., layer.py)
rsp = (rsp..., layer.rs); lpp = (lpp..., layer.lp); lxp = (lxp..., layer.lx)
end

x_plot  = plot(xp...,  Guide.xlabel("t"), Guide.ylabel("x"))
y_plot  = plot(yp...,  Guide.xlabel("t"), Guide.ylabel("y"))
c_plot  = plot(cp...,  Guide.xlabel("t"), Guide.ylabel("c"))
u_plot  = plot(up...,  Guide.xlabel("t"), Guide.ylabel("u"))
px_plot = plot(pxp..., Guide.xlabel("t"), Guide.ylabel("px"))
py_plot = plot(pyp..., Guide.xlabel("t"), Guide.ylabel("py"))
rs_plot = plot(rsp..., Guide.xlabel("x"), Guide.ylabel("x-y"))
lp_plot = plot(lpp..., Guide.xlabel("x+y"), Guide.ylabel("py"))
lx_plot = plot(lxp..., Guide.xlabel("x+y"), Guide.ylabel("x"))

draw(SVG(8inch, 16inch), vstack(hstack(x_plot, y_plot),
hstack(u_plot, c_plot),
hstack(px_plot, py_plot),
hstack(lp_plot, lx_plot)))
#hstack(rs_plot, lp_plot, lx_plot)))
end;

:

M = 1   # control bound

N = 500 # size grid
T = 10  # final time
t, c, x, y, u, p_c, p_x, p_y = ocpsolve(T, N, M);
l1 = getlayers(t, c, x, y, u, p_c, p_x, p_y, "blue");

N = 500 # size grid
T = 20  # final time
t, c, x, y, u, p_c, p_x, p_y = ocpsolve(T, N, M);
l2 = getlayers(t, c, x, y, u, p_c, p_x, p_y, "orange");

N = 2000 # size grid
T = 50   # final time
t, c, x, y, u, p_c, p_x, p_y = ocpsolve(T, N, M);
l3 = getlayers(t, c, x, y, u, p_c, p_x, p_y, "red");

Solving for N = 500, T = 10 and M = 1.

(Local) solution found
objective value = 1.0581606367508871

Solving for N = 500, T = 20 and M = 1.

(Local) solution found
objective value = 1.1503582794529834

Solving for N = 2000, T = 50 and M = 1.

(Local) solution found
objective value = 1.1530849771283864


:

# We plot the state, costate, control and cost wrt time, for different values of T.
#
# Blue:   l1 (T = 10)
# Orange: l2 (T = 20)
# Red:    l3 (T = 50)
#
plotsolution(l1, l2, l3)


On can notice that for these final times, the optimal control is Bang-Singular-Bang. Besides, along the singular arc, one can observe the turnpike phenomenon since when the final time becomes larger, all the state, costate and control curves stay close to the point $$x=y=p_x=p_y=u=0$$. In the rest of this notebook, we will exhibit this turnpike phenomenon inside the singular flow and we will show that actually, in this control problem, one encounter two situations of singular perturbations, due to the parameters $$M$$ and $$T$$, resp. the control bound and the final time.

Remark. This property of turnpike inside the singular flow has already been observed in the article microbial where a singular arc of order 2 plays a crucial role.

## Pontryagin Maximum Principle and the Turnpike phenomenon in the singular flow¶

[ ]:




thumbnail