ON NUMERICAL METHODS FOR SINGULAR OPTIMAL CONTROL PROBLEMS: AN APPLICATION TO AN AUV PROBLEM

. We discuss and compare numerical methods to solve singular optimal control problems by the direct method . Our discussion is illustrated by an Autonomous Underwater Vehicle (AUV) problem with state constraints. For this problem, we test four diﬀerent approaches to solve numerically our problem via the direct method. After discretizing the optimal control problem we solve the resulting optimization problem with (i) A Mathematical Pro- gramming Language (AMPL), (ii) the Imperial College London Optimal Control Software (ICLOCS), (iii) the Gauss Pseudospectral Optimization Software (GPOPS) as well as with (iv) a new algorithm based on mixed-binary nonlinear programming reported in [7]. This algorithm consists on converting the optimal control problem to a Mixed Binary Optimal Control (MBOC) problem which is then transcribed to a mixed binary non-linear programming problem (MBNLP) problem using Legendre-Radau pseudospectral method. Our case study shows that, in contrast with the ﬁrst three approaches we test (all relying on IPOPT or other numerical optimization software packages like KNITRO), the MBOC approach detects the structure of the AUV’s problem without a priori information of optimal control and computes the switching times accurately.

Abstract. We discuss and compare numerical methods to solve singular optimal control problems by the direct method. Our discussion is illustrated by an Autonomous Underwater Vehicle (AUV) problem with state constraints. For this problem, we test four different approaches to solve numerically our problem via the direct method. After discretizing the optimal control problem we solve the resulting optimization problem with (i) A Mathematical Programming Language (AMPL), (ii) the Imperial College London Optimal Control Software (ICLOCS), (iii) the Gauss Pseudospectral Optimization Software (GPOPS) as well as with (iv) a new algorithm based on mixed-binary nonlinear programming reported in [7]. This algorithm consists on converting the optimal control problem to a Mixed Binary Optimal Control (MBOC) problem which is then transcribed to a mixed binary non-linear programming problem (MBNLP) problem using Legendre-Radau pseudospectral method. Our case study shows that, in contrast with the first three approaches we test (all relying on IPOPT or other numerical optimization software packages like KNITRO), the MBOC approach detects the structure of the AUV's problem without a priori information of optimal control and computes the switching times accurately.
1. Introduction. This paper focuses on numerical solution of singular optimal control problems (SOCPs) via the direct method. We test and compare four different approaches applied to a problem for Autonomous Underwater Vehicle (AUV's). Our illustrative problem is a singular optimal control problem (SOCP) involving a simplified model for the path planning of an AUV's in a horizontal plane treated previously in [16] and [5].
In SOCPs the control appears linearly in the dynamics and cost and it has lower and upper bounds. It is well known that the optimal control for those problems is a concatenation of bang-bang or bang-singular type control [14,13,8]. The main difficulty when solving singular optimal control problems numerically lies in the determination of the switching structure of the optimal control, i.e., the sequence of singular and bang-bang sub-arcs composing the optimal trajectory and location of switching points [13,9].
In this paper, we confront and compare four different approaches based on the direct method to solve numerically this type of problems 1 . The first two methods make use of the nonlinear optimization package IPOPT [21] to solve the transcribed nonlinear problem. In the first method, we interface IPOPT with AMPL [10,18] using the Implicit Euler Method [2]. For the second method, we use ICLOCS [19]. The other two methods are of a different nature. One is the GPOPS [15] and the other is a recently developed method, proposed in [7], based on Mixed-Binary Non-Linear Programming (MBNLP) method.
The MATLAB based program ICLOCS is designed to solve optimal control problems with general path and boundary constraints and free or fixed final time. The non-linear programming problem (NLP) obtanied by discretization is solved by IPOPT; see [19]. The GPOPS, also a MATLAB based program, uses the Gauss pseudospectral method to transcribe the original problem to a nonlinear programming problem (NLP) which is then solved either by SNOPT or IPOPT [15].
The MBNLP algorithm described in [7] is a direct approach divided into two stages: first the original optimal control problem is converted to a mixed binary optimal control problem (MBOC) which is then transcribed as a mixed binary non-linear programming (MBNLP) using the Legendre-Gauss-Radau (LGR) pseudospectral method; see description in [7]. The control structure is not assumed a priori but it is determined automatically and it is designed so that the switching times are computed accurately. Firstly, considering the number of switching times finite and the switching times as decision variables, the SOCP is reformulated as a multi-domain optimal control problem. Next, taking advantage of feedback law provided by the Pontryagin's Minimum Principle, binary variables are introduced and the problem is rewritten as a mixed-binary optimal control problem. Finally, the Legendre-Gauss-Radau pseudospectral method is applied to discretize such problem leading to a MBNLP that is then solved by adequate optimization software packages (here we use KNITRO) yielding the structure of optimal control and the position of switching points for the original SOCP. The Legendre-Gauss-Radau pseudospectral method approximates the state and control variables by interpolating polynomials with Legendre-Gauss-Radau (LGR) points as collocation points [12]. Then, collocating the state equations and path constraints and using a differentiation matrix, the problem is transcribed to a nonlinear programming problem, which, in turn, can be solved by well-developed parameter optimization algorithms.
To discuss and compare these four different methods, we apply each of them to an AUV's problem treated in [5,16]. We show that all four methods capture the same structure for the optimal control. However, ICLOCS and GPOPS exhibit some oscillations. Our analysis also shows that the MBNLP method reported in [7] calculates the switching times with high accuracy and noteworthy, with a small number of collocation points. Here, we do not dwell on descriptions of the Implicit Euler method, ICLOCS and GPOPS methods since they are well established in the literature. In contrast, we give detailed explanation of the use of the MBNLP method for AUV problem, a singular optimal control problem with state constraints. We also refer the reader to [5], where a analytic analysis of the AUV's problem is presented as well as partial validation of the numerical solution computed using AMPL interfaced with IPOPT.
We choose the AUV problem to illustrate the MBNLP method because, as reported in [5], this is a problem with state constraints and the solution exhibits singular controls and boundary controls.
The paper is organized as follows. Section 2 is devoted to the formulation of the AUV's problem and a brief discussion of the Pontryagin's minimum principle applied to it. In Section 3, after some background necessary for understanding the LR pseudospectral methods, we describe the application of the MBNLP method to the AUV's problem. The numerical results and the accuracy of the methods are reported in Section 4. There, the results of this method are also compared to the result of ICLOCS and GPOPS and Implicit Euler Method. This paper ends with the conclusions presented in Section 5.
2. The AUV optimal control problem. To illustrate the application of different numerical schemes to SOCPs, we consider the problem of determining the path of an Autonomous Underwater Vehicle (AUV) in a horizontal plane to drive the vehicle from an initial point point to a target set T in the minimum time when currents are taken into account. This problem was first considered in [16] and later on in [5]. We observe that, as in [5], we focus our attention to the motion of the AUV in a horizontal plane using a simplified point mass model. The problem of interest is where t f is the minimum time to be determined, denotes the space of absolutely continuous functions from [0, t f ] to R 4 and L ∞ ([0, t f ]; R 2 ) the space of essentially bounded functions from [0, t f ] to R 2 . A special feature of (P ) is that it involves one state constraint , defined as h(x, y, φ, u) = u − 2.
The pair (x, y) is the position of the vehicle in the plane, φ is the heading angle of the vehicle, u is the velocity of the vehicle, r the angular velocity and f is the thruster force. The term −Ku(t)|u(t)| depicts the drag force and we assume throughout this paper that K = 1 . The aim is to determine the minimum time t f needed to take the vehicle from an initial position (x 0 , y 0 ) = (40, −2) to a target set T . We consider that the velocity v of ocean currents is known with components merely on the x but depending on the y position: v(t) = (0.8tanh(y(t)), 0).
The initial and final configurations of the vehicle are considered to be Mimicking what is done in [5], we assume the target set T to be a small ball around the origin of the horizontal plane: Observe that comparing our problem with that in [5], it is a simple matter to see that we drop the state constrains −u(t) ≤ 0. This is because since we have u(0) = 0 and u(t f ) = 0 the constraints −u(t) ≤ 0 is active at the initial and final time and, given the nature of our problem, it is forseen that [5] is redundant. Dropping such constraint from the definition of our problem does not change any of the numerical findings and, besides, save us from considering some technicalities concerning the nature of necessary conditions, as, for example, any possible degeneracy phenomena. Before engaging in the discussion of necessary conditions for (P ), it is worth mentioning that we deduce from the physical meaning of our problem the existence of an optimal solution for (P ). However, it is also possible to see that (P ) satisfies known conditions under which existence of solution is guaranteed. In this repect, a possible reference is [20].
Discussion and analysis of necessary conditions for (P ), following [20] can be found in [5], where the extra state constraint u(t) ≥ 0 is considered. The removal of this constraint does not change such analysis. Here we differ from what it is done in [5] since we use the Pontryagin's Minimum Principle as it appears in [17] (see also [14]). This is because the latter is better tailored for our purpose. However, application of the Pontryagin's Minimum Principle, requires that some regularity conditions are satisfied. Next, we check they are.
Indeed, the state constraint is of order one (see, for example, [14]). To see this, take any feasible solution of (P ) such that h(x(t), y(t), φ(t), u(t)) = 0 in a boundary interval t ∈ [t b 1 , t b 2 ], where thet b 1 is the entry time and t b 2 is the exit time, these points also known as junction points. Then, for t ∈ [t b 1 , t b 2 ], we have u(t) = 2 and proving that state constraint is or order one an, moreover, that the boundary control f (t) = 4 is a smooth control taking values in the interior of the control set: f (t) = 4 ∈]−5, 5[. Under these two regularity conditions (order one of the state constraints and interior boundary controls) we are then in position to apply the Pontryagin's Minimum Principle. Define the Hamiltonian for (P ) to be where, λ = [λ x , λ y , λ φ , λ u ] denotes the adjoint variable. Adding the multiplier µ associated to the state constraint u − 2 ≤ 0 to the Hamiltonian, we get the augmented Hamiltonian: If (x,ū, t f ) is a local optimal solution to (P ), the Pontryagin's Minimum Principle asserts the existence of adjoint functions λ x , λ y , λ φ , λ u and a multiplier µ satisfying (here we drop the dependence on t since this is clearly understood) together with the usual Minimum Condition on the Hamiltonian transversality conditions (here, for a closed set A ⊂ R n , N A (x) denotes the normal cone of convex analysis), jump conditions at contact or junctions points t 1 with the boundary of the state constraints The switching function is and, from the minimization of the Hamiltonian (5), we get the following law: where r sing (t) and f sing (t) take values in the interior of respective the control set: Let us assume that there is a time interval [t s 1 , t s 2 ] in which r takes singular values.
and λ x (t f ) = 0, (quite reasonable assumptions taking into account the physical interpretation of our problem and in agreement with the numerical findings in [5]), we get from d dt σ f (t) = 0: ) .

