ADAPTIVE TIME–MESH REFINEMENT IN OPTIMAL CONTROL PROBLEMS WITH STATE CONSTRAINTS

. When using direct methods to solve continuous-time nonlinear optimal control problems, regular time meshes having equidistant spacing are most frequently used. However, in some cases, these meshes cannot cope accurately with nonlinear behaviour and increasing uniformly the number of mesh nodes may lead to a more complex problem. We propose an adaptive time–mesh reﬁnement algorithm, considering diﬀerent levels of reﬁnement and several mesh reﬁnement criteria. Namely, we use information of the adjoint multipliers to decide where to reﬁne further. This technique is here tested to solve two optimal control problems. One involving nonholonomic vehicles with state constraints which is characterized by having strong nonlinearities and by discontinuous controls; the other is also a nonlinear problem of a compartmental SEIR system. The proposed strategy leads to results with higher accuracy and yet with lower overall computational time, when compared to results obtained by meshes having equidistant spacing. We also apply the necessary condition of optimality in the form of the Maximum Principle of Pontryagin to characterize the solution and to validate the numerical results.


1.
Introduction. Over the past decades, direct methods have become increasingly useful when computing the numerical solution of nonlinear optimal control problems (OCP) [1]. These methods directly optimize the discretized OCP without using the maximum principle and they are known to provide a robust approach for a wide variety of optimal control problems.
In a direct collocation method, the control and the state are discretized in an appropriately chosen mesh of the time interval. Then, the continuous-time OCP is transcribed into a finite-dimensional nonlinear programming problem (NLP) which can be solved using widely available software packages [13]. Most frequentely, in the discretization procedure, regular time meshes having equidistant spacing are used. However, in some cases, these meshes are not the most adequate to deal with nonlinear behaviours. One way to improve the accuracy of the results, while maintaining reasonable computational time and memory requirement, is to construct a mesh having different time steps. The best location for the smaller steps sizes is, in general, not known a priori, so the mesh will be refined iteratively.
In a mesh-refinement procedure the problem is solved, typically, in an initial coarse uniform mesh in order to capture the basic structure of the solution and 4554 LUÍS TIAGO PAIVA AND FERNANDO A.C.C. FONTES of the error. Then, this initial mesh is repeatedly refined according to a chosen strategy until some stopping criterion is attained.
Several mesh refinement methods employing direct collocation methods have been described in the recent years. In [1] and [3] a mesh refinement procedure involving an integer programming technique is developed for changing the discretization in order to improve the accuracy of the approximation. In [19] a density function is used to generate a fixed-order mesh on which the problem is solved. In [16] an approach based on varying the order of the approximation in each mesh interval and using the exponential convergence rate of a Gaussian quadrature collocation method is reported. In [14] and [15], an adaptive mesh refinement strategy based on block-structured refinement method for solving continuous-time nonlinear OCP is presented. It is a purely direct method approach and the convergence is achieved by increasing the number of nodes and by selecting their placement according the refinement criterion. In this algorithm, just one refinement criterion can be defined and it coincides with the stopping criterion.
In this paper, we propose an adaptive mesh refinement algorithm for OCP. This algorithm improves on the ideas advanced in [14,15]. There are three new main features in this algorithm: (a) we introduce several levels of refinement, obtaining a multi-level time-mesh in a single iteration -such concept allows us to implement the nodes collocation in a faster and clever way; (b) we also introduce different refinement criteria -the relative error of the primal variables, the relative error of the dual variables and a combination of both criteria can be used in the refinement procedure; (c) we consider distinct criteria for refining the mesh and for stopping the refinement procedure -the refinement strategy can be driven by the information given by the dual variables and it can be stopped according to the information given by the primal variables. Here, we show that there are advantages in choosing the error of the adjoint multipliers as a refinement criterion, namely in lowering the computational time. This is because the adjoint multipliers are solution to a linear differential equation system, easily solved with high accuracy, and because they give sensitivity information. It is still a direct method approach, although we use some information given by necessary conditions, namely, the adjoint differential equation, which is used in the indirect methods context. To decrease the overall computational time, the solution computed in the previous iteration is used as a warm start in the next one, which proved to be of major importance to improve computational efficiency.
In order to validate the results obtained for this problem, we apply the Maximum Principle of Pontryagin. This analysis allows us to characterize the optimal trajectory and control and it also provides us the differential equation system for the multipliers. This last information is needed to estimate the error in the multipliers for one of the refinement criteria.
This technique is applied to solve two different problems, both with pathwise state constraints: (P 1 ) which involves nonholonomic systems [9], namely a car-like system problem [14], and (P 2 ) involving a compartmental SEIR model [4,10,12] which describes the spreading of an infectious disease. The results show solutions with higher accuracy obtained in a overall computational time that is lower when compared to the ones computed with a mesh having equidistant spacing and to the ones computed with a mesh designed as suggested in [14].
When using this strategy, an OCP is solved using less nodes in the overall procedure which reverts in significant savings in memory and computational cost. In addition, there is no need to decide a priori an adequate mesh meeting our accuracy specifications, which is a major advantage of this procedure. With this mesh-refinement strategy, nonlinear OCP solvers can be used in real-time optimization, in particular in model predictive control [6], since approximate solutions can be produced even when the optimizer is interrupted at an early stage.
This paper is organized as follows. Section 2 states the general form of the optimal control problem considered as well as the corresponding necessary conditions of optimality in the form of a Maximum Principle. Section 3 describes the proposed time-mesh refinement strategy and the algorithm to implement it. Both problems (P 1 ) and (P 2 ) are described in section 4 along with the numerical results highlighting the features of the algorithm proposed. In section 5, we apply the necessary condition of optimality to characterize the solution and to validate the numerical results. We end with some concluding remarks.
2. The optimal control problem. Let us consider the following optimal control problem, in Bolza form, with input and state constraints [17]: • the input constraints • the pathwise state constraint

