BVPs Codes for Solving Optimal Control Problems

: Optimal control problems arise in many applications and need suitable numerical methods to obtain a solution. The indirect methods are an interesting class of methods based on the Pontryagin’s minimum principle that generates Hamiltonian Boundary Value Problems (BVPs). In this paper, we review some general-purpose codes for the solution of BVPs and we show their efﬁciency in solving some challenging optimal control problems. The continuation stiff size, convergence of nonlinear discretization schemes. as the of the codes twpbvpc_m and twpbvpc_l are the same of by the codes twpbvp_m and twpbvp_l , conﬁrming the non-necessity of these codes to use a mesh selection strategy based on conditioning for this non-stiff problem.


Introduction
Many optimal control problems arise from an interest in observing the dynamic behavior of a state variable described by a dynamic equation, namely by a differential equation, in several areas of applications such as biology, chemistry, economy, physics, and engineering. For example, we can consider the development of a specific species of animals in an ecological preserve, the dynamical behavior of a chemical process, the evolution of the selling trend of a company, or the simulation of high-performance racing vehicles. The dynamical behavior of this kind of problems is influenced by the choice of control variables, as it might be incorporating the presence of predators in the ecological preserve, moreover, both state and control variables must fulfil constraints, and minimize or maximize an objective function.
Numerical methods solving optimal control problems were considered starting from the 1950s, when Bellman introduced the dynamic programming [1], that requires solving a partial differential equation, called the Hamiltonian-Jacobi-Bellman equation. Through time the numerical approaches can be mainly divided into two classes: direct methods and indirect methods [2,3]. Perhaps the first class of direct methods is the most widely applied, it transforms the problem into a nonlinear optimization problem or nonlinear programming problem, essentially this class is focused on the use of optimization techniques. The second class of the indirect methods transforms the original optimal control problem into a twopoint boundary value problem, highlighting particular attention to numerical methods solving differential equation systems. The last strategy is often considered disadvantageous for figuring out challenging optimal control problem.
As against this last opinion, this work aims to review many of the available generalpurpose codes solving boundary value problems, able to figure out optimal control problems arising from adopting an indirect approach. The review is also devoted to some numerical strategies that are useful and sometime necessary to numerically solve the problem, such as continuation techniques associated with suitable penalty functions.
The most used solver for indirect methods has been the shooting method, based on guessing the value of the unknown boundary condition at one end of the interval, so that an initial value problem is solved to obtain the solution at the other end of the interval that is already known. Although the shooting method is simple to apply, it is not particularly advantageous to use when the boundary value problem is ill conditioned or stiff, and are illustrated and classified through different programming environments, in particular Fortran codes are allocated in Section 3.1, Matlab codes in Section 3.2 and R codes in Section 3.3. Finally, in Sections 4-8 interesting optimal control problems are solved using indirect methods and the BVPs related. In Section 9 we give some conclusions, highlighting the potentiality of the BVP codes considered.

Optimal Control Problems: Indirect Methods
Given a non-empty compact time interval [t 0 , t f ] ⊂ R, with t 0 < t f , an optimal control problem is defined as where ϕ and L are sufficiently smooth functions involved in the minimization of the objective function, x(t) ∈ R n is the state variable of the dynamical system, u(t) ∈ U ⊂ R m is the control variable and U the set of admissible controls, f is a regular function and b(x(t 0 ), x(t f )) = 0 are the general boundary conditions. Furthermore, Problem (1) might be subject to a path constraint that can be expressed by a mixed control-state constraint c(t, x, u) ≤ 0 or a pure state constraint s(t, x) ≤ 0 . There exist two main approaches solving optimal control problems (1), direct methods and indirect methods [2,29]. Direct methods suitably discretize an infinite-dimensional optimal control problem, giving back a finite-dimensional optimization problem that can be solved using appropriate nonlinear programming methods, such as sequential quadratic programming. This approach results robust and efficient if applied to several problems, besides not requiring a strong knowledge in optimal control theory, it becomes highly advantageous to use.
On the other hand, indirect methods are instead related to the Pontryagin's minimum principle [29], a necessary condition for optimality that transforms the original Problem (1) into a two-point boundary value problem for state and adjoint Lagrange multiplier functions, defined as x = f (t, x, u), where H(t, x, u, λ) = L(t, x, u) + λ λ λ · f (t, x, u) is the Hamiltonian function and the optimal control u * (t) is obtained by a local optimization of the Hamiltonian, namely u * (t) = arg min u∈U H(t, x, u, λ). Pro this approach there is the possibility to compute an accurate numerical solution; however, against we find some drawbacks, such as the necessity to have a good initial guess for the solution of the generated nonlinear boundary value problem. Now, to overcome this matter we focus on the application of some well-known two-point boundary value codes that are considered extremely efficient and robust to solve the BVP (2).

