HOMOTOPIC APPROACH FOR TURNPIKE AND SINGULARLY PERTURBED OPTIMAL CONTROL PROBLEMS

The first aim of this article is to present the link between the turnpike property and the singular perturbations theory: the first one being a particular case of the second one. Then, thanks to this link, we set up a new framework based on continuation methods for the resolution of singularly perturbed optimal control problems. We consider first the turnpike case, then, we generalize the approach to general control problems with singular perturbations (that is with fast but also slow variables). We illustrate each step with an example. Résumé. Dans cet article nous présentons le lien entre les points de vue turnpike et perturbations singulières en contrôle optimal en montrant notamment que le premier est un cas particulier du second. Ce lien important nous permet de mettre sur pieds une nouvelle approche de résolution numérique, basée sur les méthodes homotopiques, des problèmes de contrôle optimal perturbés singulièrement. Tout d’abord on s’intéresse aux cas turnpike, puis on généralise l’approche au cas perturbations singulières (marqués par la présence de variables lentes); chaque partie étant illustrée par un exemple.


Introduction
Recently, new investigations have been carried out on optimal control problems with turnpike properties [13]: the optimal solution remains for a long time close to the optimal steady-state solution associated to the static optimal control problem and at the beginning and at the end we have a transcient-short time arc. On the other hand in the 1960's lots of works have been done on optimal control problems with singular perturbations [1,[3][4][5]7,9,10,15]. The first result we present here establishes that an optimal control problem with the turnpike property is a singularly perturbed optimal control problem with only fast variables.
Then we are interested in the numerical resolution of such problems by indirect methods. First, we treat the turnpike case: we present the appropriate shooting function, see [13], in which the unknowns are the state and the costate at the middle of the time interval and for which we have a very good initial guess from the static optimal solution. We propose here to improve this procedure with two homotopies: the first one is an homotopy on the two-point condition which connects the static optimal solution, i.e. the zero-order outer solution of the problem (seen as a singularly perturbed control problem), to the original problem for a final time not too long and the second one is on the final time. At the end, it is possible to extend this method to solve singularly perturbed optimal control problems. For this we have to adapt the first homotopy on the two-point condition to take into account the state equation on the fast variables. Several exemples from the literature show the efficiency of our approach.

Turnpike optimal control problem
We consider the following formulation, excerpt from [13], of a turnpike optimal control problem: min T 0 f 0 (y(t), u(t)) dt, T > 0 large enougḣ y(t) = f (y(t), u(t)), y(t) ∈ R m , u(t) ∈ R k , y(0) = y 0 , y(T ) = y f , where y 0 and y f are given values in R n and where f 0 : R m × R k → R and f : R m × R k → R m are supposed to be regular (at least C 2 ). Note that to simplify the presentation we consider the case where y(0) and y(T ) are fixed to given values, but the results can be generalized for boundary conditions of the form R(y(0), y(T )) = 0. The associated static optimal problem is then defined as We assume that we are in the normal case for both the original and the static optimal problems and that the static optimal problem has a unique solution denoted (ȳ,ū,q), whereq is the associated Lagrange multiplier. We say that we have the turnpike property if there exists an extremal solution of (OCP T ) denoted (y T (·), u T (·), q T (·)) which remains most of the time close to the static solution (ȳ,ū,q). More precisely, it is expected that if T is large enough then the extremal is approximately made of three pieces: • the first is a short-time arc defined on [0, τ ] for a τ > 0 along which the extremal gets closer to (ȳ,ū,q); • the second is quasi static and remains close to the steady-state (ȳ,ū,q) over a long-time interval [τ, T −τ ]; • the third piece is a short-time arc defined on [T − τ, T ] which joins the final condition.
In [13], the authors prove that under suitable assumptions there exists constants C 1 > 0 and C 2 > 0 and a time T 0 > 0 such that for every T > T 0 , the optimal control problem (OCP T ) has at least one solution having a normal extremal lift (y T (·), u T (·), q T (·)) satisfying for every t ∈ [0, T ].

Links between turnpike and singularly perturbed optimal control problem
Let us reparameterize the time setting s = εt (with ε = 1/T small since T is supposed to be large). The problem (OCP T ) is then equivalent to the following problem: which has the form of a singularly perturbed optimal control problem with only fast variables [9]. The Pontryagin's Maximum Principle [11] gives us some necessary conditions of optimality in the cotangent bundle. We suppose here that the maximisation of the pseudo-Hamiltonian furnishes us the control as a smooth function of the state and costate u(y, q). This is in particular the case if we suppose the strict Legendre condition. Then the PMP implies the singularly perturbed boundary value problem of the form: εẏ(s) = g(y(s), u(y(s), q(s))), εq(s) = −∇ y H(y(s), u(y(s), q(s)), q(s)), where the pseudo-Hamiltonian H is H(y, q, u) = q · g(y, u) − f 0 (y, u). We define the following reduced problem:    0 = g(y(s), u(y(s), q(s))), 0 = −∇ y H(y(s), u(y(s), q(s)), q(s))). whose solution, denoted (ȳ,q), is the so-called zero order outer solution of (BV P T ε ). In this setting, denotinḡ u := u(ȳ,q), the triplet (ȳ,ū,q) satisfies the Karush-Kuhn-Tucker necessary optimality conditions of the static problem (SOCP T ).
We introduce the notation β := (y, q) and G := (g, −∇ y H) such that the differential equation in (BV P T ε ) may be written in the more compact form εβ(s) = G(β(s)). With this notation and under suitable assumptions, we can prove that the fast variable β ε , solution of (BV P T ε ), converges on every compact interval [a, b] ⊂ [0, 1], uniformly to the zero order solutionβ = (ȳ,q) [7,8,14]. Moreover, in order to get a uniform convergence on the fast variables on [0, 1] (of course it's trivial to see that boundary conditions on the fast variables are not necessary verified by the approximate zero order solution), one can add some boundary layer conditions on the fast variables in order to make a correction at the initial and final conditions. Roughly speaking (see [7,8,10] for more details), at the zero order one can write the fast variables in the form: with β i (resp. β f ) → 0 when τ = t/ε (resp. σ = (1 − t)/ε) → ∞ and where the zero order boundary layers β i and β f satisfy: To summarize, one can say that the boundary layers correspond to the first and third short-time pieces of the solution of the turnpike problem and that the zero order outer solution is trivially the static solution. Finally we make the following crucial observation: Proposition 2.1. A turnpike optimal control problem (i.e. an optimal control problem in the form of (OCP T )) is equivalent to a singularly perturbed optimal control problem with only fast variables (in the form of (OCP T ε )).

Numerical methods for turnpike control problems
From the numerical point of view, solving a singularly perturbed boundary value problem is difficult due to the stiffness of the dynamical system. For an initial value problem it is known that some Runge-Kutta schemes as Radau methods (which is an implicit method) are more suitable than others to integrate stiff differential equations [6]. When we are dealing with boundary value problems, one can use a collocation method on a fully discretized grid and refine the mesh iteratively where the stiffness introduces large errors [4]. The main idea is to integrate on several small time intervals and match the extremities since the stiffness amplifies errors during the integration. The multiple shooting methods hinge on this idea. Let us recall the principle of the simple shooting method, considering (BV P T ε ), before giving some details about multiple shooting. Let us ε ) with initial condition (y(s 0 , s 0 , α, β), q(s 0 , s 0 , α, β)) = (α, β) and introduce the mapping S(β) := y(1, 0, y 0 , β) − y f . The simple shooting method consists in solving S(β) = 0 with a Newton method, the point (y(1, 0, y 0 , β), q(1, 0, y 0 , β)) being computed with a numerical integrator. The idea of multiple shooting is to add intermediate steps.
We propose here a new approach based on the idea of multiple shooting to improve numerical stability and a kind of automatic mesh refining like algorithm to improve convergence of the underlying Newton method, known to be sensitive to the initial guess. The orginality of this approach is double: the multiple shooting method is a variant of the classical one that we have presented just before, which has been introduced in [13] and which is appropriate when dealing with turnpike control problems. Secondly, to improve the convergence we introduce a two-steps homotopic method, the homotopic algorithm being a differential path following technique. The homotopy method is a Predictor-Corrector method with arc length parameterization, with a high-order variable step-size Runge-Kutta scheme for the prediction and a simple modified Newton method for the correction.
The shooting function. When dealing with turnpike boundary value problems (that is boundary value problems coming from turnpike optimal control problems), it is convenient to use multiple shooting. As introduced in [13], it is better to add an intermediate step at the middle of the time interval, denoted s mid := 1/2, since we know that the extremal trajectory is close to the static solution (ȳ,q,ū) at this time. This information gives a very good initial guess if the unknown of the shooting method is the point at s mid . The shooting function is thus defined by: There are two halves integration each starting from the middle point, note that the first one is computed backward to match the initial condition while the second one is computed forward.
The shooting homotopic function on the boundary conditions. Since the static optimal point (ȳ,ū) is the equilibrium point that optimizes the cost function, if we consider the turnpike optimal control problem (OCP T ε ) where we replace the initial and final conditions by y(0) =ȳ and y(1) =ȳ, then, the optimal solution is simply the constant extremal trajectory (y, u) ≡ (ȳ,ū) with associated adjoint vector q ≡q, and this whatever the final time T . The idea is thus to connect the solution of this simple control problem to the solution of (OCP T ε ) by performing an homotopy on the boundary conditions. We thus introduce the following homotopic shooting function (still denoted S): The steps of our algorithm. Let us recall that the main goal is to solve (BV P T ε ) for an ε small enough to have the turnpike property and so possibly we have to deal with singular perturbations (the dynamical system may be stiff). To tackle this kind of optimal control problems, we propose the following approach: (1) First step: solve the static problem (SOCP T ) and get (ȳ,q), solution of S(y mid , q mid , 0) = 0.
(2) Second step: fix ε to a relatively large value to reduce stiffness and solve S(y mid , q mid , λ) = 0 by homotopy for λ ∈ [0, 1], starting from the solution (ȳ,q) at λ = 0. This step is refered as the homotopy on the boundary conditions. (3) Third step: make a homotopy on the parameter ε to the desired value (supposed to be small).

Example of turnpike control problem and numerical results
Let us consider the following turnpike optimal control problem from [13]: , with y := (y 1 , y 2 ), ε = 1/T and T = 70. The associated boundary value problem is where q := (q 1 , q 2 ) and the static optimal solution is given by (ȳ 1 ,ȳ 2 ,q 1 ,q 2 ,ū) = (2, 0, −1, −1, 1). For the second step, we fix T = 40. Figure 1 shows the evolution of the optimal extremal trajectory along the homotopy on the boundary conditions, while Figure 2 is the same but for the third step. We compare our approach with the variant of the multiple shooting method presented before taking the static solution as the initial guess and we compare also this approach with the reduced procedure where we only perform the first and second steps but with the desired value for the parameter ε. The results of these numerical experiments are given in Table 1.

Remark 2.2.
We use the software HamPath [2] to compute the paths of zeros of the homotopies and to solve the shooting equations for these numerical experiments and the followings in the rest of the paper.    Table 1. Comparison of the numerical methods (for the 3 strategies the initial guess is the solution of the static problem (step 1)): shooting from the middle point, our approach with only step 2 and our complete approach with step 2 for ε = 1/40 and step 3. The numerical integrations are done with the dopri5 code with relative and absolute local errors of 1e-8 and 1e-14. The signification of the symbols are : : the algorithm converges without difficulty; H : the algorithm converges with warnings from the numerical integrations, or is very slow; : the algorithm diverges.
3. Generalization to singularly perturbed optimal control problems
We can thus rewrite the reduced problem in the following equivalent form: This reduced boundary value problem is called the associated zero order boundary value problem of problem (BV P ε ) and will be the initialization of our homotopy, that is of the first step of our algorithm.
Remark 3.1. We can of course simplify the problem (BV P 0 ) by eliminating the β variable if it is possible to obtain the analytic expression of the implicit function Φ as it is the case in the simple example we treat in the next paragraph. Otherwise we have to numerically solve the implicit equation to obtain the boundary condition of the problem (BV P 0 ) or to modify the shooting function by replacing β 1,··· ,m (0) = Φ 1,··· ,m (ψ(0)) (respectively β 1,··· ,m (1) = Φ 1,··· ,m (ψ(1))) by the equation 0 = G(ψ(0), β(0)) (respectively 0 = G(ψ(1), β(1))). There is a second reason for using the differential equation in β : the solution of the (BV P 0 ) problem is exactly the same as the initial problem for the first homotopy we use in our method, as it is explained in the following.
Shooting function and shooting homotopic function. The aim remains to solve (BV P ε ) for ε sufficiently small and the shooting function is still defined from the middle point. The most significant difference here compared to the turnpike case may be seen in the homotopic function since we have to take into account the new boundary conditions and the new dynamics (of the fast variables) from the zero order boundary value problem. We thus define the following homotopic function: with initial condition (ψ 0 , β 0 ) at time t 0 .
The steps of the algorithm. As in the turnpike case, we proceed in three steps: (1) First step: solve the zero order boundary value problem (BV P 0 ) to get the zero order outer solution, solving S(y mid , q mid , 0) = 0. (2) Second step: fix ε to a relatively large value to reduce stiffness and solve S(y mid , q mid , λ) = 0 by homotopy for λ ∈ [0, 1], starting from the solution at λ = 0. This corresponds to make a continuation from problem (BV P 0 ) to (BV P ε ). (3) Third step: make a homotopy on the parameter ε to the desired value (supposed to be small). Let us emphasize the fact that compare to the turnpike case, the first homotopy is done on the differential equations as well as on the boundary conditions.