and
• the end-point constraints In this paper we use a smooth version of the Maximum Principle for state constrained problems which is valid under the following hypotheses. There exists a scalar ξ > 0 such that: [t 0 , t f ]; H4. G is continuously differentiable onx + ξ B; H5. U is compact; H6. h is continuously differentiable onx + ξ B.
Here B denotes the closed unit ball.
A feasible process (x * , u * ) is a W 1,1 local minimizer if there exists δ > 0 such that (x * , u * ) minimizes J(x, u) for all feasible processes (x, u) which satisfy x − x * W 1,1 ≤ δ [17]. Let us consider the Hamiltonian • the transversality condition • the Weierstrass condition • the complementary slackness condition and q : [t 0 , t f ] → R n is a normalized function with bounded total variation defined as We also use a direct method approach where the control and the state are discretized in a chosen mesh in the time interval. Then, the continuous-time optimal control problem is transcribed into a finite-dimensional nonlinear programming problem via Hermite-Simpson method [2]. To solve this nonlinear programming problem we use a numerical solver -IPOPT [18].
3. Adaptive mesh refinement algorithm. The strategy proposed here is an improvement of the one reported in [14,15] and, again, it is based on block-structured adaptive refinement method which became popular within fluid mechanics since multigrid algorithms can be used for time and space domains.
In this approach, we are concerned in solving the problem, in a first step, on a coarse mesh and, according to a decision based on some refinement criteria, the mesh is refined locally or entirely, using different levels of refinement. The discretized problem is then solved on the new refined mesh using information from the coarser mesh solution of the previous step. According to this strategy, we do not need to remove nodes.
3.1. Single-level and multi-level mesh refinement. The adaptive mesh refinement process starts by discretising the time interval [t 0 , t f ] in a coarse mesh, π 0 , containing N 0 equidistant nodes. After being transcribed into a NLP problem, the OCP is solved in this coarse mesh to catch the main structure of the problem. Then, the mesh is progressively refined. According to some refinement criteria, the mesh is divided in K mesh intervals where (τ 0 , . . . , τ K ) coincide with nodes. These mesh intervals S k form a partition of the time interval, that is, while the mesh nodes have the property that According to the algorithm reported in [14], the subintervals S k that verify the refinement criteria are refined by adding a fixed number N of equidistant nodes between each two mesh points in such way that the refined mesh π j+1 contains the nodes of the prior one π j . This property is an important feature in block-structured schemes. The procedure is repeated until the stopping criterion is achieved. We also considered a more conservative approach by refining the neighbours of S k , i.e., S k−1 and S k+1 . The resulting mesh using this strategy is denoted further ahead by π R . Now, we present an improved version of this algorithm by introducing different levels of refinement in a single iteration. After selecting the intervals S k that verify the refinement criteria, they are divided into smaller subintervals according to the user-defined levels of refinement For example, the higher levels can be defined as multiple powers of 10 of the first levelε = 1, 10, 10 2 , . . . , 10 m ε 1 .
A subinterval S k,i is at the i th level of refinement if for i = 1, . . . , m, and it will be refined by adding N i of equidistant nodes between each two mesh points. This procedure adds more node points to the subintervals in higher levels of refinement, corresponding to higher errors, and it adds less node points to those in lower refinement levels.
By defining several levels of refinement, we get a multi-level time-mesh in a single iteration. The resulting mesh using this strategy is denoted further ahead by π M L .
3.2. Refinement and stopping criteria. In order to proceed with the mesh refinement strategy, we have to define some refinement criteria and a stopping criterion. We consider three refinement criteria: 1. the estimate of the relative error of the trajectory (primal variables) (ε x ), 2. the estimate of the relative error of the adjoint multipliers (dual variables) (ε q ), 3. a combination of both criteria. For the stopping criterion, we consider a threshold for the relative error of the trajectory.
In the first refinement criterion, the relative error estimate is, at each time, the difference between the obtained state trajectory and an higher order approximation of the solution of the dynamic differential equation. In this case, the solution is given by piecewise cubic polynomials using Hermite interpolation and, then, it is integrated using the Romberg quadrature. At each refinement iteration, the local error (ε x ) is computed and this information is taken into account when deciding if the refinement procedure should continue.
In the second case, we consider the multipliers q MP which are solution of the differential equation system (AS), (T) and (1) given by the Maximum Principle 2.1, and we also consider the multipliers q KKT obtained by applying the Kuhn-Tucker conditions to nonlinear optimization problem which results from the transcription of the optimal control. The relative error estimate is, at each time, the difference between the multipliers q KKT computed by the numerical solver and q MP computed by integrating numerically the adjoint equation given by the Maximum Principle. This criterion is chosen because these multipliers give sensitivity information. Furthermore, q MP are solution to a linear differential equation system, which can be easily solved in a faster way and with higher accuracy. At each refinement step, the local error of the multipliers (ε q ) is evaluated for each time instant t ε q (t) = ||q MP (t) − q KKT (t)|| and the procedure selects which time intervals should be further refined.
In the last case, we use both refinement criteria simultaneously and the procedure will continue until the stopping criterion is satisfied.
As stopping criterion, we consider the L ∞ norm of the relative error of the primal variables (ε x ). Even when using the relative error of the adjoint multipliers (ε q ) as refinement criteria, we still need to estimate the error on the trajectory since it is the stopping criterion. However, we do not need to compute ε x (t) for all t -as in the case of the refinement criteria -but just an estimate of ε which is much faster to obtain.
3.3. Warm start. Since the proposed procedure increases the number of nodes, more computational time would be expected. To decrease the CPU time, when going from a coarse mesh to a refined one progressively, the previous solution is used as a warm start for the next iteration. To create this warm start, the solution obtained in the coarse mesh is projected in the refined one using the cubic Hermite interpolation. This action proved to be vital in the decreasing of the overall computational time.
In particular, we notice that the number of iterations of the NLP solver remains within the same order of magnitude when we considerably increase the number of nodes.
3.4. Algorithm implementation. The overall procedure is described in Algorithm 1. To test the algorithm, the proposed procedure was implement in MATLAB R2008a combined with the Imperial College London Optimal Control Software -ICLOCS -version 0.1b [5]. ICLOCS is an optimal control interface and it uses the IPOPT -Interior Point OPTimizer -solver, which is an open-source software x ; estimate the error on the multipliers -ε (0) q ; while stopping criteria not met do select the mesh subintervals S j,i to be refined ; apply the discretization scheme according to the multi-level refinement criteria; transcribe the OCP ; create warm start; solve the NLP; estimate the discretization error -ε (j) x ; estimate the error on the multipliers -ε (j) q ; end Algorithm 1: Adaptive time-mesh refinement algorithm considering both refinement criteria 4. Application. To test and to validate the proposed algorithm, a problem involving a nonholonomic system [9,14] and another one regarding the SEIR model [4,10,12] are solved considering, in both cases, a pathwise state constraint. 4.1. Car-like system. A car-like system that moves in a plane generally has three degrees of freedom: translation along the two axes in the plane and rotation about the axis perpendicular to the plane. Nevertheless, these vehicles cannot move freely in all three degrees of motion due to their steering constraints. These kind of models are highly nonlinear and, for that reason, it is expected that a refined mesh having non equidistant spacing is more suitable. Such problems belong to the class of nonholonomic systems [9]. It is known that they require discontinuous controls. So, methods in which the parameterizations of the control are based on polynomial, such as pseudo-spectral methods, are less adequate.
In Fig. 2, the geometry of a car-like system is given. For a given time t, (x(t), y(t)) is the position of mid-point of the axle connecting the rear wheels, ψ(t) is the yaw angle, δ(t) is the steering angle and l is the wheelbase of the vehicle, i.e., the distance between its front and rear wheels. We consider also the curvature c(t) which relates to the steering angle δ and the minimum turning radius by tan(δ(t)) l and R min = 1 |c max | .
x y l ψ δ Figure 2. Car-like system geometry Let us consider t ∈ [0, 10], in seconds, x(t) = (x(t), y(t), ψ(t)) and u(t) = (u(t), c(t)). Aiming minimum energy, the car-like system problem (P 1 ) can be stated as: subject to • the dynamic constraintṡ a.e. t ∈ [0, 10] (4) where u(t) is the speed and c(t) is the curvature, • the input constraints • the end-point constraints where r 2 = 0.1 and x f = (x f , y f , ψ f ) = (10, 0, 0) is a user-defined target point, and • the pathwise state constraint where (x,ȳ) = (5, 1) and k = 10. The goal is to drive this car-like system with minimum energy from x 0 to some point near x f according to the terminal condition (6) while overcoming the state constraint (7).