Codes for BVPs
Boundary value problems arise in many fields of application, so in the last 40 years a great effort has been done to develop efficient methods solving this kind of problems.
Among them many are methods applied to two-point boundary value problems, i.e., to systems of first-order ordinary differential equations with boundary conditions, others can be applied directly to second or high-order boundary value problems without any transformation of the original problem. Moreover, these codes are available in different programming environment, so in the following we will give information about their characteristics.

Fortran Codes
The code colsys was written by U. Ascher, R. Matteij and R. Russell [5] and it is based on method of spline collocation at Gaussian points and solves mixed-order systems of multipoint BVPs, high-order equations, problems with non-separated boundary conditions and problems with singularity. The code computes the solution on a sequence of meshes that are refined using the equidistribution of error to satisfy the required input tolerance. The error estimate is obtained roughly at each step halving the mesh. The components of the collocation solution are expressed by B-spline basis, which are evaluated by the de Boor's algorithms. Indeed, the damped Newton's method of quasilinearization is used for solving the nonlinear problems.
The code colnew [6,30] is the descendant of colsys and, contrary to this last, it uses a Runge-Kutta monomial representation for the piecewise polynomial solution, instead of B-spline basis. This change returns a code faster than the native version colsys.
The codes twpbvp,twpbvpl and acdc were written by J. R. Cash and his collaborators. The code twpbvp [7], differently from colsys, uses mono-implicit Runge-Kutta formulae and a deferred correction method for solving two-point boundary value problems. The mono-implicit Runge-Kutta formulae are implemented applying the deferred correction procedure, which allows discovery of the solution of a high-order method using only low order schemes. The code guarantees to construct a mesh refinement that is very suitable for singular perturbation problems.
The code twpbvpl, differently from twpbvp, is based on three Lobatto Runge-Kutta formulae of order 4, 6, 8, which are implemented using a suitable deferred correction scheme, solved with a damped Newton iteration scheme. The code is devoted in solving efficiently nonlinear stiff two-point boundary value problems.
The code acdc [9] has been developed from twpbvpl including an automatic continuation strategy, implemented to suitably solve linear and nonlinear singular perturbation problems characterized from a small parameter . The parameter often brings about stiffness in the problem, so that for a nonlinear problem a good initial solution is required to reach the convergence of the Newton method. The continuation strategy arises to overcome these matters, specifically it consists of selecting an initial perturbation parameter 0 , chosen to compute a solution of a problem not particularly stiff, usually for 0 ≈ 1, and satisfying a certain exit tolerance tol. The idea is to obtain an initial rough profile of the solution of the problem for a desired perturbation parameter . Then, chosen an integer N the interval [ 0 , ] is discretized in N subintervals, so that 0 > 1 > · · · > · · · > N −1 > N . Now, N + 1 boundary value problems satisfying an exit tolerance tol are iteratively computed, so that the solution of the problem obtained at iteration i = 0, . . . , N − 1, for i on a mesh π i , is the initial solution of the next problem with perturbation parameter i+1 . A crucial point of this strategy is the selection of the initial parameter 0 and the value of discretizazion N , both depend on the problem. In codes such as acdc 0 is set equal to 0.5 by default; however the suggestion is to consider 0 as a value not extremely small allowing the obtaining of an accurate solution of the problem for that value of perturbation; The code acdc chooses the sequences of parameters and the total number of continuation steps automatically. It is however possible to implement a continuation strategy for the other codes, in this case for N it would be convenient to start with a small integer and then double or increment it, if the procedure does not converge.
The code colmod [9] is a modified version of the code colsys using the same continuation strategy adopted in acdc.
The codes twpbvpc, twpbvplc and acdcc [14] are the modified version of the codes twpbvp, twpbvpl and acdc that implement a mesh selection strategy based on the estimation of the local error and of two conditioning parameters [31]. This hybrid mesh strategy has first been used in the Matlab code TOM, described in the next section.
The code mirkdc written by W. Enright and P. Muir [11] uses MIRK method and controls the defect, also BVP_M-2 written by J.J. Boisvert, P. Muir and R. Spiteri [12] is based on MIRK methods, but this last controls both the defect and/or the global error, giving, moreover, information about the conditioning constant.
Detailed information about all the numerical schemes and techniques related to the Fortran codes in this subsection can be found in [32] where a review of global methods for solving BVPs is presented.

Matlab Codes
The BVP codes available officially in the Matlab environment are bvp4c [33] and bvp5c [34]. The code bvp4c [16] is based on a collocation method with a C 1 piecewise cubic polynomial, or equivalently on an implicit Runge-Kutta formula with a continuous extension, namely the collocation method is equivalent to a three-stages Lobatto IIIa implicit Runge-Kutta formula. This code implements a method of order four and solves a large class of BVP, such as equations with non-separated boundary conditions, singular problems, Sturm-Liouville problems. An advantage of this code is being able to compute numerical partial derivatives and use a vectorized finite difference Jacobian. Differently from the other codes the error estimation and the mesh selection are based on the residual estimation. We recall that if S(x) approximates the solution y(x), then the residual control in the differential equation The code bvp5c is based on the four-stages Lobatto IIIa formula, giving a method of order five. Contrarily to bvp4c, bvp5c controls the residual and the true approximate error. It is clear that if the BVP is well-conditioned a small residual implies a small true error, but this is not satisfied if the BVP is ill-conditioned, hence the strategy to control the residual and the true error is more efficient than the one applied in bvp4c.
The next two codes TOM and HOFiD_bvp belong to the class of Boundary Value Methods [35], especially suitable for solving BVPs.
The code TOM [19], based on the TOP Order Methods and the BS method of order four, six, eight and ten distinguishes for the use of conditioning in the mesh selection strategy. In [36] the authors analyzed how the conditioning and the stiffness of a problem depend on the estimation of the following conditioning parameters: κ conditioning constant with respect to all type of perturbation, computed using the maximum norm; κ 1 conditioning constant with respect to a perturbation of the boundary conditions, computed using the maximum norm; κ 2 conditioning constant with respect to a perturbation of the differential problem, computed using the maximum norm; γ 1 conditioning constant with respect to a perturbation of the boundary conditions, computed using the one norm; σ the stiffness ratio.
Specifically, the problem is: well-conditioned if κ, κ 1 , γ 1 and σ are of moderate size; stiff if σ 1; ill-conditioned if κ 1 and γ 1; ill posed if κ 2 > κ 1 . A complete description of the parameters and the algorithms used to compute their approximation is presented in [37]. The hybrid mesh selection algorithm controls the approximation of conditioning parameters and chooses the mesh points to have an estimation of those discrete quantities close to the continuous ones. Meanwhile, the code controls that the error of the solution computed is less than a prescribed tolerance. The error approximation is computed using a deferred correction technique with a higher order method, moreover a quasi-linearization technique is implemented to solve nonlinear problem. The release of May 2021, which has been used for the numerical tests in this paper, has the possibility to choose two different mesh selections, one suitable for regular problem and the other one for stiff or singular perturbation problems.
The code HOFiD_bvp [20] is based on high-order finite difference schemes (HOFiD) of order four, six, eight and ten, and an upwind method. Each derivative in the high-order boundary value problem is approximated directly by these schemes, hence it is not required any transformation of the problem in a system of first-order differential equations. The error estimation is computed applying the deferred correction technique to two consecutive order methods. The mesh selection is based on the error equidistribution. For nonlinear problems, the code uses a continuation strategy, as explained previously, and also combines an order variation strategy, this means that a solution of the problem obtained with a lower order and tolerance can be considered to be initial solution to run the code with higher order and tolerance. The strategy adopted returns a code suitable to solve high-order boundary value problems that can be singularly perturbed, singular, with discontinuous terms and multipoint. Other versions of the code solve singular second order initial value problems [38], Sturm-Liouville problems [39] and multi-parameters spectral problems [40].
An interesting code for solving high-order BVPs is the bvpSuite2.0 package, based on collocation methods. The collocation points could be chosen by the users among Gauss, Lobatto, uniform or user defined points. The code solves implicit BVPs, eigenvalue problems, differential algebraic problems of index 1 and it is particularly suited for singular problems. BvpSuite2.0 [21] is the evolution of two previous versions of the code with improved usability. The mesh selection strategy used is described in [41].
Finally, we consider the Matlab code bvptwp [18] based on an efficient translation of the Fortran codes twpbvp, twpbvpl and acdc in the Matlab environment, which are named twpbvp_m, twpbvp_l, acdc. Moreover, the Matlab package also contains the translation of the Fortran version of the same codes that use a hybrid mesh selection based on conditioning, similar to the one used in the code TOM, called twpbvpc_m, twpbvpc_l, acdcc. The code bvptwp is available on the calgo website and on the web-page called Test Set for BVP Solvers [15]. The version used in this paper is the release of May 2021.

R Codes
In recent years, the use of the open-source software R is upward among the problemsolving environments (PSEs) , and although it is mainly used as a software for statistics and visualization, several powerful methods solving differential equations have been developed. In this regard we highlight the package bvpSolve [23], which, using an interface, implements all the Fortran codes introduced in Section 3.1.

Experiments
Since our aim is to show the suitability and the efficiency of the BVP solvers in computing the solution of the Hamiltonian boundary value problems deriving from the application of the indirect method to optimal control problems, in the following sections we carry out some interesting numerical tests. We run experiments using the Matlab codes bvp4c, bvp5c, and bvptwp. For the last solver we consider all the codes available, i.e., twpbvp_m, twpbvp_l, twpbvpc_m, twpbvpc_l, acdc, acdcc. We also add the results obtained with the new release of the code TOM (May 2021). This code allows the choice of a boundary value method of specific order and a mesh variation strategy. For all the examples we choose the BS method of order 4 and we denote by tom the code run using a mesh variation for regular problems and by tomc the one implementing a mesh variation suited for stiff problems. For R-users all the examples could be solved applying all the codes included in the bvpSolve package. Since bvpSolve runs the Fortran codes by an interface, the obtained results are similar to those computed by means of the original Fortran codes. We also observe that some of the codes considered here for the numerical tests are also present in the R package bvpSolve rel. 1.4.2. The R version of these codes on the same examples show comparable results.
In our tests we use an initial mesh with 16 equidistant points and an initial solution with zero elements, except in some examples where specified. Moreover, the maximal mesh allowed has been set to 10 4 and the function evaluations have been vectorized. In the tables we report the number of points in the final mesh f M (in reading this value we recall that the code TOM does not use any auxiliary steps but all the others codes needs also several intermediate steps depending on the order of the methods used), the total number of vectorized function evaluation NVF and the mixed relative error on some significant components of the solution defined for a generic component x by the following formula where x i is the numerical approximation of x(t i ). If the exact solution of the test problem is not available, the error is computed by running the code twpbvpc_l using a doubled mesh and a halved input tolerance. For all the codes we give in input equal absolute and relative tolerances. If the codes twpbvp_m/twpbvpc_m, twpbvp_l/twpbvpc_l, acdc/acdcc give the same results we report only one result in the tables. If a code cannot solve the problem, we put * in the tables.