Using this expression in
and ∂ ∂r Recall that we require that (x(t f ),ȳ(t f )) ∈ T . Given the current and the initial position of the vehicle we expect x(t f ) < 0 again, this is in agreement with the numerical findings reported in [5]). From the transversality conditions we deduce that λ x (t f ) ≤ 0 and it follows from (7) that the generalized Legendre-Clebsch condition is satisfied. It follows from (6) that on any possible singular interval [t 1 , t 2 ] the singular value of r is given by r = r sing (x) = −0.8 cos 2 (φ) cos 2 (φ) cosh 2 (y) .
We now turn to f . We know that if the state constraint is active on a boundary interval, then f is a boundary control taking singular values: and the boundary interval is also a singular interval. We can not exclude from the analysis the existence of any singular control f outside a boundary interval. However, the numerical findings as reported in [5] indicate that there is no such singular interval for f . Moreover, the physical interpretation of our problem does also point out to the exclusion of such situations. We then proceed considering that only boundary controls for f are present.
Here, and as in most of the practical problems, we assume that the number of switching points of the controls is finite. Although we do not prove analytically that this is the case, this is, once again, to be expected by the physical meaning of our problem.
3. Mixed-binary non-linear methods. Before presenting the MBNLP method applied to the AUV's problem, some preliminaries are called for.