4.2.
The SEIR model. The SEIR model is a compartmental model that describes the spreading of an infectious disease among a population (N ) by dividing it into four different compartments: susceptible (S), exposed but not yet infectious (E), infectious (I) and recovered (R) [4,10,12]. SEIR models can represent many human infectious diseases such as measles, pox, flu or dengue. According to [4], we can add to the dynamical system given in [12] an extra variable (W ), which stands for the number of vaccinated people and which is governed by the differential equatioṅ W (t) = u(t)S(t). Then, the ordinary differential equation governingṘ can be replaced with the one forṄ .

4.3.
Numerical results. Both problems, (P 1 ) and (P 2 ), are solved considering four meshes: a) π ML : obtained by the multi-level time-mesh refinement strategy; b) π R : obtained by the (single-level) time-mesh refinement strategy; c) π F : considering equidistant spacing with the smallest time step of π ML ; d) π S : considering equidistant spacing with the same number of nodes of π ML . initial susceptible population 1000 E 0 initial exposed population 100 I 0 initial infected population 50 R 0 initial recovered population 15 N 0 initial population 1165 W 0 initial vaccinated Population 0 For problem (P 1 ), we consider As it can be seen in Fig. 3, the car-like system successfully overcomes the obstacle and it stops when the terminal condition (6) is satisfied.  Fig. 4b, where the controls associated to (P 1 ) are shown, the constraint for the curvature c(·) becomes active at the start, in the middle and at the end of the trajectory.  The local errors of the trajectory for all meshes are shown in Fig. 5a using a logarithmic scale. The subintervals that need refinement are at the start, in the middle and at the end of the time interval, since the local errors are greater than the user-specified threshold, coinciding with the subintervals where the curvature is nonzero and the constraint for the curvature becomes active. We also notice that there are different subintervals belonging to different levels of refinement which indicates that the procedure to generate π ML is quite distinct from the one used to generate π R .
The adjoint multipliers are presented in Fig. 4c and the local error associated to them is shown in Fig. 5b. Comparing Fig. 5a with Fig. 5b, we see that the errors in the trajectory and in the adjoint multiplier have a similar structure along time, further validating the use of the adjoint multiplier in the refinement criteria.
The numerical results concerning the four meshes are shown in Table 2, which shows information about the number of nodes, the smallest time step, the number of iterations needed to solve the NLP problem, the objective functional, the maximum absolute local errors of the trajectory and the q multipliers, and the CPU times for solving the OCP problem and for computing the local error as well.  Table 2, the mesh π ML has only 32.4% of the nodes of π F , nevertheless both meshes have maximum absolute local errors of the same order of magnitude. Since the procedure to obtain π ML uses a warm start at each refinement iteration, the OCP can be solved three consecutive times and it is still much faster than to solve this problem with the mesh π F . In fact, computing the solution using π ML takes only 4% of the time needed to get a solution using π F , causing significant savings in memory and computational cost. The use of multi-level refinement algorithm in real time optimization problems, such as MPC, has benefits    Table 2, if the procedure is interrupted after 12 seconds, a solution with local error lower than 4.096 × 10 −4 is obtained.
The mesh π S has the same number of nodes of π ML but considering equidistant spacing. The analysis of the solution obtained using this mesh allows us to verify the importance of nodes collocation, i.e., having a mesh with non-equidistant spacing nodes produces a solution with higher accuracy than the one obtain using a mesh with equidistant spacing nodes. Moreover, the CPU time spent to compute solution using π ML is 31, 3% lower then the one spent to obtain a solution using π S , emphasizing the relevance of using meshes with non-equidistant spacing nodes.
When comparing both refined meshes, we notice that the process to compute the mesh π ML having several refinement levels in a single iteration took only 2 refinement iterations, producing a mesh that has 95% of the nodes of π R . The solution is obtained 6% faster when compared to the CPU time spent to compute the solution using π R . We could expect to be even faster but, as shown in the Table  2, the second refinement iteration takes 50 IPOPT iterations to converge to the optimal solution. This event occurs because the solution of the previous refinement iteration, which is used to create a warm start, has only about 10% of the nodes, having less information about the structure of the solution.
In terms of CPU time, we can see that it is much faster to compute ε q than ε x . In the all procedures, the use of ε q in the refinement criterion reduces the computational time, making the refinement algorithm faster.
With respect to IPOPT and according to Table 2, even if the number of nodes is increasing fast at each refinement step, the number of IPOPT iterations are of the same order of magnitude. This is another advantage of the mesh refinement strategy.  The optimal trajectory is shown in Fig. 6a and 6b and the control is presented in Fig. 6c. The local errors of the trajectory for all meshes are shown in Fig. 6d using the logarithmic scale. Regarding the latter figure, we see that the time domain needs to be entirely refined in the first step and there are different subintervals belonging to different levels of refinement.
The numerical results concerning the four meshes are shown in Table 3, which, as before, shows information about the number of nodes, the smallest time step,  Figure 6. Results for problem (P 2 ) the number of iterations needed to solve the NLP problem, the objective functional, the maximum absolute local error of the trajectory and the CPU times for solving the OCP problem and for computing the local error as well.
According to Table 3, the mesh π ML has 85% of the nodes of π F and computing the solution using π ML takes only 30% of the CPU time needed to get a solution using π F . Nevertheless, we get solution with the same accuracy.
The solution obtained using the mesh π S , which has the same number of nodes of π ML but with equidistant spacing, has less accuracy than the one computed on π ML . Moreover, computing the solution on π S took 3 times the CPU time spent to get the solution using π ML .
When comparing both refined meshes, we notice that the process to compute the mesh π ML having several refinement levels in a single iteration took only 2 refinement iterations, producing a mesh that has 95% of the nodes of π R . The use of the refinement levels brings a significant improvement in the overall computing time, since the solution is obtained 19% faster when compared to the time spent to compute the solution using π R .
With respect to CPU time, once again, we see that it is much faster to compute ε q than ε x . Also in this application, the use of ε q as refinement criterion reduces the computational time, making the refinement algorithm quicker.
In terms of IPOPT and according to Table 3, the number of IPOPT iterations are of the same order of magnitude at each refinement step, even if the number of nodes is increasing.

5.
Characterization of the solution using the necessary conditions of optimality. In this section we establish necessary conditions of optimality for the car-like system problem (P 1 ) considering the pathwise state constraint (7). Then, we verify that the solution obtained numerically satisfies the necessary conditions. Similar analysis to the second problem (P 2 ) involving the SEIR model is reported in [4].