Example and results
Let us consider the following optimal control problem [9]: The associated boundary value problem obtained from the maximum principle and its associated zero order boundary value problem are given by: , y(0) = 0, q(t) = 2x 4 (t), . Figure 3 and 4 show the results we obtained for the first, respectively second, homotopy. In Table 2 we make a comparison of three different approaches: (i) simple shooting from the middle point fixing λ = 1, (ii) simple shooting followed by the homotopy on λ (that is steps 1 and 2 of our algorithm) (iii) simple shooting followed by the two homotopies on λ and then ε (that is steps 1, 2 and 3).

Conclusion and perspectives
In this work, we have presented a new approach for the numerical resolution of singularly perturbed problems in optimal control, this has been possible thanks to the link that we have established between the turnpike and singular perturbations frameworks. It must be said that the idea of this homotopic method is well justified by the fact that the two problems (original problem and reduced problem) are close (in the sense of ε) and the reduced problem is much easier to solve. This paper has been prepared with the aim of laying the groundwork for a more thorough study. We therefore propose to carry out a detailed analysis of the different integrators adapted to stiff problems; then, to apply  Algorithms ε = 0.09 ε = 0.08 ε = 0.06 ε = 0.027 ε = 0.025 Shooting from the middle point H Step 2 Steps 2-3 Table 2. Comparison of the numerical methods (for the 3 strategies the initial guess is the solution of the static problem (step 1)): shooting from the middle point, our approach with only step 2 and our complete approach with step 2 for ε = 1/40 and step 3. The numerical integrations are done with the dopri5 code with relative and absolute local errors of 1e-8 and 1e-14. The signification of the symbols are given in Table 1.
our approach with the most efficient integrator on several case studies and make a comparison with existing methods for solving stiff problems (considering direct and indirect methods).