Hypersensitive Optimal Control Problems
The first class of examples we consider is the class of hypersensitive optimal control problems. Problems in this class are stiff, and need a suitable mesh variation strategy when solved using both direct and indirect methods. Usually, they are considered extremely difficult to be solved by indirect methods, because the solution is sensitive to changes in the initial conditions. In [42] the authors describe a dichotomic basis method which is inspired to the computation of the solution of singular perturbation problems for stiff initial value problems. In the following examples we show that general-purpose finite differences codes can solve very efficiently this class of problems. The codes can be applied for the numerical solution of completely hypersensitive problems whose solution has fast rates in all directions and partially hypersensitive problems, with the fast rate in only one direction.

Nonlinear Mass Spring System with Quadratic Cost
As first example we consider a hypersensitive nonlinear mass spring system [43], where the mass position x is defined such that the spring is unstretched when The control is exerted on the mass by an external force denoted by F(t), hence the control input is u(t) = F(t). The equation of motion of the mass is mx = F s (x) + F(t). We assume that k 1 = 1, k 2 = 1 and m = 1.
The optimal control problem needs to determine the control u on the fixed time

The associated Hamiltonian is
and the optimal control, obtained by computing ∂H ∂u = 0, is given by u * = −µ. Therefore, applying the indirect method the optimal control problem is equivalent to solve the following BVP In Figure 1 we show the solution for T = 20 and T = 40. In Table 1 we present some results obtained increasing the value of T from 20 to T = 2 × 10 6 . First, we choose an initial mesh of 16 equidistant points and try to run all the codes, except the codes acdc and acdcc, since for this formulation of the problem there is not a parameter to be used for continuation. If on one hand, for T = 20 all the methods converge to the solution, and for T = 2 × 10 4 only the codes bvp4c and bvp5 fail, on the other hand for T = 2 × 10 6 no one goes to convergence except the codes tom and tomc (see Table 2). Essentially, there are some troubles with a singular Jacobian for bvp4c and bvp5c, or a drawback with the maximum number of mesh points allowed with the other codes. In the last case we could increase the maximum value of mesh points; however, we will try to differently overcome this matter and to debunk the idea that the indirect methods are not as competitive as direct ones.   First, we point out that the results presented in Table 2 clearly show that the mesh selection based on conditioning allows the solution of the problem using a reduced number of mesh points and vectorial function evaluations. To gain the convergence for the other codes, it can be sufficient, in some cases, to increase the number of points in the initial mesh. To this aim, in Table 3 we show the numerical results obtained for bvp5c using 501 or 1001 initial equidistant points and T = 2 × 10 4 . This strategy is advantageous for bvp5c, not yet for bvp4c, that needs an initial mesh of 2501 mesh points to reach the convergence. However, we observe that bvp5c is not able to reach convergence if we use an initial mesh of 2501 mesh points. For the other classes of methods increasing the number of mesh points is not advantageous in terms of computational cost and time execution. To improve the performance of all considered codes, the BVP (3) is reformulated using a variable transformation. Let τ = t/T with τ ∈ [0, 1], we solve the following BVP Now, we set the perturbation parameter = 1/T, so that we can run for parameters less than 1 the codes acdc and acdcc that use an automatic continuation strategy. For all the other codes, we can adopt a continuation strategy starting with an initial value of that we discretize the interval with N = 3 and N = 5 respectively for T = 2 × 10 4 and T = 2 × 10 6 . In Table 4 we show the results obtained applying this successful continuation strategy. All the methods converge for all the values of T using a low computational cost. In this case, the codes based on automatic continuation strategy are very efficient, using acdc/acdcc the users do not need to decide how to change the continuation parameters, even if in some cases the automatic continuation could fail to reach the final desired value.
More information about this problem could be obtained by analyzing the conditioning parameters given in output by the codes twpbvpc_m and tomc, reported in Table 5. As we can see the stiffness parameter σ grows with the width of the interval, and depends on this last, moreover κ 2 > κ 1 shows that the problem could be ill posed, and γ 1 tending to zeros shows the presence of different time scales. The transformation of time interval in [0, 1] does not change the stiffness of the problem, but the problem is well posed (see Table 6).

