Notebook source code:
examples/turnpike/2d/turnpike.ipynb
Run the notebook yourself on binder
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
considering that the evolution of the state \((x,y) \in \mathbb{R}^2\) is governed by the control system
fixing the initial condition to
and considering that the system has to reach the final condition given by
The turnpike property in the singular flow: numerical evidence¶
In this part, we fix the control bound
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:
where \(\alpha\) is the regularization parameter. We set \(\alpha = 10^{-3}\).
[1]:
using JuMP # NLP modelling
using Ipopt # NLP solver
using Gadfly # Plots
[6]:
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[1] - c0 == 0)
@constraint(turnpike, con_x0, x[1] - x0 == 0)
@constraint(turnpike, con_y0, y[1] - 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[1] > 0) # We prefer pc[1] < 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;
[7]:
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;
[8]:
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
[9]:
# 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¶
[ ]: