INEXACT RESTORATION AND ADAPTIVE MESH REFINEMENT FOR OPTIMAL CONTROL

. A new adaptive mesh reﬁnement algorithm is proposed for solving Euler discretization of state-and control-constrained optimal control problems. Our approach is designed to reduce the computational eﬀort by applying the inexact restoration (IR) method, a numerical method for nonlinear programming problems, in an innovative way. The initial iterations of our algorithm start with a coarse mesh, which typically involves far fewer discretization points than the ﬁne mesh over which we aim to obtain a solution. The coarse mesh is then reﬁned adaptively, by using the suﬃcient conditions of convergence of the IR method. The resulting adaptive mesh reﬁnement algorithm is convergent to a ﬁne mesh solution, by virtue of convergence of the IR method. We illustrate the algorithm on a computationally challenging constrained optimal control problem involving a container crane. Numerical experiments demonstrate that signiﬁcant computational savings can be achieved by the new adaptive mesh reﬁnement algorithm over the ﬁxed-mesh algorithm. Conceivably owing to the small number of variables at start, the adaptive mesh reﬁnement algorithm appears to be more robust as well, i.e., it can ﬁnd solutions with a much wider range of initial guesses, compared to the ﬁxed-mesh algorithm.


1.
Introduction. Optimal control problems have many applications in science and engineering and, due to their importance, finding solutions of these problems accurately and efficiently is of considerable interest. Because of the difficulty in dealing with the continuous-time optimal control problem, much interest focuses on approximating these problems and applying a numerical method to obtain an approximate numerical solution. The continuous-time optimal control problem is often approximated by a finite dimensional optimization problem employing Euler, midpoint and more generally Runge-Kutta discretization schemes (see, for example, [2,7,8,9,12,15,16,18,21,28]). Euler discretization is the simplest of these schemes and is often the preferred approximation, because the problem can then be easily manipulated and the conditions of optimality easily formulated. The Euler discretized problem is numerically solved by using a finite-dimensional nonlinear optimization method. In this process, two issues need to be addressed: (i) efficiency of the finite-dimensional optimization method and (ii) convergence of the approximate for optimal control problems has also been studied in [17], where refinement of the mesh is based on a local error produced by applying Runge-Kutta discretization and the resolution of an integer linear programming problem. Yet another mesh refinement technique is introduced in [30] with the refinement based on a mesh density function, which generates a nonuniform mesh by suitably allocating the grid points over the entire time interval. In all of the aforementioned studies, mesh refinement is based on an estimate of some kind of local error.
The sufficient conditions of convergence of IR, which are used as a criterion for the refinement of mesh in the AMIR algorithm, plays a similar role to that of the local error described in the AMR techniques existing in the literature. To the best of our knowledge, using IR type of conditions for mesh refinement is new.
The paper is organized as follows. In Section 2, the optimal control problem and its Euler discretization is stated. In Section 3, the authors' earlier convergence result in [2] on the IR method for the fixed-mesh Euler discretization of constrained optimal control problems is cited. In Section 4, the AMIR algorithm is described in detail. In Section 5, numerical experiments involving a container crane are carried out in order to illustrate the functionality and robustness of the algorithm, as well as to demonstrate computational savings.
2. The optimal control problem and its approximation. Consider an optimal control problem subject to box constraints on the state and control variables as follows. (P) subject toẋ(t) = f (x(t), u(t)) for a.e. t ∈ [t 0 , t f ], x(t 0 ) = x 0 , ξ(x(t f )) = 0 , a i ≤ x i (t) ≤ a i , , for all j = 1, . . . , n, and all t ∈ [t 0 , t f ] , b k ≤ u k (t) ≤ b k , for all k = 1, . . . , m, and a.e. t ∈ [t 0 , t f ] , where the state variable x ∈ W 1,∞ (t 0 , t f ; R n ),ẋ = dx/dt, the control variable u ∈ L ∞ (t 0 , t f ; R m ). The functions f 0 : R n × R m → R, f : R n × R m → R n , ξ : R n → R n and φ : R n → R are twice continuously differentiable. The initial state is specified as x 0 . In this problem t 0 and t f are fixed and L ∞ (t 0 , t f : R n ) corresponds to the space of essentially bounded, measurable functions equipped with the essential supremum norm. Furthermore, W 1,∞ (t 0 , t f ; R n ) is the Sobolev space consisting of functions x : [t 0 , t f ] → R n whose first derivative lies in L ∞ . In order to approximate Problem (P) by the Euler scheme, we first subdivide the time horizon [t 0 , t f ] into N pieces. Subdivision points t i , i = 0, . . . , N, can be non-uniform with the mesh size ∆t i = (t i+1 − t i ). In this case, t N = t f . Define the partition π = {t 0 , t 1 , . . . , t N } , and the vectors x π = (x T 0 , x T 2 , . . . , x T N ) ∈ R n(N +1) , u π = (u T 0 , u T 2 , . . . , u T N −1 ) ∈ R mN , where x i and u i are the approximation of the state and control variables at time t i , respectively, so x i ≈ x(t i ) and u i ≈ u(t i ) for i = 0, . . . , N . Introduce the forward difference

