Notebook source code: examples/turnpike/2d/turnpike.ipynb
Run the notebook yourself on binder Binder badge

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}\).

[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)
x+y 0.0 0.5 1.0 0.0 0.5 1.0 x x+y 0.0 0.5 1.0 -1 0 1 2 3 py t 0.0 0.5 1.0 -1 0 1 2 3 py t 0.0 0.5 1.0 -2 -1 0 1 2 3 px t 0.0 0.5 1.0 0.0 0.5 1.0 1.5 c t 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 u t 0.0 0.5 1.0 -0.5 0.0 0.5 1.0 y t 0.0 0.5 1.0 0.0 0.5 1.0 x

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