Completely Hypersensitive Control Problem
This example is a hypersensitive optimal control problem implemented in ICLOCS2, defined as a problem "extremely difficult" to solve using an indirect method [42,44] and given by , the first-order necessary conditions for optimality leads to the following boundary value problem (BVP) where the optimal control is u * = − λ 2 . We choose T = 10 4 , T = 10 6 and an initial mesh of 11 equidistant points, the solution is plotted in Figure 2. Numerical results shown in Table 7 point out good performance of all the codes except bvp4c and bvp5c, which are not suitable for stiff problems, indeed we underline as they converge to the solution respectively up to T = 38 and T = 29.
In Table 8 the approximations of the conditioning constants show the dependence of the stiffness on the width of the interval. Moreover, the numerical results underline the necessity of adopting a good mesh selection strategy for computing the solution.      For the purpose of improving the performance and overcoming some drawbacks, we propose, as already done for the test problem in Section 4.1, to use the transformation of the variable τ = t/T, such that the BVP (6) can be reformulated for τ ∈ [0, 1] as x(0) = 1, x(1) = 1.5.
The advantage of this formulation is that considering = 1/T as a perturbation parameter, we can apply the continuation strategy on that parameter. In Table 9 we report the results using as starting value 0 = 1/10 and changing the continuation parameters in the interval [ 0 , ] among the value of the set 10 −2 , 10 −3 , 10 −4 . We remember that acdc and acdcc, using an automatic continuation strategy, needs only to insert the desired value of and uses as 0 the default value 0.5. The numerical tests and the conditioning parameters in Tables 8 and 10 clearly show that for this class of problems, if we cannot use a continuation of parameters, the codes able to give a solution are the ones suited for stiff problems that work still better if also the mesh selection is appropriate for this class of problems. Table 9. Hypersensitive problem using the variable τ = t/T, initial mesh with 11 equidistant points and continuation strategy on T: final mesh (fM), total number of vectorized function evaluation (NVF) and mixed errors for x,v,u.

Bang-Bang Optimal Control Problem
The bang-bang optimal control problem [45] is among the more challenging ones. It arises from a model in which a point unit mass m subjects to a limited force in onedimensional space, i.e., mx (t) = u(t) and u(t) ≤ 1. The main feature of optimal control problem of moving the mass from x = 0 to the maximum distance x in one second can be formulated as follows The associated Hamiltonian function is defined as H(x, v, λ, µ, u) = −v + λv + µu and the optimal control is given by Now, by applying the indirect method the solution of the optimal control Problem (8) is equivalent to solve the following BVP problem We observe that the optimal control is defined as and the exact solution is given by The discontinuity of the switching function is overcome by a smoothing technique that can be executed by different strategies. We choose two of them in particular. The first strategy, given a small parameter , consists of using the approximation The exact bang-bang solution is better approximated when becomes smaller; however, this for value around smaller than 10 −4 can give ill-conditioning problems. Table 11 contains all the results obtained using the Matlab codes, the solution is plotted in Figure 3. Only bvp5c fails, and for getting the solution is necessary to use the continuation strategy. To this regard we consider as initial perturbation parameter 0 = 1 and then we change it choosing N = 10 logarithmically equispaced points between 1 and the value required . When tol = 10 −4 bvp5c converges using 19 points for both equal to 10 −3 and 10 −6 , instead when tol = 10 −4 bvp5c gets the solution with 36 and 28 points respectively for = 10 −3 and = 10 −6 . Table 11. Bang-Bang optimal control Problem (9) For the second smoothing technique we can add a barrier or a penalty function. In this regard, we consider a piecewise quadratic penalty function defined as in [45] where the parameter σ gives the distance from the border where the penalty changes fast. Consequently, the Problem (8) is reformulated without inequality constraint as follows The optimal control u, obtained as a solution of the equation In Table 12 we show the numerical results obtained for σ = 10 −4 and = 10 −4 , 10 −6 , starting with an initial mesh of 16 equidistant points and a null initial solution. It is clear that all the codes have a good performance, we do not report the results for bvp4c and bvp5c because they fail. To overcome this drawback in Table 13 we consider the continuation strategy, this means that the codes bvp4c and bvp5c are run for different values of starting from 0 = 10 up to the desired value . In particular, we choose N = 10 values logarithmically equispaced. Table 12. Bang-Bang optimal control problem-solving (8) Table 13. Bang-Bang optimal control problem-solving (8)  We also report the results of acdc and acdcc that use an automatic continuation strategy. The results point out the suitability and efficiency of the strategy in solving this kind of problems, also for bvp4c and bvp5c when the nonlinear solution is approximated using a continuation strategy. The conditioning parameters reported in Tables 14 and 15 show that the problem is not stiff since σ is of moderate size, indeed the main difficulty is caused by the convergence of the nonlinear discretization schemes. In this regard we highlight as the results of the codes twpbvpc_m and twpbvpc_l are the same of those gained by the codes twpbvp_m and twpbvp_l, confirming the non-necessity of these codes to use a mesh selection strategy based on conditioning for this non-stiff problem. Table 14. Bang-Bang optimal control problem: conditioning parameters computed using tol = 10 −6 .

