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

Bang-Singular-Bang problem and turnpike

  • Author: Olivier Cots

  • Date: November 2021

I) Description of the optimal control problem

We consider the following optimal control problem:

\[\begin{split}\left\{ \begin{array}{l} \displaystyle J(u) := \displaystyle \int_0^{t_f} x^2(t) \, \mathrm{d}t \longrightarrow \min \\[1.0em] \dot{x}(t) = f(x(t), u(t)) := \displaystyle u(t), \quad |u(t)| \le 1, \quad t \in [0, t_f] \text{ a.e.}, \\[1.0em] x(0) = 1, \quad x(t_f) = 1/2. \end{array} \right.\end{split}\]

To this optimal control problem is associated the stationnary optimization problem

\[\min_{(x, u)} \{~ x^2 ~ | ~  (x, u) \in \mathrm{R} \times [-1, 1],~ f(x,u) = u = 0\}.\]

The static solution is thus \((x^*, u^*) = (0, 0)\). This solution may be seen as the static pair \((x, u)\) which minimizes the cost \(J(u)\) under the constraint \(u \in [-1, 1]\). It is well known that this problem is what we call a turnpike optimal control problem. Hence, if the final time \(t_f\) is long enough the solution is of the following form: starting from \(x(0)=1\), reach as fast as possible the static solution, stay at the static solution as long as possible before reaching the target \(x(t_f)=1/2\). In this case, the optimal control would be

\[\begin{split}u(t) = \left\{ \begin{array}{lll} -1 & \text{if} & t \in [0, t_1], \\[0.5em] \phantom{-}0 & \text{if} & t \in (t_1, t_2], \\[0.5em] +1 & \text{if} & t \in (t_2, t_f], \end{array} \right.\end{split}\]

with \(0 < t_1 < t_2 < t_f\). We say that the control is Bang-Singular-Bang. A Bang arc corresponds to \(u \in \{-1, 1\}\) while a singular control corresponds to \(u \in (-1, 1)\). Since the optimal control law is discontinuous, then to solve this optimal control problem by indirect methods and find the switching times \(t_1\) and \(t_2\), we need to implement what we call a multiple shooting method. In the next section we introduce a regularization technique to force the control to be in the set \((-1,1)\) and to be smooth. In this context, we will be able to implement a simple shooting method and determine the structure of the optimal control law. Thanks to the simple shooting method, we will have the structure of the optimal control law together with an approximation of the switching times that we will use as initial guess for the multiple shooting method that we present in the last section.

Main goal

Find the switching times \(t_1\) and \(t_2\) by multiple shooting.

Steps:

  1. Regularize the problem and solve the regularized problem by simple shooting.

  2. Determine the structure of the non-regularized optimal control problem, that is the structure Bang-Singular-Bang, and find a good approximation of the switching times and of the initial co-vector.

  3. Solve the non-regularized optimal control problem by multiple shooting.

Remark 1. See this page for a general presentation of the simple shooting method.

Remark 2. In this particular example, the singular control does not depend on the costate \(p\) since it is constant. This happens in low dimension. This could be taken into consideration to simplify the definition of the multiple shooting method. However, to stay general, we will not consider this particular property in this notebook.

II) Regularization and simple shooting

We make the following regularization:

\[\begin{split}\left\{ \begin{array}{l} \displaystyle J(u) := \displaystyle \int_0^{t_f} (x^2(t) - \varepsilon\ln(1-u^2(t))) \, \mathrm{d}t \longrightarrow \min \\[1.0em] \dot{x}(t) = f(x(t), u(t)) := \displaystyle u(t), \quad |u(t)| \le 1, \quad t \in [0, t_f] \text{ a.e.}, \\[1.0em] x(0) = 1, \quad x(t_f) = 1/2. \end{array} \right.\end{split}\]

Our goal is to determine the structure of the optimal control problem when \((\varepsilon, t_f) = (0, 2)\). The problem is simpler to solver when \(\varepsilon\) is bigger and \(t_f\) is smaller. It is also smooth whenever \(\varepsilon>0\). Hence, we will start by solving the problem for \((\varepsilon, t_f) = (1, 1)\). In a second step, we will decrease the penalization term \(\varepsilon\) and increase the final time \(t_f\) to the final value \(2\).

This problem is still a turnpike control problem but for \(\varepsilon>0\), the optimal control is smooth since the strict Legendre-Clebsch condition will be satisfied under this assumption. For \(\varepsilon>0\), the stationnary solution is still \((x, u) = (0, 0)\) and the associated Lagrange multiplier is \(p=0\). Hence, we are going to start by solving the regularized optimal control problem with the limit conditions \(x(0)=x(t_f)=0\). This particular case amounts to solve the stationnary optimization problem for which the optimal solution is known.

Preliminaries

[1]:

# import packages
import nutopy as nt
import nutopy.tools as tools
import nutopy.ocp as ocp
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [16, 4]
plt.rcParams['figure.dpi'] = 100

[2]:

# Finite differences function
# Return f'(x).dx
def finite_diff(fun, x, dx, *args, **kwargs):
    v_eps = np.finfo(float).eps
    t = np.sqrt(v_eps) * np.sqrt(np.maximum(1.0, np.linalg.norm(x))) / np.sqrt(np.maximum(1.0, np.linalg.norm(dx)))
    j = (fun(x + t*dx, *args, **kwargs) - fun(x, *args, **kwargs)) / t
    return j

[3]:

# Compilation of the hamiltonian and associated derivatives up to order 2

!python -m numpy.f2py -c hfun.f90       -m hfun       > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d.f90     -m hfun_d     > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d_d.f90   -m hfun_d_d   > /dev/null 2>&1

from hfun       import hfun       as hf
from hfun_d     import hfun_d     as hf_d
from hfun_d_d   import hfun_d_d   as hf_d_d

[4]:

# Parameters

t0        = 0.0
x0        = np.array([1.0])
xf_target = np.array([0.5])

x_steady  = np.array([0.0+0.001]) # solution of stationnary problem + epsilon

e_init    = 1.0
e_final   = 0.001 #

tf_init   = 1.0 # With this value the problem is simpler to solver since the trajectory stay
                # less time around the static solution
tf_final  = 2.0

Maximized Hamiltonian and its derivatives

The pseudo-Hamiltonian is (in the normal case)

\[H(x,p,u,\varepsilon) = pu - x^2 + \varepsilon ln(1-u^2).\]

Note that we put the parameter \(\varepsilon\) into the arguments of the pseudo-Hamiltonian since we will vary it.

[5]:

# Fortran Hamiltonian code
!pygmentize hfun.f90

subroutine hfun(x, p, e, h)

    double precision, intent(in)  :: x(1), p(1), e
    double precision, intent(out) :: h

    double precision :: u

    if(abs(p(1))<1e-8)then
        u = 0d0
    else
        u = (-e + sqrt(e**2+p(1)**2))/p(1)
    end if
        
    h = p(1)*u - x(1)**2 + e*(log(1d0-u**2))

end subroutine hfun

[6]:

# Control in feedback form u[p,e].
#
@tools.vectorize(vvars=(1,))
def ufun(p, e):
    if(np.abs(p[0])<1e-8):
        u = 0.
    else:
        u = (-e + np.sqrt(e**2+p[0]**2))/p[0]
    return u

We give next the maximized Hamiltonian with its derivatives. This permits us to define the flow of the associated Hamiltonian vector field.

[7]:

# Hamiltonian and derivatives
hfun   = lambda t, x, p, e                : hf(x, p, e)
dhfun  = lambda t, x, dx, p, dp, e, de    : hf_d(x, dx, p, dp, e, de)
d2hfun = lambda t, x, dx, d2x, p, dp, d2p, e, de, d2e : hf_d_d(x, d2x, dx, p, d2p, dp, e, d2e, de)

hfun   = nt.tools.tensorize(dhfun, d2hfun, tvars=(2, 3, 4), full=True)(hfun)
h      = nt.ocp.Hamiltonian(hfun)

f = ocp.Flow(h)             # The flow associated to the Hamiltonian object is
                            # the exponential mapping with its derivative
                            # that can be used to define the Jacobian of the
                            # shooting function

Shooting function and its derivative

The shooting function is

\[\begin{split}S(x_\mathrm{mid}, p_\mathrm{mid}) := \begin{pmatrix} x(t_0, (t_0+t_f)/2, x_\mathrm{mid}, p_\mathrm{mid}, \varepsilon) - x_0 \\ x(t_f, (t_0+t_f)/2, x_\mathrm{mid}, p_\mathrm{mid}, \varepsilon) - x_f \end{pmatrix},\end{split}\]

where \(t_f = 1\), \(\varepsilon = 1\), \(x_0=x^*\), \(x_f=x^*\), \(x^* = 0\) (the steady state) and where \(z(t_f, t_0, x_0, p_0, \varepsilon)\) is the solution of the associated Hamiltonian system with the initial condition \(z(t_0) = (x_0, p_0)\).

Procedure

First solve \(S=0\) then use homotopy techniques to converge towards the solution of the original problem.

Remark: This first shooting is useless since we know the solution. We do it for illustration and validation.

[8]:

# Definition of the shooting function
#
def shoot(zmid):
    s      = np.zeros((2,))
    x12    = zmid[0:1]
    p12    = zmid[1:2]
    xf, pf = f((t0+tf_init)/2., x12, p12, tf_init, e_init)
    xi, pi = f((t0+tf_init)/2., x12, p12,      t0, e_init)
    s[0]   = xf[0] - x_steady[0]
    s[1]   = xi[0] - x_steady[0]
    return s

[9]:

# Derivative of the shooting function
#
def dshoot(zmid, dzmid):
    ds     = np.zeros((2,))
    x12    = zmid[0:1]
    p12    = zmid[1:2]
    dx12   = dzmid[0:1]
    dp12   = dzmid[1:2]
    (xf, dxf), _ = f((t0+tf_init)/2., (x12, dx12), (p12, dp12), tf_init, e_init)
    (xi, dxi), _ = f((t0+tf_init)/2., (x12, dx12), (p12, dp12),      t0, e_init)
    ds[0]  = dxf[0]
    ds[1]  = dxi[0]
    return ds

shoot = nt.tools.tensorize(dshoot, tvars=(1,))(shoot)

[10]:

# Function to plot the solution
def plotSolution(zmid, e, tf):

    x12     = zmid[0:1]
    p12     = zmid[1:2]

    N       = 100

    tspanf  = list(np.linspace((t0+tf)/2, tf, N+1))
    tspani  = list(np.linspace((t0+tf)/2, t0, N+1))

    xf, pf = f((t0+tf)/2., x12, p12, tspanf, e)
    xi, pi = f((t0+tf)/2., x12, p12, tspani, e)

    uf     = ufun(pf, e)
    ui     = ufun(pi, e)

    fig = plt.figure()

    ax  = fig.add_subplot(131);
    ax.set_xlabel('t'); ax.set_ylabel('$x$');
    ax.axhline(0, color='k', linewidth=0.5); ax.axvline( 0, color='k', linewidth=0.5); ax.axvline( tf, color='k', linewidth=0.5)
    ax.plot(np.concatenate((np.flip(tspani), tspanf)), np.concatenate((np.flip(xi), xf)))

    ax  = fig.add_subplot(132);
    ax.set_xlabel('t'); ax.set_ylabel('$p$');
    ax.axhline(0, color='k', linewidth=0.5); ax.axvline( 0, color='k', linewidth=0.5); ax.axvline( tf, color='k', linewidth=0.5)
    ax.plot(np.concatenate((np.flip(tspani), tspanf)), np.concatenate((np.flip(pi), pf)))

    ax  = fig.add_subplot(133);
    ax.set_xlabel('t'); ax.set_ylabel('$u$');
    ax.axhline(0, color='k', linewidth=0.5); ax.axvline( 0, color='k', linewidth=0.5); ax.axvline( tf, color='k', linewidth=0.5)
    ax.axhline(-1, color='k', linewidth=0.5)
    ax.axhline( 1, color='k', linewidth=0.5)
    ax.plot(np.concatenate((np.flip(tspani), tspanf)), np.concatenate((np.flip(ui), uf)))

Resolution of the regularized problem

[11]:

# Shooting for (tf, e) = (tf_init, e_init)
p0_guess = np.array([x_steady[0], 0.0])
sol_nle  = nt.nle.solve(shoot, p0_guess, df=shoot)


     Calls  |f(x)|                 |x|

         1  1.804903493012030e-04  1.000000000000000e-03
         2  8.904124663136239e-13  8.868188925481338e-04
         3  2.168404344971009e-19  8.868188931064897e-04

 Results of the nle solver method:

 xsol    =  [ 8.86818893e-04 -4.17259936e-19]
 f(xsol) =  [-2.16840434e-19  0.00000000e+00]
 nfev    =  3
 njev    =  1
 status  =  1
 success =  True

 Successfully completed: relative error between two consecutive iterates is at most TolX.


[12]:

# Plot solution for (tf, e) = (tf_init, e_init)
plotSolution(sol_nle.x, e_init, tf_init)
../../_images/turnpike_1d_bsb_turnpike_24_0.png 
[13]:

# Definition of the homotopic function and its first order derivative
# This function is used to solve S=0 for different values of e=epsilon and tf.
def dhomfun(x12, dx12, p12, dp12, x0, dx0, xf_target, dxf_target, e, de, tf, dtf):
    s      = np.zeros((2,))
    ds     = np.zeros((2,))
    (xf, dxf), _ = f(((t0+tf)/2., dtf/2.), (x12, dx12), (p12, dp12), (tf, dtf), (e, de))
    (xi, dxi), _ = f(((t0+tf)/2., dtf/2.), (x12, dx12), (p12, dp12), t0, (e, de))
    s[0]   = xf[0] - xf_target[0]
    s[1]   = xi[0] - x0[0]
    ds[0]  = dxf[0] - dxf_target[0]
    ds[1]  = dxi[0] - dx0[0]
    return s, ds

@tools.tensorize(dhomfun, tvars=(1, 2, 3, 4, 5, 6), full=True)
def homfun(x12, p12, x0, xf_target, e, tf):
    s      = np.zeros((2,))
    xf, _  = f((t0+tf)/2., x12, p12, tf, e)
    xi, _  = f((t0+tf)/2., x12, p12, t0, e)
    s[0]   = xf[0] - xf_target[0]
    s[1]   = xi[0] - x0[0]
    return s

#
def dfoo(zmid, dzmid, par, dpar):
    x12 = zmid[0:1]
    p12 = zmid[1:2]
    x0  = par[0:1]
    xf_target = par[1:2]
    e   = par[2]
    tf  = par[3]
    dx12 = dzmid[0:1]
    dp12 = dzmid[1:2]
    dx0  = dpar[0:1]
    dxf_target = dpar[1:2]
    de   = dpar[2]
    dtf  = dpar[3]
    return homfun((x12, dx12), (p12, dp12), (x0, dx0), (xf_target, dxf_target), (e, de), (tf, dtf))

@tools.tensorize(dfoo, tvars=(1, 2), full=True)
def foo(zmid, par):
    x12 = zmid[0:1]
    p12 = zmid[1:2]
    x0  = par[0:1]
    xf_target = par[1:2]
    e   = par[2]
    tf  = par[3]
    return homfun(x12, p12, x0, xf_target, e, tf)

[14]:

# Homotopy on x0, xf_target, e and tf
z0         = sol_nle.x
par0       = np.array([x_steady[0], x_steady[0], e_init, tf_init])
parf       = np.array([x0[0], xf_target[0], e_final, tf_final])
sol_path   = nt.path.solve(foo, z0, par0, parf, df=foo)


     Calls  |f(x,pars)|     |x|                Homotopic param    Arclength s     det A(s)        Dot product

         1  0.00000000e+00  8.86818893106e-04  0.00000000000e+00  0.00000000e+00  9.01434919e-01  0.00000000e+00
         2  7.43775334e-17  9.33573984338e-04  6.70794276397e-05  1.03075172e-04  9.03010300e-01  9.99998882e-01
         3  4.84869952e-19  9.64837418738e-04  1.08962612366e-04  1.67427190e-04  9.03036211e-01  9.99999997e-01
         4  6.85709662e-19  1.27184593066e-03  4.60774245338e-04  7.07783579e-04  9.03254085e-01  9.99999757e-01
         5  8.67361738e-19  5.18648967128e-03  3.99102264582e-03  6.11134747e-03  9.05476558e-01  9.99975544e-01
         6  3.46944695e-18  2.57896654845e-02  2.24054686124e-02  3.37688763e-02  9.18153664e-01  9.99337654e-01
         7  1.25092689e-17  5.74148021184e-02  5.27886841712e-02  7.76322374e-02  9.43189242e-01  9.98226605e-01
         8  3.46944695e-17  1.00646620601e-01  9.93762075234e-02  1.41259096e-01  9.92182643e-01  9.96012280e-01
         9  1.38777878e-17  1.52659546116e-01  1.65529599945e-01  2.25704323e-01  1.08596134e+00  9.92744888e-01
        10  5.55111512e-17  1.89293795722e-01  2.21306894398e-01  2.92870221e-01  1.18969049e+00  9.95556923e-01
        11  1.00074151e-16  2.18461572047e-01  2.73519445643e-01  3.53238945e-01  1.30989319e+00  9.96657657e-01
        12  2.77555756e-17  2.41814215487e-01  3.22496552056e-01  4.08135926e-01  1.44567721e+00  9.97441347e-01
        13  1.38777878e-16  2.60988921198e-01  3.69796792418e-01  4.59857244e-01  1.60115434e+00  9.97854101e-01
        14  2.23772605e-16  2.76736507048e-01  4.16150609217e-01  5.09507420e-01  1.78074245e+00  9.98054248e-01
        15  6.20633538e-17  2.89491443456e-01  4.62247997775e-01  5.58016980e-01  1.99148644e+00  9.98072736e-01
        16  1.57009246e-16  2.99389449225e-01  5.08573401511e-01  6.06030397e-01  2.24332160e+00  9.97929456e-01
        17  1.24126708e-16  3.06363655493e-01  5.55653637571e-01  6.54212747e-01  2.55217391e+00  9.97598998e-01
        18  1.57009246e-16  3.10106365685e-01  6.04127293989e-01  7.03356588e-01  2.94529167e+00  9.97009300e-01
        19  1.11022302e-16  3.09951152206e-01  6.54832535901e-01  7.54526068e-01  3.47368121e+00  9.95992054e-01
        20  3.37661151e-16  3.04981726150e-01  7.06376513332e-01  8.06702825e-01  4.20251291e+00  9.94722894e-01
        21  2.98936698e-16  2.95772672540e-01  7.52377517350e-01  8.53915367e-01  5.12876819e+00  9.94525749e-01
        22  2.77555756e-16  2.82969437321e-01  7.93494978840e-01  8.97222890e-01  6.35437842e+00  9.94280757e-01
        23  5.55111512e-17  2.66701502275e-01  8.30915482421e-01  9.38245761e-01  8.09398753e+00  9.93778186e-01
        24  2.77555756e-16  2.50326740271e-01  8.59883587997e-01  9.71699346e-01  1.02459775e+01  9.95178590e-01
        25  2.48253415e-16  2.35598909619e-01  8.81246696360e-01  9.97791662e-01  1.27104593e+01  9.96748182e-01
        26  4.57756680e-16  2.21947659717e-01  8.98179916350e-01  1.01967322e+00  1.56477541e+01  9.97556034e-01
        27  1.57009246e-16  2.09183431385e-01  9.12045930899e-01  1.03864441e+00  1.92079750e+01  9.98082496e-01
        28  3.55444798e-16  1.96024229003e-01  9.24698400764e-01  1.05703438e+00  2.40694199e+01  9.98148434e-01
        29  7.55033286e-16  1.83687053464e-01  9.35269594126e-01  1.07341555e+00  3.02427271e+01  9.98508073e-01
        30  4.44089210e-16  1.73295748220e-01  9.43324165245e-01  1.08668379e+00  3.72360385e+01  9.99015222e-01
        31  8.45520665e-16  1.62688475397e-01  9.50827985298e-01  1.09980930e+00  4.68939848e+01  9.99035494e-01
        32  4.47545209e-16  1.52921712963e-01  9.57151157024e-01  1.11157521e+00  5.91569052e+01  9.99226123e-01
        33  7.30135199e-16  1.43719347126e-01  9.62633789353e-01  1.12241931e+00  7.52309816e+01  9.99343505e-01
        34  1.55431223e-15  1.35189004722e-01  9.67330713384e-01  1.13228817e+00  9.61761743e+01  9.99455478e-01
        35  5.55111512e-16  1.27251257720e-01  9.71386435477e-01  1.14133127e+00  1.23821652e+02  9.99539627e-01
        36  5.11787527e-16  1.19869086209e-01  9.74898198049e-01  1.14963281e+00  1.60658599e+02  9.99606164e-01
        37  1.05471187e-15  1.12978365286e-01  9.77958464273e-01  1.15729606e+00  2.10475541e+02  9.99655592e-01
        38  5.11787527e-16  1.06517129710e-01  9.80642830523e-01  1.16441274e+00  2.79057555e+02  9.99691006e-01
        39  1.97123217e-15  1.00412994694e-01  9.83018182365e-01  1.17107910e+00  3.75768105e+02  9.99713386e-01
        40  4.96506831e-16  9.45766536101e-02  9.85146048437e-01  1.17740422e+00  5.16855528e+02  9.99722229e-01
        41  7.94798730e-16  8.88666370281e-02  9.87094124432e-01  1.18354839e+00  7.34427115e+02  9.99712048e-01
        42  2.50046775e-15  8.29078471906e-02  9.88987300284e-01  1.18991539e+00  1.11508711e+03  9.99651642e-01
        43  7.85046229e-16  7.84233956841e-02  9.90318361358e-01  1.19467750e+00  1.59188256e+03  9.99779068e-01
        44  8.84371215e-15  7.47257788893e-02  9.91355440424e-01  1.19858522e+00  2.20470266e+03  9.99834169e-01
        45  6.78054567e-15  7.11226362510e-02  9.92313465141e-01  1.20237684e+00  3.12841540e+03  9.99827145e-01
        46  3.29204997e-15  6.79371102338e-02  9.93117209516e-01  1.20571573e+00  4.39616231e+03  9.99852425e-01
        47  3.04203135e-14  6.50440251475e-02  9.93811970164e-01  1.20873739e+00  6.16179608e+03  9.99868378e-01
        48  2.78119732e-14  6.24294384108e-02  9.94411010546e-01  1.21145948e+00  8.58549779e+03  9.99884880e-01
        49  1.03417513e-13  6.00387019144e-02  9.94934794344e-01  1.21394137e+00  1.19196569e+04  9.99897895e-01
        50  1.90958393e-13  5.78386876521e-02  9.95396572127e-01  1.21621933e+00  1.65022217e+04  9.99909051e-01
        51  2.19935209e-13  5.57985872554e-02  9.95807484200e-01  1.21832671e+00  2.28139829e+04  9.99918351e-01
        52  5.06095361e-13  5.38955550458e-02  9.96175806190e-01  1.22028824e+00  3.15255051e+04  9.99926315e-01
        53  1.12854197e-12  5.21110162432e-02  9.96508085370e-01  1.22212397e+00  4.35828504e+04  9.99933184e-01
        54  2.06684670e-12  5.04303690917e-02  9.96809450937e-01  1.22384966e+00  6.03212566e+04  9.99939194e-01
        55  4.49212896e-12  4.88417894702e-02  9.97084023846e-01  1.22547807e+00  8.36345237e+04  9.99944500e-01
        56  1.05765397e-11  4.73356684135e-02  9.97335148746e-01  1.22701953e+00  1.16217803e+05  9.99949226e-01
        57  2.61988209e-11  4.59040738333e-02  9.97565583056e-01  1.22848261e+00  1.61922303e+05  9.99953463e-01
        58  6.33618158e-11  4.45403959797e-02  9.97777626609e-01  1.22987441e+00  2.26272331e+05  9.99957282e-01
        59  1.59496194e-10  4.32390667198e-02  9.97973219856e-01  1.23120093e+00  3.17223226e+05  9.99960740e-01
        60  4.06610579e-10  4.19953517352e-02  9.98154016544e-01  1.23246725e+00  4.46275664e+05  9.99963882e-01
        61  1.05659032e-09  4.08051886917e-02  9.98321439241e-01  1.23367775e+00  6.30120516e+05  9.99966745e-01
        62  2.76687706e-09  3.96650582821e-02  9.98476722454e-01  1.23483621e+00  8.93077332e+05  9.99969361e-01
        63  7.27046351e-09  3.85718770121e-02  9.98620946871e-01  1.23594594e+00  1.27072832e+06  9.99971755e-01
        64  1.90688532e-08  3.75229046521e-02  9.98755067057e-01  1.23700987e+00  1.81536901e+06  9.99973948e-01
        65  4.95636168e-08  3.65156610860e-02  9.98879934251e-01  1.23803065e+00  2.60425218e+06  9.99975958e-01
        66  1.27181858e-07  3.55478507865e-02  9.98996315217e-01  1.23901074e+00  3.75218983e+06  9.99977799e-01
        67  3.20787500e-07  3.46172954765e-02  9.99104907646e-01  1.23995244e+00  5.43107987e+06  9.99979483e-01
        68  8.17312655e-07  3.37218784375e-02  9.99206352235e-01  1.24085799e+00  7.90141477e+06  9.99981019e-01
        69  1.90668617e-06  3.28600169848e-02  9.99301185942e-01  1.24172907e+00  1.15550460e+07  9.99982437e-01
        70  4.56128184e-06  3.20285447670e-02  9.99390076535e-01  1.24256896e+00  1.70160360e+07  9.99983678e-01
        71  1.06303662e-05  3.12260254927e-02  9.99473462203e-01  1.24337917e+00  2.52540460e+07  9.99984822e-01
        72  2.45205324e-05  3.04501689382e-02  9.99551837429e-01  1.24416207e+00  3.78561440e+07  9.99985841e-01
        73  5.60902572e-05  2.96989933916e-02  9.99625630555e-01  1.24491970e+00  5.74849855e+07  9.99986756e-01
        74  1.28378074e-04  2.89704814963e-02  9.99695242826e-01  1.24565415e+00  8.88517207e+07  9.99987570e-01
        75  2.96757532e-04  2.82627470030e-02  9.99761034692e-01  1.24636736e+00  1.40903482e+08  9.99988296e-01
        76  7.02924299e-04  2.75740304576e-02  9.99823329714e-01  1.24706112e+00  2.32635077e+08  9.99988944e-01
        77  1.04964897e-03  2.69027503221e-02  9.99882413270e-01  1.24773707e+00  3.81903592e+08  9.99989522e-01
        78  4.25047489e-03  2.62475521286e-02  9.99938531828e-01  1.24839659e+00  8.15859970e+08  9.99990043e-01
        79  1.87028511e-02  2.56073834174e-02  9.99991890566e-01  1.24904077e+00  3.17920129e+09  9.99990519e-01
        80  9.99715086e-02  2.55085957196e-02  1.00000000000e+00  1.24914015e+00  3.15924546e+07  9.99999765e-01

 Results of the path solver method:

 xf            =  [ 0.02544342 -0.00182239]
 parsf         =  [1.e+00 5.e-01 1.e-03 2.e+00]
 |f(xf,parsf)| =  0.09997150862297144
 steps         =  80
 status        =  1
 success       =  True

 Homotopy successfully completed.


[15]:

# Plot final solution
if(np.abs(sol_path.parsf[2]-e_final)<1e-5):
    plotSolution(sol_path.xf, e_final, tf_final)
else:
    print('Homotopy not complete')
../../_images/turnpike_1d_bsb_turnpike_27_0.png

III) Resolution of the optimal control problem by multiple shooting

We come back to the original optimal control problem:

\[\begin{split}\left\{ \begin{array}{l} \displaystyle J(u) := \displaystyle \int_0^{t_f} x^2(t) \, \mathrm{d}t \longrightarrow \min \\[1.0em] \dot{x}(t) = f(x(t), u(t)) := \displaystyle u(t), \quad |u(t)| \le 1, \quad t \in [0, t_f] \text{ a.e.}, \\[1.0em] x(0) = 1, \quad x(t_f) = 1/2. \end{array} \right.\end{split}\]

We have determined that the optimal control follows the strategy:

\[\begin{split}u(t) = \left\{ \begin{array}{lll} -1 & \text{if} & t \in [0, t_1], \\[0.5em] \phantom{-}0 & \text{if} & t \in (t_1, t_2], \\[0.5em] +1 & \text{if} & t \in (t_2, t_f], \end{array} \right.\end{split}\]

with \(0 < t_1 < t_2 < t_f=2\).

Goal

The goal is to find the values of the switching times \(t_1\) and \(t_2\) together with the initial covector \(p_0\) (see Remark 2).

Maximized Hamiltonian and its derivatives

We define first the three control laws \(u \equiv \{-1, 0, 1\}\).

[16]:

# Controls in feedback form
@tools.vectorize(vvars=(1, 2, 3))
def uplus(t, x, p):
    u = +1.0
    return u

@tools.vectorize(vvars=(1, 2, 3))
def uminus(t, x, p):
    u = -1.0
    return u

@tools.vectorize(vvars=(1, 2, 3))
def using(t, x, p):
    u = 0.0
    return u

The pseudo-Hamiltonian is

\[H(x,p,u) = pu - x^2.\]
[17]:

# Definition of the Hamiltonian and its derivatives for u = 1
#
def dhfunplus(t, x, dx, p, dp):
    # dh = dh_x dx + dh_p dp
    u  = uplus(t, x, p)
    hd = u*dp - 2.0*x*dx
    return hd

def d2hfunplus(t, x, dx, d2x, p, dp, d2p):
    # d2h = dh_xx dx d2x + dh_xp dp d2x + dh_px dx d2p + dh_pp dp d2p
    hdd    = -2.0 * d2x * dx
    return hdd

@tools.tensorize(dhfunplus, d2hfunplus, tvars=(2, 3))
def hfunplus(t, x, p):
    u = uplus(t, x, p)
    h = p*u - x**2
    return h

hplus = ocp.Hamiltonian(hfunplus)
fplus = ocp.Flow(hplus)

We give the Hamiltonians for \(u=-1\) and \(u=0\) with their derivatives.

[18]:

# Definition of the Hamiltonian and its derivatives for u = -1
def dhfunminus(t, x, dx, p, dp):
    # dh = dh_x dx + dh_p dp
    u  = uminus(t, x, p)
    hd = u*dp - 2.0*x*dx
    return hd

def d2hfunminus(t, x, dx, d2x, p, dp, d2p):
    # d2h = dh_xx dx d2x + dh_xp dp d2x + dh_px dx d2p + dh_pp dp d2p
    hdd    = -2.0 * d2x * dx
    return hdd

@tools.tensorize(dhfunminus, d2hfunminus, tvars=(2, 3))
def hfunminus(t, x, p):
    u = uminus(t, x, p)
    h = p*u - x**2
    return h

hminus = ocp.Hamiltonian(hfunminus)
fminus = ocp.Flow(hminus)

Shooting function

The multiple shooting function is given by

\[\begin{split}S(t_1, t_2) := \begin{pmatrix} x(t_0, t_1, 0, 0, u_-) - 1 \\ x(t_f, t_2, 0, 0, u_+) - 1/2 \end{pmatrix},\end{split}\]

where z(t, s, a, b, u) is the solution at time \(t\) of the Hamiltonian system associated to the control u starting at time \(s\) at the initial condition \(z(s) = (a,b)\).

We have introduced the notation \(u_-\) for \(u\equiv -1\) and \(u_+\) for \(u\equiv +1\).

Remark: We know that \((x(t_1), p(t_1))=(x(t_2), p(t_2))=(0,0)\).

[19]:

# Multiple shooting function
#

tf = tf_final # we set the final time to the value tf_final

def shoot_multiple(y):
    t1 = y[0]
    t2 = y[1]

    zero = np.zeros([1,])
    xi, pi = fminus(t1, zero, zero, t0)  # on [t1, t0]
    xf, pf =  fplus(t2, zero, zero, tf)  # on [t2, tf]

    s = np.zeros([2])
    s[0] = xi[0] - x0[0]
    s[1] = xf[0] - xf_target[0]  # x(tf, x0, p0) - xf_target
    return s

Resolution of the shooting function

[20]:

# Initial guess for the Newton solver
t1_guess = 0.9
t2_guess = 1.6

[21]:

# Resolution of the shooting function
y_guess  = np.array([t1_guess, t2_guess])
sol_nle_mul  = nt.nle.solve(shoot_multiple, y_guess)


     Calls  |f(x)|                 |x|

         1  1.414213562373096e-01  1.835755975068582e+00
         2  5.058608532060305e-10  1.802775638152897e+00
         3  1.110223024625157e-16  1.802775637731995e+00

 Results of the nle solver method:

 xsol    =  [1.  1.5]
 f(xsol) =  [0.00000000e+00 1.11022302e-16]
 nfev    =  3
 njev    =  1
 status  =  1
 success =  True

 Successfully completed: relative error between two consecutive iterates is at most TolX.


[22]:

# function to plot solution

def plotSolutionBSB(t1, t2, tf):

    N      = 20

    tspan1  = list(np.linspace(t1, t0, N+1))
    tspan2  = (t1, t2)
    tspanf  = list(np.linspace(t2, tf, N+1))

    zero    = np.zeros([1,])

    x1, p1 = fminus(t1, zero, zero, tspan1)  # on [ 0, t1]
    x2     = np.array([zero, zero])
    p2     = np.array([zero, zero])
    xf, pf = fplus(t2, zero, zero, tspanf)  # on [t2, tf]

    u1     = uminus(tspan1, x1, p1)
    u2     = np.array([0., 0.])
    uf     =  uplus(tspanf, xf, pf)

    fig = plt.figure()

    ax  = fig.add_subplot(131);
    ax.set_xlabel('t'); ax.set_ylabel('$x$'); ax.axhline(0, color='k', linewidth=0.5);
    ax.axvline( 0, color='k', linewidth=0.5); ax.axvline( tf, color='k', linewidth=0.5)
    ax.plot(np.concatenate((np.flip(tspan1), tspan2, tspanf)), np.concatenate((np.flip(x1), x2, xf)))

    ax  = fig.add_subplot(132);
    ax.set_xlabel('t'); ax.set_ylabel('$p$'); ax.axhline(0, color='k', linewidth=0.5);
    ax.axvline( 0, color='k', linewidth=0.5); ax.axvline( tf, color='k', linewidth=0.5)
    ax.plot(np.concatenate((np.flip(tspan1), tspan2, tspanf)), np.concatenate((np.flip(p1), p2, pf)))

    ax  = fig.add_subplot(133);
    ax.set_xlabel('t'); ax.set_ylabel('$u$'); ax.axhline(0, color='k', linewidth=0.5);
    ax.axvline( 0, color='k', linewidth=0.5)
    ax.axvline( tf, color='k', linewidth=0.5)
    ax.axhline(-1, color='k', linewidth=0.5)
    ax.axhline( 1, color='k', linewidth=0.5)
    ax.plot(np.concatenate((np.flip(tspan1), tspan2, tspanf)), np.concatenate((np.flip(u1), u2, uf)))

[23]:

# plot solution
t1 = sol_nle_mul.x[0]
t2 = sol_nle_mul.x[1]
plotSolutionBSB(t1, t2, tf)
../../_images/turnpike_1d_bsb_turnpike_44_0.png

thumbnail

previous

SIR model

next

Turnpike property and the singular flow