Preliminary considerations.
For the Legendre-Gauss-Radau pseudospectral method applied to optimal control problems [11] the first n nodes ξ i are Legendre-Gauss-Radau (LGR) nodes and the last node is selected as ξ n = +1. The Legendre-Gauss-Radau nodes are the roots of P n−1 (τ ) + P n (τ ), where P n (τ ) is the well-known Legendre polynomial of degree n.
Define j (t), j = 0, · · · , n to be where j (t) are the Lagrange polynomials based on LGR nodes. It is an easy task to see that the following Kronecker property holds Consider any differentiable function g : [−1, 1] → R approximated by with derivative approximated by: Since˙ i (τ ) is a polynomial of degree n, we have: Using the above two equations we get: is the (i, j)-th component of a matrix D, called the differentiation matrix in [6]. According to (11), the entries of differentiation matrix D are computed by taking the analytical derivative of i (τ ) and evaluating it at the collocation points ξ j for i, j = 0, . . . , n. More computationally practical methods for deriving these entries, in an accurate and stable manner, can be found in [23,22,1].

3.2.
Application of the MBNLP to (P ). We now turn to the description of the MBNLP method as it appears in [7]. This is done in two stages. First, using the information provided by the Pontryagin's Minimum Principle on the structure of the optimal controls and considering that the number of switching points s of the controls to be finite, we convert (P ) to a fixed time binary optimal control problem.
To fully clarify this procedure, we do such conversion in three steps. In stage one, using the structure of the controls, we divide the time interval [0, t f ] into a finite number of subintervals and we consider a multi-domain optimal control problem. Then we introduce binary variables and we get a multi-domain mixed binary optimal control. This problem is then converted into a fixed time binary optimal control problem using the known time transformation, a conversion needed because of the nature of the Legendre Guass Radau Pseudospectral method applied in the stage two, where the discretization of resulted binary optimal control problem following the LGR Pseudospectral method described above. This discretization leads to a MBNLP problem. The solution of MBNLP problem is an approximate solution of the original problem (for a fully description see [7]). Next we apply such a procedure to the AUV problem (P ).