Longitudinal Dynamics of a Vehicle
We consider an example of nonlinear optimal control problem derived from a model of the longitudinal dynamics of a vehicle with the aerodynamic down-force [2]. In particular, a vehicle, supposed to be a point mass, is moved in a fixed time T from an initial zero velocity to a final zero velocity The Hamiltonian function associated with this problem is and the optimal control is given by Now, applying the indirect method the global optimal control problem is reduced to the boundary value problem We observe that the optimal control problem has a theoretical solution given by that can be approximated using a barrier function defined as π .
Let g + = g + k 0 and g − = g − k 0 , if t s = 1 k 1 ln g − + g + e k 1 T 2g is the switching time, then the solution for the optimal control is defined as Moreover, the exact solution for the space and the velocity is expressed by while the multipliers assume the form In Table 16 are shown all the numerical results obtained using all the Matlab codes considered starting with an initial mesh of 11 equispaced points and an initial approximation with null elements, the solution is plotted in Figure 4. For this problem only the codes of the bvptwp package are able to give a solution for = 10 −6 , so for the other codes we have used a continuation strategy with a starting value 0 = 10 −3 and N = 10 logarithmic equispaced intermediate points. In Table 16 all the results obtained are shown in order that the symbol c in bracket labels those computed using the continuation strategy. Moreover, the results emphasize that not always the automatic continuation is advantageous and cheaper from a computational cost of view, since it is evident that the total number of vectorial functions evaluation is much greater for acdc than for twpbvp_m and twpbvp_l. Remember that they use the same numerical scheme. The conditioning parameters in Table 17 are all moderate size, hence the problem is not stiff. Table 16. Longitudinal dynamics of a vehicle T = 10, g = 9.81, k 0 = 0.02 g, k 1 = 10 −5 g, k 2 = 0, k 3 = 0: final mesh (fM), total number of vectorized function evaluation (NVF) and mixed errors for x, v, u.  Table 17. Longitudinal dynamics of a vehicle T = 10, g = 9.81, k 0 = 0.02 g, k 1 = 10 −5 g, k 2 = 0, k 3 = 0: conditioning parameters computed using tol = 10 −6 .