NAHID BANIHASHEMI AND C. YALÇ IN KAYA
and approximate the integral The Euler discretization of Problem (P) is then simply given by We denote by x ij the jth component of x i , j = 1, . . . , n, and by u ik the kth component of u i , k = 1, . . . , m, i = 0, . . . , N . Problem (PE) is a finite dimensional optimization problem and by solving it we can find an approximate solution for Problem (P). In the following section we apply the IR method to this problem.
3. IR for optimal control. In this section, we will formulate the IR method for solving optimal control problems. For the sake of simplicity in appearance, we define the functions for the equality constraints as follows.
Problem (PIR) is a box-constrained optimization problem. Birgin and Martínez [5] present the IR method to solve this class of problems. Their method is local and iterative with two distinctive phases, namely feasibility and optimality . It starts from an initial guess which is close enough to a local solution and converges to a local solution of the problem under conditions which are checked during the iterations. The idea of their method can be translated into our setting of the Euler-discretized state and control constrained optimal control problems as follows.
The IR iterations start with a given point (x π , u π ) ∈ Ω and first find a "more feasible" point (y π , v π ) ∈ Ω. This process is referred to as the feasibility phase and is responsible for improving the feasibility of the current point (x π , u π ). The next phase is the optimality phase which can be described as follows.
The optimality phase is responsible to find a "more optimal " point (z π , w π ) ∈ Ω in the plane tangent to the constraints and passing through (y π , v π ). In order to formulate these two phases we define the Lagrangian for Problem (PIR) which is given, for all (x π , u π ) ∈ Ω, by where σ ∈ R n and λ π = (λ T 0 , . . . , λ T N ) ∈ R (N +1)n are the Lagrange multipliers. We where ∇ is the gradient operator with respect to (x π , u π ), and P is the projection operator on Ω with respect to the Euclidian norm · . Now we start describing the feasibility phase of the IR method. In order to improve the feasibility of the (x π , u π ) ∈ Ω, one should find a point (y π , v π ) ∈ Ω which satisfies h(y π , v π ) ≤ θ h(x π , u π ) , θ ∈ [0, 1). ( The new point (y π , v π ) is required to be not too far from (x π , u π ) within the set Ω, in some sense. The condition of "not being too far" is expressed by the following inequality.
where K 1 is a positive constant. An IR iteration continues with the optimality phase for finding a point (z π , w π ) ∈ Ω in the plane tangent to the constraints described by such that "optimality of the point (y π , v π )" is improved. The tangent plane given in (5) can be expressed equivalently by where we have used the short-hand notation A i = ∂f /∂x(y i , v i ) and B i = ∂f /∂u(y i , v i ). As in the case of the feasibility phase, the "more optimal" point (z π , w π ) is also required to be "not so far", in some sense, from the point (y π , v π ) within the set Ω. A condition to meet this requirement will be given later. This condition will involve not only the "more optimal" point (z π , w π ) but also the Lagrange multipliers, which are defined below. In order to find (z π , w π ), one minimizes the Lagrangian L(z π , w π , λ π , σ) constrained to the tangent plane (6)-(8) as well as to Ω. Now for all (z π , w π ) ∈ Ω, we define the Lagrangian function for the optimality phase by L(z π , w π , µ π , σ) = L(z π , w π , λ π , σ) where µ π = (µ 0 , . . . , µ N ) ∈ R n(N +1) and σ ∈ R n . The σ − σ and (µ i − λ i ) are the Lagrange multipliers corresponding to the linear constraints (6)- (8). The necessary condition of optimality for this subproblem can be written as G(z π , w π , σ, µ π , λ π , v π , y π , σ) := P (z π , w π ) − ∇ L(z π , w π , µ π , σ) − (z π , w π ) = 0, (9) which is expanded in [2]. Replace ∆t with ∆t i in that expansion.
To check the improvement of the optimality of (z π , w π , µ π , σ), the following IR conditions should be satisfied and are given in [16] and [2].
Theorem 3.1 and the convergence theory of Dontchev et al. [8] presented in [2] guarantees convergence of a local solution of Problem (PIR) obtained by the IR method to that of the original problem (P), as stated in the corollary below. Note that Dontchev et al. [8] use the indirect adjoining approach with continuous adjoint variables [13] for applying the minimum principle. This approach is used for the box-constrained optimal control problem by Banihashemi and Kaya [2] and the Lagrange multipliers they have used are denoted by ψ, ν 1 , ν 2 , λ, σ and κ. It is worth mentioning that, by applying the IR method to Problem (PIR), only the Lagrange multipliers of the direct adjoining approach [13] are obtained. Direct and indirect adjoining approaches differ in the way the corresponding Hamiltonian functions are defined. However, the Lagrange multipliers of the indirect adjoining approach can be obtained through the multipliers of the direct adjoining approach. Refer to [2] for the details about the multipliers of these approaches and their transformation to one another. [2]) hold and Conditions (3), (10)- (12) and v π −u π ≤ K c h(x π , u π ) , where K c is a positive constant, are satisfied. Then, as ∆t → 0,a local solution (x * π , u * π , σ * , λ * π ) of Problem (PIR) found by the IR method will converge to a local solution (x * , u * , σ * , λ * ) of the original problem (P), where λ * is the adjoint variable for the direct adjoining approach. One obtains the adjoint variable for the indirect adjoining approach as ψ * = λ * + ν 1 * − ν 2 * .