3.2.1.
Conversion of AUV's problem to a binary optimal control. As stated before, let us sssume that the optimal control problem AUV's has a finite number of switching points and take s to be the upper bound of such number. We consider t 1 , . . . , t s as candidates for switching points to be determined, where The s points t 1 , . . . , t s break the total interval [0, For k = 0, . . . , s, let be the restriction of x(t) and u(t) to the k-th subinterval [t k , t k+1 ]. Then, the control and state functions can be expressed as: To do so, we introduce binary variables as follows. In each subinterval [t k , t k+1 ], we consider two unknown vector with binary components. The control functions are then expressed as Indeed, if µ 2 = 1 to be 1 at the same time, we add the following conditions: Putting all together we now convert our problem into the following multi-domain mixed-binary optimal control problem min ( all the above hold for t ∈ [t k , t k+1 ] and k = 0, . . . , s) and, for all k = 0, . . . , s − 1 Observe that (16i)-(16l) enforce the continuity of state functions at the switching points.
Taking into account the conclusions of the Pontryagin's Minimum Principle of the previous section (recall that we assume the number of switching times is finite and bounded by s) it is a simple matter to see that this new problem is equivalent to (P ).
Before engaging into the discretization using the LGR Pseudospectral method, we need to convert the above problem into a fixed time optimal control problem. We do so by mapping each time sub-interval [t k−1 , t k ], k = 0, . . . , s to [−1, 1] using the time transformation: Note that dt With the time transformation (17) we get a fixed time optimal control problem in the time domain [−1, 1]: 2 )), 2 ) ≤ π, (18j) We observe that the time transformation (17) is a known procedure to convert free time problems (see, for example, [20]) into an equivalent fixed time problem.