Gottard Rocket
Now, we consider an example of optimal control problem with a singular arc [46]. A rocket of mass m lifts off vertically at time t = 0 with (normalized) altitude h(0) = 1 and velocity v(0) = 0. Known the initial mass, the fuel mass and the drag characteristics of the rocket, the aim is to choose the thrust u(t) and the final time T to maximize the altitude h(T) at the final time T. The optimal control problem is given by Given the constants D c and h c , the aerodynamic drag is defined by . Moreover, if g 0 is the gravitational force at the earth's surface, then the gravitational force is given by The equation is scaled choosing the model parameters m(0), h(0) and g 0 , which allows management of dimension-free equations. As in [46], we consider where g 0 = 1, h c = 500, m c = 0.6 and v c = 620.
Since the problem (13) has a free final time, we fix the time interval using the variable transformation t(τ) := τT, with τ ∈ [0, 1]. A new state variable T satisfying the differential constrainṪ = 0 is added to the problem and a penalty function P(u; ,σ) is used as smoothing technique, so that the problem can be reformulated as follows min T,v,u 1 0 (−Tv + TP(u; ,σ)) dτ, As in Section 5, P(u; ,σ) is a piecewise quadratic penalty function defined as Now, the Hamiltonian formulation of the problem (14) gives as a result the following BVP where the thrust u, computed by solving the equation Since the problem is highly nonlinear, it is chosen as starting approximation of the solution h = λ 1 = λ 2 = λ 3 = 1, λ 4 = 0, v(τ) = τ(1 − τ), m(τ) = (m(1) − m(0))τ + m(0) and T = 0.01. The choice of a good initial approximation is the main matter when the parameters of the penalty functionσ and become extremely small. In this case, it is helpful to apply a continuation strategy for the parameter , changing the value of this parameter from 0 = 10 −1 to the desired value of . To highlight the advantages of this strategy, we solve the optimal control problem (15) choosingσ = 10 −4 and two different values of = 10 −3 , 10 −6 , the solution is plotted in Figure 5.
In Table 18 the results are computed without applying the continuation strategy, hence we observe that if on one hand only the codes bvp5c, tom and tomc fail for = 10 −3 , on the other all the codes do not converge for = 10 −6 . Consequently, in Table 19 we run the codes using the continuation strategy. All the numerical tests use an initial mesh of 16 equidistant points. For the continuation strategy in Table 19, except for acdc and acdcc, the parameter is initially set to 0 = 10 −1 ( 0 = 1 for tom and tomc), and then it is changed using N = 10 logarithmically equispaced values up to reach the value required . However, to obtain the convergence of bvp4c for = 10 −6 , we put the value of N = 100 when tol = 10 −4 and N = 20 when tol = 10 −6 and for tom/tomc we put the value of N = 55. The conditioning parameters reported in Table 20 show that the problem is not stiff, but it is ill conditioned since κ 1 > κ 2 .

Minimization of the Fuel Cost in the Operation of a Train
As in [2,47] an optimal control problem in transportation is to minimize fuel cost in the operation of a train. To simplify the track is supposed to be straight. Let x be the position along the track measured from a fixed reference point and v the velocity of the train, such that the minimization problem is equivalent to solve the optimal control problem min v,u a 4.8 0 u a v dt, where F(v(t)) models the friction due to the rolling of the wheels and the air resistance and h(x) is the active component of the gravitational force due to hill slopes that are respectively defined as Moreover, the control variables u a and u b represent respectively the acceleration provided by the engine and the deceleration from applying the brakes.

Conclusions
In this paper, after a review of general-purpose codes for solving boundary value problems we have solved some challenging optimal control problems derived using the indirect method. The presented results show that this approach could be a good alternative to the direct methods for the solution of this kind of problems, especially if the mesh selection strategy adopted is suitable for stiff problems in the case of hypersensitive problems, or an appropriate initial condition is computed for the nonlinear iteration using a continuation strategy. All these techniques can sometimes require the application of some regularization procedure, as in the presence of singular arc. Our goal with this paper is to give some indications useful to handle the input parameters of a BVP code to achieve an accurate solution, since the default values assigned usually works for very simple regular problems. Moreover, some codes give in output information about the stiffness and the conditioning of the problems, which could be used in choosing the correct solution method.