Corollary 1 (Corollary 4.1 in [2]). Suppose Assumptions (A1)-(A6) (as listed in
In the IR algorithm presented in [2], the subdivision points are considered equidistantly and the mesh size N is fixed throughout the whole IR iterations. We call this IR algorithm the fixed-mesh IR (FMIR) algorithm where the procedure is described in broad terms as follows.
(1) Take an initial iterate for the state and control variables and the multipliers.
Take the solution of the optimality phase as the current iterate and go to (2) above. Banihashemi and Kaya [2] show by means of numerical experiments that the FMIR algorithm using the optimization software Ipopt in each phase is more robust than using Ipopt on its own directly in finding a local solution of Problem (PIR). A more robust method can handle the optimization process starting from a wider range of initial guesses of the local solution. This ability is an important motivation for applying the IR method in solving box-constrained optimal control problems.
In terms of elapsed CPU-time for finding a local solution, the FMIR algorithm using Ipopt in general takes a longer time compared to Ipopt on its own according 528 NAHID BANIHASHEMI AND C. YALÇ IN KAYA to [2], although FMIR is more robust than Ipopt on its own. Decreasing the running time required by the FMIR algorithm is an important research question since this algorithm is known to enjoy local convergence, as well as robustness. It was mentioned earlier that an IR iteration is completed if the feasibility and optimality conditions (3), (4) and (10)-(12) are satisfied. As can be seen, these conditions are stated inexactly, which allows flexibility in implementation. In this paper, we improve on the FMIR algorithm by exploiting this flexibility in such a way that the CPU-time for finding a solution is decreased without compromising the accuracy of the solution. We develop an Adaptive Mesh Refinement (AMR) approach by employing the flexible features of the IR method, which we present in the next section.
4. Adaptive mesh refinement. In this section, we convey the basic idea of the AMR approach, intertwined with the IR method to solve Problem (PIR). Furthermore, the algorithm arisen from combining the AMR approach and IR, called the adaptive mesh IR (AMIR) algorithm, is described at the end of the section. Before explaining the AMR approach we need to define the so-called coarse and finest mesh, first.
The mesh for which we wish to obtain a solution is called the finest mesh. For this mesh, N f denotes the mesh size. The finest mesh is represented by the partition where N c is the size of the coarse-mesh. Usually, N c is much smaller than N f . The process of increasing the number of subdivisions in the discretization of the time horizon is called refinement and can be implemented locally in the mesh or over the entire mesh. Subsequent refinements eventually lead to the finest mesh.
By going from a coarse mesh to the finest mesh progressively, we aim to decrease the CPU-time elapsed in solving Problem (PIR) by the IR method. Generally, having a finer mesh (i.e., a larger N c ) for solving Problem (PIR) results in an improvement in the accuracy of the solution. However, the finer the mesh is, the more computational time is required. As mentioned in the introduction, AMR can reduce the computational time needed to find a local solution of the problem on the finest mesh. The idea (or concept) of AMR is implemented in various ways by the practitioners in [4,14,17,30] to solve optimal control problems. The common feature in these approaches is that the discretized problem is first solved on a coarse mesh and, according to a decision based on some local truncation error, the mesh is refined locally or entirely. The discretized problem is then solved on the new mesh. As described in the survey given in the introduction, these approaches differ in the manner in which the refinement is carried out. The way we define a refinement is different from these approaches: we employ some features of the IR method, as we explain below. 4.1. AMR using IR. The main purpose of this section is to present the AMIR algorithm to solve Problem (PIR) on the finest mesh in less computational time than that of the FMIR algorithm, previously introduced in [2].
Consider the initial iterations of the FMIR algorithm as broadly explained in the previous section. Both feasibility and optimality phases of these iterations are normally done on the finest mesh with size N f . We would hope that having a coarsemesh with N c << N f in the initial iterations would help us achieve computational savings. Therefore, we change the structure of the initial iterations of the FMIR algorithm in such a way that the subproblems of both phases in these iterations are accomplished over a coarse mesh. We call these iterations coarse-mesh iterations in our new algorithm AMIR.
In the AMIR algorithm, the most important feature is the mechanism with which we do the refinements. In this regard, we plan to employ the IR functions G and G and Conditions (3), (10) and (11) for a refinement. The coarse-mesh iterations will presumably result in a very good initial guess for the final AMIR iterations that are to be carried out with the finest mesh. It should be noted that the AMIR iterations with the finest mesh are exactly the same as those of the FMIR algorithm.

The Coarse Mesh
Before going into further details of the new algorithm, we need to define the optimization variables for Problem (PIR) over the coarse mesh, π c = {t c 0 , t c 1 , . . . , t c Nc }, as follows.
Here, we note the approximations . , N f . While the vector (x πc , u πc , λ πc , σ) is defined over the coarse-mesh, the vector (x π , u π , λ π , σ) is defined over the finest mesh. For ease of bookkeeping, we consider N f = 2 M and N c = 2 p , where M and p are positive integers, with p < M . We rewrite Problem (PIR) on a coarse mesh as follows. (PIR-C) subject to h(x πc , u πc ) = 0 , (x πc , u πc ) ∈ Ω .
As mentioned before, the initial iteration of the AMIR algorithm which is called the coarse-mesh iteration has a different structure from that of the FMIR algorithm. This iteration starts with the feasibility phase on a uniform mesh with N c number of subdivisions. In this phase, (y πc , v πc ), which is a more feasible point than (x πc , u πc ), will be found. Note that the elements of (x πc , u πc ) are chosen from (x π , u π ) in such a way that, for i = 0, . . . , N c , where k = iN f /N c . Obviously, the computational time needed for the feasibility phase with N = N c , where N c is small, is expected to be much less than the case with N = N f .

Interpolation From Coarse to Finest Mesh
The process in the feasibility phase continues with generating (y π , v π ) on the finest mesh by means of linear interpolation using the coarse mesh vector (y πc , v πc ). An approximation of the control at the terminal time, v c Nc , has not been computed in the feasibility phase yet but, its value is needed to generate (y π , v π ). For this terminal control value, we do an extrapolation: Due to a uniform distribution of the discretization points on the coarse mesh, , for i = 0, . . . , N c − 1, need to be generated for the finest mesh. Linear interpolation is applied to find approximate values of these intermediate points as follows.
where N g i is the number of the intermediate points that need to be generated for the ith subdivision. Note that if the distribution of points in the coarse mesh is not uniform, e.g., in a subsequent coarse-mesh iteration, N g i will be different for each

Feasibility Phase Refinement
Having obtained the interpolated (y π , v π ), Condition (3) is required to be satisfied to achieve improvement in the feasibility. If Condition (3) is not satisfied by (y π , v π ), the current coarse mesh is refined. Refinement in this stage is called the feasibility phase refinement. For refinement, each subdivision of the current mesh is halved to obtain a finer mesh with 2N c number of equidistant subdivisions, so N c is updated by 2N c and the feasibility phase is started with the new (x πc , u πc ), which is obtained as in (13). Note that the feasibility phase is carried out over a mesh with uniformly distributed points. Therefore, the feasibility phase refinement results in another uniform mesh. Consequently, a new (y π , v π ) is generated by interpolation and Condition (3) is checked again. The process of solving the subproblem of the feasibility phase and refining the mesh continues until Condition (3) is satisfied. The feasibility phase is followed by checking the stopping criteria, also given in [2], which, for a given small > 0, are stated as If (y π , v π , λ π , σ) satisfies (17), the AMIR process stops, otherwise, another refinement process is done before starting the optimality phase. Before explaining this additional refinement, for a given point (y r , v r ), r = 0, . . . , N f , we define its associated subdivision on the current coarse mesh. Recall that (y π , v π ) was generated by interpolation from (y πc , v πc ), so for a given (y r , v r ) of the finest mesh, there is some k ∈ {0, 1, . . . , N c − 1} such that t r ∈ (t c k , t c k+1 ). Here, [t c k , t c k+1 ] is called the associated subdivision of (y r , v r ) on the coarse mesh. In the case when k = 0 or N c , and t r = t c k , we assign two associated subdivisions of (y r , v r ) which are [t c k−1 , t c k ] and [t c k , t c k+1 ]. Now, we describe the post-feasibility refinement process as follows.

Post-Feasibility Refinement
The function G(y π , v π , λ π , σ), as the representative for the gradient of the Lagrangian L(y π , v π , λ π , σ), can be used to "distinguish" a point (y r , v r , λ r , σ) of the finest mesh which is rather "far" from the local solution (x * r , u * r , λ * r , σ * ) of Problem (PIR) at t r . Therefore, the function G serves as an "error estimator" for the refinement process. Note that σ is constant with respect to t so it is not discretized. However, we include it in the argument of G. The threshold, above which one classifies a point to be "far" from a solution, is defined by G(y π , v π , λ π , σ) ∞ /3. For refinement, first, we find a point (y r , v r , λ r , σ) in the finest mesh for which G r (y π , v π , λ π , σ) ∞ ≥ G(y π , v π , λ π , σ) ∞ /3, where G r (y π , v π , λ π , σ) is the rth component of the vector G(y π , v π , λ π , σ) given by (2). Note that G r (y π , v π , λ π , σ) is mainly associated with the approximate variables at the node t r . The associated subdivisions of the triplet (y r , v r , λ r ) on the current coarse mesh must be found for refinement. Suppose that t r = t c i ∈ π c . Then [t c i−1 , t c i ] and [t c i , t c i+1 ] are the two associated subdivisions, which should be refined. These subdivisions are further subdivided into N r subdivisions. By refining each associated subdivision, N r − 1 new time nodes are added to the set π c . The new nodes also have to be in the set π.
In the case when t r ∈ (t c i , t c i+1 ), we set r = i and only refine [t c i , t c i+1 ]. As N c is updated, the value of the state, control and co-state variables (y ic j , v ic j , λ ic j ) for j = 1, . . . , N r − 1, at the new time nodes in [t c i , t c i+1 ] need to be be calculated. We obtain them simply as where k = r+jL i and L i = N i g /N r . Recall that N i g is defined immediately after (16) for the interpolation formulas. The superscript ic stands for the new points added to the ith subdivision of the current coarse mesh. The post-feasibility refinements of [t c i−1 , t c i ] and [t c i , t c i+1 ] are illustrated in Figure 1 for N r = 4. The refinement process is continued until no points remained in the finest mesh at which Condition (18) holds. Note that distribution of the discretization points in the coarse mesh after this refinement will in general be non-uniform, because the coarse mesh is refined locally. Finally, we update (y πc , v πc , λ πc ) by including the new points (y ci j , v ci j , λ ci j ) and start the optimality phase.

Optimality Phase Refinement
In the optimality phase, a more optimal point (z πc , w πc , µ πc , σ c ) in the tangent plane passing through (y πc , v πc ) is found. As in the feasibility phase, (z π , w π , µ π , σ) for the finest mesh is generated by interpolating (z πc , w πc , µ πc , σ c ) using (15) and (16). Also w c Nc and µ c Nc are needed and estimated via the extrapolation in (14). Note that σ = σ c and µ π is given by the following equation.

Post-Optimality Refinement
The next step after satisfying the optimality conditions (10) and (11) is the preparation for the next feasibility phase, in which (z π , w π , µ π , σ) is the current iterate: we set x π = z π , u π = w π , σ = σ and λ π = µ π . The post-optimality refinement is carried out simply as follows: update N c in such a way that it is larger than the current value and N f (modN c ) = 0, i.e., N f is divisible by N c . Furthermore, choose N c + 1 number of equidistant subdivision points from π and (x π , u π , λ π ) to generate π c and (x πc , u πc , λ πc ) as given by (13). Note that the post-optimality refinement is designed to produce a uniform coarse mesh because we want to solve the next feasibility phase over a uniform coarse mesh. Completion of the post-optimality refinement is the end of a coarse-mesh iteration of the AMIR algorithm. We introduce a parameter N I to be the number of the coarse-mesh iterations for the AMIR algorithm. Even if N c is not equal to N f at the end of the N I th iterate, we switch to the finest mesh and take the approximate solution (z π , w π , µ π , σ) generated through linear interpolation as the initial guess to start the next AMIR iteration. The next iteration with the finest mesh is exactly the same as the one in the FMIR algorithm.
Local refinement is carried out only before and during the optimality phase (as post-feasibility and optimality refinements) because of the higher cost of solving the optimality phase subproblem. On the other hand, the feasibility and post-optimality refinements, which are related more to the the feasibility phase are done uniformly, in view of accelerating refinement and the fact that the feasibility phase subproblem is relatively less costly to solve. It is worth mentioning that the post-feasibility and post-optimality refinements are also designed to avoid accomplishing the feasibility and optimality refinements as much as possible since they are time consuming.

Adaptive-mesh IR algorithm.
Step 0: (Initialization) Set the inexact feasibility parameter, 0 < θ < 1, inexact optimality parameter, 0 < η < 1, the other algorithmic parameters, K 1 ,  Step 1: (Feasibility phase) Step 1.0: If k ≤ N I or N c = N f , then solve approximately (with accuracy f ) the minimization problem (Pfeas) to obtain (y πc ). If k > N I , then solve Problem (Pfeas) on the finest mesh to obtain (y If not able to obtain an approximate solution, stop and declare "Failure at the feasibility phase." Step 1.1: If N c < N f then obtain (v c Nc ) (k) by extrapolation (14) and generate the approximate solutions (y πc ) of Problem (Pfeas) by linear interpolation given by (15) and (16).
Step 1.2: If the approximate solution (y then go to Step 2. Step 1.3: (Feasibility phase Refinement) If k ≤ N I and Condition (21) is not satisfied then halve each subdivision of the current coarse mesh to generate a finer mesh with size 2N c and do the feasibility phase refinement as explained in Section 4.1. Then go to Step 1.0. If not able to find an approximate solution which satisfies (21), stop and declare "Failure at the feasibility phase." Step 2: (Stopping Criteria) If the approximate solution (y π , σ (k) ) ∞ ≤ , then stop the execution of the algorithm and declare the quadruplet (y π , σ (k) ) as an approximate solution of the problem.
Step 4: (Optimality Phase) Step 4.0: If k ≤ N I or N c = N f then solve approximately (with accuracy o ) the problem π and σ = σ (k) , check (20). Do the optimality phase refinement as explained in Section 4.1 and go to Step 4.0. If not able to find an approximate solution which satisfies Conditions (22) and (23), stop and declare "Failure at the optimality phase." Step 5: ) and go to Step 1.0. Otherwise, set N c = N f and go to Step 1.0.

Numerical experiments.
In this section, we solve a nontrivial optimal control problem from the literature by using the proposed AMIR algorithm we have described in the preceding section. For coding the algorithm we use AMPL [11], which is a modelling language for general finite-dimensional optimization problems. For solving the subproblems (Pfeas) and (Popt) in the new AMIR algorithm, we employ (interfaced with AMPL) the optimization software Ipopt [29], version 3.2.4s. We also solve the same optimal control problem by the FMIR algorithm in order to make comparisons between the computational performances of these two algorithms. Note that the FMIR algorithm is also implemented in APML, and Ipopt is used to solve the subproblems of the feasibility and optimality phases. The implementation of both algorithms is accomplished on a system with an Intel Core2 Quad CPU, at 2.83 GHz each, with a 3.25-GB RAM, using Windows XP OS Professional version 2002 service pack 3.
Container cranes are used to handle cargo in ship or rail road terminals. Optimal control of container cranes have been widely studied by many authors; see, for example, [1,22,25,27]. This test problem is also studied by Banihashemi and Kaya in [2]. They implement the FMIR algorithm for this problem in AMPL using Ipopt for solving the feasibility and optimality subproblems. Via numerical experiments, they observe that the FMIR algorithm using Ipopt in its subproblems is more robust than Ipopt on its own. Also, Assumptions (A1)-(A6) listed in [2] are checked and all satisfied for the container crane problem. Furthermore, the multipliers of indirect adjoining approach [13] and those obtained by the FMIR algorithm which are related to the direct adjoining approach [13] are obtained. We consider the same container crane problem here for our AMIR algorithm, to illustrate that further computational savings can be achieved.
The dynamical model, boundary conditions, control and state constraints of the container-crane are considered over the time horizon [0,9], as follows.
In order to solve this problem with the AMIR algorithm, we set the algorithmic parameters to be η = 0.9, θ = 0.9, K 1 = 10 3 , K 1 = 10 3 , K 3 = 0.3, K 2 = 10 6 and = 10 −6 . We generate initial guesses uniformly at random in given intervals for the state and control variable components, x si and u rj , such that and for s = 1, . . . , 6, r = 1, 2, i = 0, . . . , N f , j = 0, . . . , N f − 1. We set N c = 2 4 = 16 and N f = 2 10 = 1024. Because we want to compare the numerical results of the AMIR algorithm for this example with that of the FMIR algorithm, we explain the FMIR process and its solution along with the demonstration of the AMIR process.
In the AMIR algorithm, the first iteration for this application starts with the feasibility phase with 16 subdivisions to obtain (y πc , v πc ). Using the same algorithmic parameters and initial guesses in (28) for the FMIR algorithm, Problem (Pfeas) is solved with 1024 subdivisions to find (y π , v π ). For demonstrating the approximate solutions of these two algorithms in the feasibility phase of the first iteration, we show the graphs of y As can be seen from Figure 2, the iterates found by the AMIR algorithm (y c 4 ) (0) and (v c 2 ) (0) with just 16 subdivisions are already close to the optimal solution (shown by a thick curve), although they are obtained inexactly. It should be noted that the 16 AMIR iterates here are much closer to the solution than the 1024 FMIR iterates. This observation is a major motivation in designing the AMIR algorithm. We have run the codes for these two algorithms 50 times to see the differences between the first iterates of (y πc , v πc ) and (y π , v π ) with initial guesses generated as in (28). For visualization, we show the iterate components (y c 4i   The reason why we present the first iterates of both algorithms in the feasibility phase for 50 runs is that we want to show, starting with random initial guesses given by (28), that the approximate solutions obtained by both algorithms in this phase are still relatively close to the optimal solution. It also seems that all these iterates for a given variable have a similar pattern. However, as the mesh size increases, the FMIR algorithm starts having difficulties in approximating the solution in the feasibility phase of the first iterate. This observation shows that when the dimension of the feasibility subproblem increases in the FMIR algorithm, finding a approximate solution is hard and time consuming. In some cases, it is even impossible to find a local solution of Problem (Pfeas) using FMIR. This points to AMIR being more robust.
In order to report the CPU-time elapsed for finding (y πc , v πc ) and (y π , v π ) by the two algorithms in the first iteration, we run both codes for 50 times with an additional set of randomly generated initial guesses and different values for N f . We generate initial guesses uniformly in the following intervals:   for i = 0, . . . , N f , j = 0, . . . , N f − 1. These ranges are roughly consistent with the ranges of the minimum and maximum values of the control and state variables. They are also considered in [2] for numerical experiments. Note that these intervals are wider than those of (28). Therefore, solving Problem (Pfeas) starting with initial guesses from these intervals is more challenging. Table 1 shows the average CPU-time that Ipopt needs to find (y πc , v πc ) and (y π , v π ) in the first iteration of the AMIR and FMIR algorithm. Note that the CPU-times of the successful runs starting from the initial guesses (28) and (29) are averaged. The maximum number of iterations allowed in Ipopt is 1000 in this paper while this number is set to be 4000 in [2]. Decreasing this number is planned for comparing the performance of both AMIR and FMIR algorithms in a more challenging situation.
As can be seen from Table 1, the AMIR algorithm is significantly more successful in completing the feasibility phase right in the first iteration than the FMIR algorithm. Moreover, it requires much less time in completing the feasibility phase starting from initial guesses in (28). Note that (29) generates initial guesses which are quite far from a solution for the container crane example. The AMIR algorithm still seems to work very well starting with these bad initial guesses, while the FMIR algorithm failing in every single run. The average CPU-time reported for the feasibility phase of the first iteration of the AMIR algorithm in Table 1 includes the time needed for doing the linear interpolation (15) and (16) (15) and (16) for 50 runs satisfy the feasibility conditions (21).
So far we have made an attempt to illustrate the importance and effectiveness of the AMIR algorithm in the feasibility phase of the first iteration. The numerical experiments have shown that (y π ) obtained through the linear interpolation (15) and (16) has passed Condition (21) in all 50 runs. Therefore, they are now ready to be used in the next step, the post-feasibility refinement. As mentioned before, Condition (18) is checked for each (y π , σ (0) ) in the finest mesh for the post-feasibility refinement. We observe that with (y π ) at hand, this process requires refining the whole coarse mesh. Therefore, N c is updated to 64 (N r = 4 for the refinement) and the optimality phase is started with a uniform coarse mesh of 65 points.
For this example, N I = 2, so after two coarse-mesh iterations we switch the mesh from the coarse to the finest. For most of the test examples we studied, we took, N I as 2. We observed that, in this problem, two iterations with the same structure of those in the FMIR algorithm with the finest mesh were sufficient for the AMIR algorithm to be successful.
In order to compare the CPU-time elapsed for the whole process of the AMIR and the FMIR algorithms, CPU-times have been averaged over the successful runs of the 50 attempts. We have run both algorithms 50 times with different values of N f and initial guesses in (28) and (29). The maximum number of iterations allowed is set to be 1000 for Ipopt in each phase for both algorithms. This means that the run is deemed successful if Ipopt can solve each subproblem in both phases in less than 1000 iterations and the IR conditions are satisfied. Table 2 lists the numerical results with the initial guesses as in (28). As can be seen from this table, we obtain significant savings in the CPU-time elapsed in finding a local solution of the problem using the AMIR algorithm.  Table 3. Numerical results for the FMIR and AMIR algorithms with initial guesses generated uniformly as in (29).
We claim that the AMIR algorithm is in general more robust than the FMIR algorithm because it has a better success rate in terms of finding a local solution (especially when the initial guess is poor). Next we justify this claim numerically.
The importance and effectiveness of the AMIR algorithm is further justified by the numerical results presented in Table 3. This table presents the numerical results of 50 runs of both algorithms starting from the initial guesses as in (29) and shows the impressive performance of the AMIR algorithm in terms of robustness. We observe that in this case, the FMIR algorithm is not able to find the solution in the feasibility phase of the first iteration (see Table 1). Therefore, it fails right at the beginning. Our observation shows that the FMIR algorithm is more likely to find a solution if it passes through the first iteration, so choosing an AMR strategy is very effective because it reduces the dimensions of the feasibility and optimality subproblems significantly in the early iterations. It is common experience of the practitioners that a small-scale optimization problem is more likely to be successfully solved even with bad initial guesses, compared to a large-scale problem. 6. Conclusion. We have proposed the AMIR algorithm for solving Euler discretization of state-and control-constrained optimal control problems. The AMIR algorithm starts the initial iterations with a coarse mesh, which is refined progressively using the IR conditions, which guarantee convergence to a local minimum over a (much) finer mesh. Since the initial iterations are usually the most computationally costly ones, the new algorithm, taking advantage of having a coarse mesh in the initial iterations, seems to decrease the overall computational time considerably, compared to the fixed-mesh, FMIR, algorithm. By having a coarse mesh in the initial iterations, the dimension of the discretized problem is decreased significantly, too. Therefore, finding an approximate solution in these iterations becomes easier. We applied the AMIR algorithm to the computationally challenging container-crane problem, and illustrated both significant computational savings and robustness, compared to the FMIR algorithm.
Future work should include optimal control problems involving partial differential equations for which mesh refinement strategies are of paramount importance. It should be noted that, in these problems, not only time, but also space, is discretized, giving rise to much larger-scale problems, which makes mesh refinement techniques mandatory to use.