3.2.2.
Discretization of the mixed binary optimal control. We can now proceed with the discretization. Considering (9), the states x [k] , y [k] φ [k] and u [k] are approximated by Lagrange polynomials as where for k = 0, . . . , s, i = 0, . . . , n, the coefficients a k i , b k i , c k i and e k i are unknown and By replacing approximations (19) and (20) in the dynamics equations of the fixed time optimal control problem and collocating it in LGR points ξ j , j = 0, . . . , n − 1, we get: 2 ) − e k j .|e k j |) = 0 Finally, the fixed time optimal control problem is converted to the following finitedimensional MBNLP problem corresponding to the Radau pseudospectral method: 4. Numerical results. We now compare the numerical results given by the method described in Section 3.2 with those provided by classical methods using AMPL with implicit Euler method, ICLOCS and GPOPS. We start by reporting on the numerical solution of AUV's problem given by AMPL with implicit Euler method. We consider N = 10000 grid nodes. AMPL (see [10]) calls the solver IPOPT (see [21]) with acceptable tolerance of 10 −8 . The control function obtained by this method is shown in Figure 4.
We also use ICLOCS Version 1.2.1 [19] with N = 10000 grid points interfacing with IPOPT with the same acceptable tolerance 10 −8 . Moreover, we use GPOPS Version 4.x [15] with three phases and in each phase N = 40, again, with the same acceptable tolerance 10 −8 . The control functions computed by these two methods are shown in Figure 2 and Figure 3.
As clearly shown in Figures 2 and 3, the control functions exhibit the Gibbs phenomenon (with respect to Gibbs phenomenon we refer the reader to, for example, [4]). This is due to the fact that the functions are discontinuous and both ICLOCS and GPOPS present some difficulties dealing with such behavior.   We now turn to the the method in Section 3.2. We implement it in Matlab on a personal computer. The final BINLP (21) is solved with the help of AMPL, implemented in Matlab, interfaced with Knitro (for a description of Knitro see [3]). Using the method in Section 3.2 and to detect the structure of the optimal control, we first set s = 5 and n = 20. The resulting switching points and the  control function are  and the controls are shown in Figure 4. Analyzing these data, we see that t 1 = 4.913684e − 02 is reported as switching point but for the control function r, the switching point t 2 = 4.922503e − 02 is very close to t 1 and f does not switch in t 1 .
Then we conclude that this point t 1 is not switching point and it can be removed. Accordingly, the exact number of switching points is s = 4 and not s = 5. We present the computed switching times for various values of n with s = 4 in the Table 1, indicating that there is no need to go for values above n = 12.
The control and state functions, computed with s = 4 and n = 14, are plotted in Figure 5 and 6 .
In Table 2, we summarize the main results of the AUV's problem obtained by the four numerical methods.  The MBNLP method described in section 3.2 detects the structure of the AUV's problem without a priori information of the optimal control and computes the switching times accurately using a small number of collocation points. ICLOCS, GPOPS and AMPL with Implicit Euler method can find the structure of optimal control. However, the optimal control computed by ICLOCS and GPOPS exhibits oscillations and low precision. Noteworthy, MBNLP method gives direct information about the switching times whereas the other three methods tested fail to do so.

5.
Conclusion. In this paper, we compared various known numerical methods to solve SOCPs with the recent developed MBNLP method [7] using an AUV's problem with state constrains. We see that all of them detected the structure of the AUV's problem without a priori information of optimal control but the computed control functions in ICLOCS and GPOPS have oscillations. Remarkably, the optimal control calculated by the MBNLP method and AMPL with implicit Euler method have higher accuracy than the other two methods. The main advantage of the the MBNLP method is its ability to compute accurately the switching times and, moreover, it does that with a small number of collocation points. The other methods fail to provide direct information about such points, so needed when dealing with SOCP. This is because the MBNLP method described in section 3.2 is designed for singular optimal control problems while the other three methods can be applied to more general problems.
The AUV's problem has controls with singular values given by the feedback laws derived in section 2. It is however known that the singular values for SOCP may not be known or, if known, they may depend on the multipliers. The MBNLP method reported in [7] can be generalized to cover such cases; illustration of this method to problems with singular arcs that depend on the multipliers will be reported elsewhere.