An Integrated Multicriterion hp-Adaptive Pseudospectral Method for Direct Optimal Control Problems Solving

Pseudospectral methods PMs for solving general optimal control problems OCPs attract an increasing amount of research and application in engineering. It is challenging to improve the convergence rate, the solution accuracy, and the applicability of PMs, especially for nonsmooth problems. Existing hp-adaptive PMs consider only one heuristic criterion, which cannot produce satisfactory performance inmany cases. In this paper, we propose a novel methodwhich integrates multicriterion to hp-adaptive PM, in order to further improve the performance. For this purpose, we first devise an OCP solving framework of hp-adaptive PM. We then design a multicriterion hpadaptive strategy which introduces prior knowledge, intermediate error and curvature as useful criterions for adaptive refinement. We last present an iterative procedure for solving general nonlinear OCPs. Results from two examples show that our method significantly outperforms competitors on the convergence rate and the solution accuracy. The method is practical and effective for direct solving of various OCPs in a broad range of engineering.


Introduction
There exist a large number of optimal control problems OCPs in the engineering area, such as aircraft trajectory optimization, robotic control, space rendezvous and docking, lunar landing trajectory design, wave maker, energy production, fatigue test, and so on 1-9 .Numerical solving methods for OCPs fall into two categories which are indirect methods and direct methods, respectively.Indirect methods require rigidly necessary conditions which often cannot be obtained in engineering practice.In direct methods, an OCP is transcribed into a finite-dimensional nonlinear programming problem NLP , in which the necessary conditions from the calculus of variations are not rigorous.Direct methods are becoming more and more popular for solving various OCPs with the well development of computer technology 1 .As a major class of direct methods, pseudospectral method PM attracts more and more attention from academic research and engineering application 2 .For example, SIAM news reported that PM was successfully applied to optimal control and motion planning of the international space station in 2007.The application helped NASA to save about one million dollars in three hours 3 .
The basic principle of PM is to approximate the state using a set of basis functions while a set of differential algebraic constraints are enforced at a finite number of collocation points.How to divide the element intervals and adjust the degree of basis function polynomial are two key components.Three kinds of approach are usually employed in order to achieve a specified error tolerance.The first one is named h-method, which uses many low-degree approximating subintervals.The terminology "h" denotes the method in which the degree of the elements is fixed and convergence is achieved by refining the mesh size h.The name of second one is p-method, which utilizes a few fixed numbers of intervals with a variable degree polynomial in each interval.The terminology "p" denotes the method fixes the mesh and achieves convergence by increasing the polynomial degree p of the elements.And the third approach is classified as hp-adaptive method, which allows for a variable number of approximating intervals with a variable degree approximation within each interval.The terminology "hp" denotes the method which allows for refinement in both the element size h, and the polynomial degree p simultaneously.The three approaches derive from the finite element method FEM in mechanics and fluid dynamics 10, 11 , which are h-FEM, p-FEM, and hp-adaptive FEM, respectively.
Gauss PM GPM 12 , Radau PM RPM 13 , and Lobatto PM LPM 14 are three most well-developed PMs.The collocation points of them are based on accurate quadrature rules.Their basis functions are typically Chebyshev or Lagrange polynomial.They all belong to p-methods and they achieve convergence by increasing the degree of basis function polynomial.p-methods have a simple structure and converge at an exponential rate for problems which their solutions are infinitely smooth 15, 16 .However, many OCPs in engineering have either nonsmooth solutions or nonsmooth problem formulations.The nonsmooth feature results in low convergence in p-methods and leads to an extremely large NLP in h-methods.
It is natural to combine the strongpoint of h-methods and p-methods to form a class of hp-adaptive methods to improve the applicability of PM.One typical hp-adaptive method was proposed by Darby et al. in 17 , their method is based on the relative curvature.If the ratio between the maximum curvature and the average curvature is sufficiently large, then the error is decreased by refining the mesh; the process is called h-refinement.Otherwise, the accuracy is improved by increasing the degree of approximating polynomial within an element; the process is named p-refinement.Another adaptive strategy decides h-refinement or p-refinement according to the convergence rate in each element 18 .The h-refinement is adopted only when exponential convergence is lost in an existing element.Based on the observation, single criterion is used by the aforementioned hp-adaptive methods.There are various widely used criterions, heuristic mechanisms, and adaptive strategies in the FEM area, which are of great help for us to design a more efficient hybrid adaptive algorithm.
hp-adaptive methods attract increasing amount of attentions in the research of FEM 11, 19-25 .Babuška discussed the mathematical properties of pand hp-adaptive methods 10 .Galvão et al. addressed a hp-adaptive least-squares spectral element method LS-SEM for solving hyperbolic partial differential equations 26 .Ben Dhia et al. proposed a new adaptive method based on an optimal control approach for adaptive modeling of an atomic-to-continuum coupling method constructed from the Arlequin framework 27 .In Mathematical Problems in Engineering 3 28 , Oden and Prudhomme expanded the adaptive control of modeling error to include ideas of statistical calibration, validation, and uncertainty quantification.Dorao and Jakobsen gave a use case of applying hp-adaptive LS-SEM for the population balance equation 29 .An hp-adaptive spectral solver for reactor modeling was given in 30 .Lu investigated the adaptive mixed FEMs for parabolic optimal control problems and introduced an adaptive algorithm based on posteriori error estimates to guide the mesh refinement 31 .Mitchell and McClain summarized several hp-adaptive strategies in FEM 11 .These strategies include prior knowledge of solution regularity, regularity estimation, type parameter, Texas three steps strategy, prediction error, nonlinear programming, and so on.Texas three steps include initialization, adaptive h-refinement, and adaptive p-refinement 32 .The adaptive strategy uses intermediate error as a criterion, which reduces the number of iterations efficiently.An adaptive strategy based on the local regularity of the solution was proposed in 33 , which includes a method for estimating the local regularity.The successful applications of various hp-adaptive strategies in the FEM research give us a good deal of enlightenment.The heuristic criterions and the adaptive strategies are valuable for developing an integrated multicriterion hp-adaptive PM for direct OCPs solving.
In this paper, we intend to integrate multicriterion to hp-adaptive PM, for the purpose of further enhancing the performance of computation and approximation.We first devise an OCP solving framework of hp-adaptive PM.We investigate the approximation method of PM and the assessment of approximation error.In the framework, we focus on the algorithm of adaptive strategy and refinement here, which is called hp-adaptive algorithm.We then design a multicriterion hp-adaptive strategy which introduces prior knowledge, intermediate error, and curvature as useful criterions.These heuristic criterions support both adaptive strategy and adaptive refinement.The criterions of intermediate error and curvature are complementary for adaptive strategy in various OCPs solving.After that strategy design, we present an iterative algorithmic procedure.Our method converges by increasing the polynomial degree in the smooth segments of a solution; meanwhile it adaptively refines the mesh for the nonsmooth segments in an efficient way.Our evaluation results show that the proposed method significantly outperforms hp-adaptive PM based on curvature and effectively improves the convergence rate of computation and the accuracy of solution.
The remainder of this paper is organized as follows.Section 2 defines multiple-interval nonlinear continuous OCP in Bolza form.Section 3 describes the OCP solving framework of hp-adaptive PM.Section 4 designs a novel multicriterion hp-adaptive strategy.Section 5 details our integrated multicriterion hp-adaptive PM.Section 6 provides two examples for illustrating our method.Finally, Section 7 concludes the paper.

Multiple-Interval Nonlinear Continuous Optimal Control Problem
Consider the following multiple-interval nonlinear continuous optimal control problem in Bolza form.Minimize the cost function τ 0 ≡ −1, and τ f ≡ 1.Furthermore, the inverse transformation is given as The nonlinear continuous-time optimal control problem on element t 0 , t f , which is defined by 2.1 -2.4 , can be expressed as a K element intervals problem.First, the cost function of 2.1 can be written as Next, the dynamics constraints, the boundary conditions, and the inequality path constraints are given, respectively, as and continuity constraints at element interfaces

Flow Chart of hp-Adaptive Pseudospectral Method
Pseudospectral method transcribes a continuous OCP into a nonlinear programming problem NLP .The flow chart of hp-adaptive PM is described in Figure 1.The major steps include pseudospectral transcription, solving NLP, error estimate, and adaptive strategy and refinement.In this paper, we choose Sparse Nonlinear OPTimizer SNOPT as the NLP solver.
Here we focus on the algorithm of adaptive strategy and refinement, which is called hpadaptive algorithm.hp-adaptive algorithm refines mesh and basis function adaptively, which is the key component of hp-adaptive PM. hp-adaptive algorithm involves hp-adaptive strategy, adaptive h-refinement, and adaptive p-refinement.During an iterative procedure, the algorithm assesses the approximation error of solution and analyses the curvature of state curve; and then hp-adaptive strategy decides which refinement h-refinement or p-refinement is more suitable.Adaptive h-refinement improves the accuracy of solution by refining the distribution of meshes and collocation points.And the goal of adaptive p-refinement is to achieve an exponential convergence rate by adjusting the degree of basis functions adaptively.In general, hp-adaptive algorithm integrates the strongpoint of both h-methods and p-methods to optimize the mesh and the polynomial degree of each element, which can simultaneously improve accuracy and computational efficiency.

Transcribe OCP Using Pseudospectral Method
In this paper, we choose the RPM 12, 34 as the foundation for transcribing a OCP into a NLP.For an element interval t k−1 , t k , t 0 ≡ t k−1 , t f ≡ t k , our method approximates the states and controls via the local interpolating polynomials with variable low degree.The state variables and the control variables in each element can be approximated by variable low-order Lagrange interpolating polynomials, which are expressed as where τ 1 , . . .τ N are collocation points, plus the start point τ 0 ≡ −1, N is the number of collocation points.L i τ , i 0, 1, . . ., N , is a basis of Lagrange polynomial with a variable degree.Then differential constraint of 2.8 can be transcribed into the algebraically form The initial state X 0 ≡ X τ 0 , and terminal state X f ≡ X τ f satisfies Gauss quadrature formula Differentiating X τ in 3.1 with respect to τ, we get The continuous objective function of 2.7 can be approximated via LGR quadrature where ω i is the weight coefficients of Gauss quadrature.And the boundary conditions of 2.9 is Combining the state and control variances, the trajectory is constrained by Through the above process, a continuous Bolza OCP is transcribed into a limited dimensional NLP.

Assessment of Approximation Error
In the iterative procedure, assessment of approximation error is not only a stopping criterion, but also useful information for adaptive refinement.If local solution on a mesh element has not satisfy the specified accuracy, then the distribution of collocation points needs to be modified, either by increasing the polynomial degree on the element, or by redividing the element.
The state and control variables are approximated by basis function polynomials.For the kth element interval t k−1 , t k in mesh, k 1, 2, . . ., K , we firstly transform the interval of variable t into the interval of τ ∈ τ 0 , τ f .Then the approximation polynomials of terminal state are given as where N is the number of collocation points in the kth element, i 0, 1, . . ., N .Hence the local error of the element is given as Furthermore, we can define the absolute error of the kth element interval as where h k t k − t k−1 is the length of the kth element interval.And the maximum relative error over all state variable components j of the differential equations is given by where x j,i .

3.13
ω j is the even value of the jth state variable on all N collocation points.η k is the maximum relative error of all state variables on the kth element interval.
Let η be a specified accuracy tolerance for the discretized differential algebraic constraints.If the maximum violation of the differential algebraic equations in the kth element η k is less than η, then the approximation of the kth element is considered to satisfy the accuracy tolerance.

Multicriterion hp-Adaptive Strategy Design
The local approximation error estimator indicates the element which should be refined, but it does not indicate whether the element is better to refine by an h-refinement or by a prefinement.It is no longer sufficient to guide the adaptive refinement.A method for making that determination is called hp-adaptive strategy.We learn the lesson of adaptive strategy in FEM.Based on the observation, we propose three classes of useful criterion for hp-adaptive strategy: prior knowledge, intermediate error, and curvature.
First, the priori knowledge of solution regularity is useful for hp-adaptive strategy.As one criterion, it guides a p-refinement on the smooth segments and an h-refinement near singularities or flections of state curve.For example, the linear elliptic partial differential equations with smooth coefficients only have point singularities near corners of the boundary and where boundary conditions change.
The priori knowledge-based criterion requires an estimator, which evaluates the regularity value of intermediate approximate solution.Here we define the regularity value as r.A method for estimating the local regularity of a solution on a given mesh was proposed in 33 , it indicated that where α is the parameter of a priori estimate for the error.α can be obtained from the previous computation of the error estimator.At the initial step of hp-adaptive strategy, the regularity value r k of approximate solution on the kth element is readily available by initialization.Let p 0 k and p k , respectively, denote the initial degree and current degree of basis polynomial in the k element.If p k 2 ≤ r k , then p-refinement is performed, otherwise h-refinement is adopted.Furthermore, let p 0 k r 0 k − 2 be an optimal choice for initial p-refinement.
Second, intermediate error can be adopted as a criterion for hp-adaptive strategy.Define an intermediate error where γ is a parameter generally ranging from 5 to 20, η is a specified error tolerance, for example let η 0.01 1% error .An appropriate value of γ should make tradeoff between the criterion of intermediate error and the criterion of curvature.From a lot of experimentations, we find that an adopted γ between 5 and 20 effects good complement of the multicriterions.
For the element k in mesh, if the estimated error η k > η I , it needs detailed hrefinement.The criterion based on intermediate error should lead to distribute the error more even on every subinterval.It has been shown that the error is equidistributed over all the subintervals for an optimal mesh using h-refinement 11 .Hence, an optimal h-adaptive mesh distribution should simultaneously reduce the global error and distribute the local error of each element equably.Assume κ k max and κ k be the maximum and mean value of κ k τ .Furthermore, the ratio of the maximum curvature to the mean curvature is defined

4.4
For the element k in mesh, if r k < r max where r max > 0 is a user-defined parameter , then p-refinement is adopted and a larger degree polynomial is used to obtain a better approximation.If r k ≥ r max , then h-refinement is adopted in this element.
Furthermore, we design an hp-adaptive strategy via combining prior knowledge, intermediate error and curvature as heuristic criterions.The flowchart of the integrated multicriterion hp-adaptive strategy is described in Figure 2. The priori knowledge of singularities or flections is adopted as heuristic criterion for initial h-refinement.The assessment of solution regularity is used to determine the initial degree of p-refinement.In the iterative procedure, the criterions of intermediate error and curvature are complementary for adaptive strategy in various OCPs solving.

Integrated Multicriterion hp-Adaptive Algorithm Procedure
In this section, we depict an integrated multicriterion hp-adaptive algorithm for a nonlinear continuous OCP with K intervals.The algorithm mainly includes six procedures.
Step 1 initialization .Firstly, the error tolerance η is specified, for example η 0.01 1% error .Then the intermediate error η I is specified, for example, η I 0.1.And r max is specified, for example r max 8.
Consider the initial h-refinement for each state equation, if the equation is a piecewise function, the initial h-refinement is located at the switch points of piecewise function; if there is a priori knowledge of singularities or flections of state curve, it guides the meshes division also.If there is no information about switch points or singularities, a coarse initial h-refinement is made with uniform mesh element size h 0 .The discrete approximation in each element is defined by the set of Legendre-Gauss-Radau points.And for the initial prefinement, the polynomial-order of kth element is determined via the assessment of local solution regularity, it is given as Step 2 solve the NLP using the current mesh .The problem is transcribed into the NLP on the current mesh, and then solving the NLP using sequential quadratic programming SQP .A posteriori error estimator is used to estimate the approximation errors for the intermediate solution on each element, which is described in Section 3.3.The accuracy of the approximation in each element is assessed by calculating the differential-algebraic constraints on the collocation points.
Step 3. If η k ≤ η, then continue proceed to next k 1 element .η k is the approximation error of the kth element in mesh.
Step 4 adaptive h-refinement .If either η k > η I , or r k ≥ r max , then h-refine the kth element.h-refinement keeps p p 0 fixed and adaptively constructs elements with varying size.It maintains refinement until the solution achieves the intermediate error η I and curvature ratio r k < r max .The new number of element intervals is computed by formula where η k / η is the ratio between current estimative maximum error and the specified error tolerance, η k /η I is the ratio between current estimative maximum error and the intermediate error, and • is the operator that rounds to the next integer.c h is an integer constant as a parameter for adjusting iterative speed, for example, c h 2. If there is some prior knowledge about state curve, c h can be adjusted for flections.Formula 5.1 determines a dissatisfactory element interval should be subdivided into how many subelements by h-refinement.The formula is based on two ratios.The first part depicts the approximation degree from the current error to the specified error tolerance; and the second part indicates the approximation degree from the current error to the intermediate error.They are calculated by the log operation, because h-refinement is hoped to attain an exponential convergence for the nonsmooth segment of function.The choices of logarithm are based on the quantity relation of the ratios.c h controls the growth in the number of element intervals.If c h is sufficiently large, the algorithm will use less iteration to converge to an acceptable solution, but the number of collocation points may increase quickly between iterations.If c h is small, the mesh will increase slowly, but the algorithm may require many iterations.An appropriate value of c h should make tradeoff between h-refinement and prefinement.
The locations of the new required-elements are determined using the integral of a curvature density function in a manner similar to that given in 35 .Specifically, let ρ τ be the density function given by where c is a constant to satisfy 1 −1 ρ ζ dζ 1.

5.3
The density function expressed by 5.2 is a probability density function pdf with fractional power-laws.Here we adopt 1/3 power in the function.In fact, from the point of view of power-laws in stochastic processes and/or fractal time series, 1/3 power is a typical heavytailed case.By heavy-tailed we mean that the density ρ τ decays slowly.And heavy-tailed pdfs imply that τ is in wild randomness due to infinite or very large variance 36 .Let F τ be the cumulative distribution function given by The n k new mesh points are chosen by Finally, it is noted that if n k 1, then no subintervals are created.Therefore, the minimum value for n k is 2.
Step 5 adaptive p-refinement .p-refinement keeps mesh and the size of elements fixed.It adjusts the degree of polynomials on each element until the final error is satisfied.The degree of polynomial p k on the kth element is determined by where p k is the initial degree of polynomial on this element.The growth in the polynomial degree is related to the log of the error ratio; because of a smooth function always exhibits exponential convergence in pseudospectral method.
Step 6 return to Step 2 .Our hp-adaptive algorithm integrates prior knowledge, intermediate error, and curvature for adaptive strategy and refinement, which is fit for smooth or nonsmooth problem.The illustration and comparison of our algorithm will be given by the examples in next section.
Remark 5.1.The accuracy and efficiency of the PMs for solving numerical optimal control problems have motivated the forefront of research.For continuous and smooth problems, pseudospectral knots and Gaussian quadrature rules are used to generate a natural spectral mesh that is dense near the points of interest.For nonsmooth problems, how to locate switches, kinks, corners, and other discontinuities is a challenge for solving practical OCPs.
In our hp-adaptive algorithm procedure, adaptive h-refinement is devised to locate the knots where the solution is subject to sudden changes, as in points of discontinuity.In fact, from the point of view of stochastic processes, the curve of solution is higher order discontinuities in the control time history.According to its statistic properties, it is fractal time series.In general, fractal time series have a heavy-tailed probability distribution function, a slowly decayed autocorrelation function ACF , and a power spectrum function PSD of 1/f type 37 .In our approach, we design the density function expressed by 5.2 .We adopt the 1/3 power, the density decays slowly, accordingly, heavy-tailed density.In our research of the practical OCP in engineering, such as aerospace engineering or mechanical engineering, the nonsmoothness and discontinuities are attributed to the differential dynamic and the constraint of system.The output or response of a differential system can be considered as a fractal time series 37, 38 .The density functions for mesh refinement in numerical optimal control are an instance of power-law-type pdf.The key problem of h-refinement is to find an appropriate density function.We mean that an appropriate choice of density function may help increase the accuracy of the solution and improve numerical robustness.Hence the theory of fractal

Numerical Examples
In this section, we present two numerical examples to illustrate the convergence rate and the applicability of the proposed method.The first example is a hypersensitive OCP, which solution has dramatic change on the "take-off" and "landing" segments.And the second example is a reusable launch vehicle re-entry trajectory optimization.The experimentations compare three algorithms: GPM and the boundary conditions where t f is fixed.This problem has a "take-off", "cruise", and "landing" structure where all of the interesting behaviors occur in the "take-off" and "landing" segments.In this example, t f 1000.
The specification of initial grid is listed in Table 1.The efficiency and accuracy of various methods are listed in Table 2. Figure 3    on the final mesh for hp 2 -3 method, at the accuracy tolerance grade η 1e − 3. The rings denote the location of collocation points.From analysis of these records, some results are illustrated.
First, hp-adaptive methods hp 1 and hp 2 outperform p-method in computational efficiency and approximation error.The results are better at all three accuracy tolerance grades.p-method even cannot achieve 1e − 3 accuracy tolerance forever.Because p-method does not refine mesh, which results in the transcribed NLP is dense for this kind of nonsmooth problem.It incurs dramatic increase of computation time but low accuracy of approximation.
Second, the proposed hp 2 method is better than hp 1 method in computational efficiency.At all three accuracy tolerance grades, the count of collocation points for hp 2 -3, hp 2 -5 is less than the count of hp 1 -3, hp 1 -5, respectively.The use of prior knowledge can get an accurate initial mesh.The criterion of intermediate error optimizes the distribution of collocation points.
Third, the count of collocation points for hp-3 is less than the count of hp-5, but hp-3 method needs more iterations.From the comparison of experimental records, hp 2 -3 is the best method for solving this problem, which attains the fastest computational efficiency and the least approximation error.
Fourth, as shown in Figure 3, the state and control curves rapid change on the "takeoff" and "landing" segments, where the relatively dense collocation points are distributed via the hp-adaptive algorithm.Meanwhile, the distribution of collocation points on the smooth segments is sparse.Compared with p-method and hp 1 method, our hp 2 method improves computational efficiency and reduces approximation error for this kind of nonsmooth OCP solving.
Example 6.2.Reusable launch vehicle re-entry trajectory optimization is an OCP of maximizing the cross range during the atmospheric entry 1 .Minimize the cost function where r is the geocentric radius, θ is the longitude, φ is the latitude, v is the speed, γ is the flight path angle, ψ is the azimuth angle, and σ is the bank angle.The radius of the earth R e 6371203.92 m.
The specification of initial grid is listed in Table 3.The efficiency and accuracy of various methods are listed in Table 4. Figure 4 shows the state and control curve of solution on the final mesh for hp 2 -6 method, at the accuracy tolerance grade η 1e − 4. The state and control curves are relatively smooth for this example.From analysis of these records, some results are illustrated.First, hp-adaptive methods hp 1 and hp 2 are similar or better than p-method in this example.p-method can achieve 1e − 5 accuracy tolerance and has the fewest number of collocation points, but it needs much iteration.The computational time of p-method is longer than hp-adaptive methods.The low computational efficiency of p-method is due to the fact that the transcribed NLP is dense.
Second, our hp 2 method outperforms hp 1 method in computational efficiency.The count of collocation points for hp 2 is less than the count of hp 1 at three accuracy tolerance grades.The comparison illustrates that the distribution of collocation points via hp 2 is better than the distribution via hp 1 .
Third, the count of collocation points for hp-6 is less than the count of hp-4.Based on the comparison of experimental records, hp 2 -6 is the best method for solving this problem, which has the fastest computational efficiency and the least approximation error.In general, our method is feasible and improves computational efficiency for this kind of smooth OCP solving.
Remark 6.3.The results of examples highlight several points of p-methods and hp-adaptive methods.p-methods are fit for solving problems which their solutions are smooth.hpadaptive methods perform similarly or better to p-methods for smooth problems.For nonsmooth problems, p-methods cannot accurately capture the discontinuities and rapid  changes in the trajectory.They only achieve low-accuracy even with much iteration.When a high-accuracy tolerance is required to satisfy, hp-adaptive methods are the likeliest candidates.They transcribe the OCPs into the NLPs with computational sparsity, which leads to the reduction of execution time.Compared with hp-adaptive PM based on curvature, our method is more efficient in the convergence rate; meanwhile the solutions achieve a comparable or better accuracy.The prior knowledge and the criterion of intermediate error leverage the enhancement of performance.The prior knowledge guides an accurate initial mesh, and the criterion of intermediate error optimizes the distribution of collocation points.

Conclusions
PMs become increasingly popular for direct OCPs solving in the engineering area.It is challenging to improve the convergence rate, the solution accuracy, and the applicability of algorithms, especially for nonsmooth problems.In this paper, we proposed a novel method which integrates multicriterion to hp-adaptive PM, in order to further improve the performance of computation and approximation.For this purpose, we first devised an OCP solving framework of hp-adaptive PM.We then designed a multicriterion hp-adaptive strategy which introduces prior knowledge, intermediate error, and curvature as useful criterions for adaptive refinement.The criterions of intermediate error and curvature were complementary for adaptive strategy in various OCPs solving.We last proposed an iterative procedure for solving general nonlinear OCPs.Our method converged by increasing the polynomial degree in the smooth segments of a solution; meanwhile it adaptively refined the mesh for the nonsmooth segments in an efficient way.Results from examples showed that our method significantly outperforms hp-adaptive PM based on curvature and effectively improves the convergence rate of computation and the accuracy of solution.Meanwhile, it enhances the applicability of hp-adaptive PMs for solving various OCPs.
For the future work, we will evaluate the performance of the proposed method using more OCPs in the engineering area.We will also try to exploit other valuable criterions for hp-adaptive method, which may further improve the convergence rate and the applicability of PMs.In addition, the choices of initial grid of PMs and the values of parameters γ, r max , and c h will be further researched for various problems.Different problem types will need dissimilar parameters for fast solving.Another topic is applying fractal time series and power-laws to guide the choice of density function for various problems.This point itself is an issue worth discussing in our future work.the functions are approximated by Lagrange interpolating polynomials of order N in which interpolation occurs at Gaussian quadrature points.There are three different families of Gauss quadrature points known as Legendre-Gauss, Legendre-Gauss-Radau, and Legendre-Gauss-Lobatto. Gauss PM GPM 12 , Radau PM RPM 13 , and Lobatto PM LPM 14 are the three most well-developed PMs.The collocation points of them are based on accurate quadrature rules.Their basis functions are typically Chebyshev or Lagrange polynomial.
Here we describe the procedure of GPM as an illustrative method.The GPM approximates the state using a basis of global interpolating Lagrange polynomials.These global polynomials are based on a set of discrete Legendre-Gauss LG points across the interval.
The standard interval considered here is denoted as τ ∈ −1, 1 .By using a linear transformation, the actual time t can be expressed as a function of τ via where t 0 and t f stand for the initial and final time, respectively.The state is approximated using a basis of N 1 Lagrange interpolating polynomials, L i τ i 0, . . ., N , as follows: where where τ i , i 0, . . ., N are the LG points belonging to the set.An approximation to the control, U τ , is formed with a basis of N Lagrange interpolating polynomials, L i τ i 1, . . ., N , as follows: A.4 The differential dynamic constraints should be posed to ensure that the entire discrete state satisfies the dynamic equations.The left-hand side of the dynamic equations is approximated by differentiating the state approximation of A.1 at the LG points as follows: A.5 The differentiation matrix, D ∈ R N× N 1 , can be determined offline from dynamic equations: A.6 The N collocation equations require A.5 to be equal to the right-hand side of the dynamic equations at the collocation points: As it is shown shortly, collocating strictly on the interior of the interval leads to a unique mathematical equivalence used to approximate the costate.
Combining the state and control variances, the constraints are discretized as follows: Lastly, the integral term in the cost functional can be approximated with a Gauss quadrature as follows: ω k g X k , U k , τ k ; t 0 , t f , A.9 where ω k is the weight coefficient of Gauss quadrature.Through the above process, a continuous OCP is transcribed into a limited dimensional NLP.The formulation is as follows: min F y s.t.
Initial h-refinement: meshes division based on the priori knowledge of singularities of state curve Adaptive h-refinementAdaptive p-refinementInitial p-refinement: for each mesh interval kSolve the NLP and assess the approximation error

b
shows the state and control curve of solution Mathematical Problems in Engineering Control u versus time t s

Figure 3 :
Figure 3: State and control curve on the final of Example 6.1 for hp 2 -3 and η 1e − 3.

bef
Longitude θ deg versus time t s Flight path angle γ deg versus time t s Azimuth angle ψ deg versus time t s

Figure 4 :
Figure 4: State and control curve on the final mesh of Example 6.2 for hp 2 -6 and η 1e − 4.
∈ R m is the state, u t ∈ R n is the control, C x t , u t , t ∈ R p is the control and state constraints, t is time, t 0 is the start time, and t f is the terminal time.Next, consider dividing the nonlinear continuous optimal control problem defined above into K element intervals.Let t 0 < t 1 < t 2 < • • • < t K , where t K t f is the terminal time.An elements k is defined to begin at the mesh point t k−1 and end at the mesh point t k .The time domain in each element k, t ∈ t k−1 , t k is transformed to τ ∈ −1, 1 by the formula time series and power-laws would guide the choice.Li gives a tutorial review about fractal time series 37 .The locally weak stationarity of modified multifractional Gaussian noise is discussed in 39 , which presents a set of stationary ranges.One-dimensional random function with long memory is introduced to model the sea level fluctuation 40 .Topics in power-laws are paid attention too, see for example, the contribution of approximating ideal filters by systems of fractional order in 41 , the work of addressing power laws in cyber-physical networking systems in 36 , the contribution of simplicial approach to fractal structures by Cattani et al. in 42 , and the work ofLi et al. in 38 .
12, hp-adaptive PM based on curvature 17 , and our integrated multicriterion hp-adaptive PM.For conciseness, the three algorithms are denoted by p, hp 1 , hp 2 , respectively.All CPU time recorded in table is whole time of iteration, which includes the time required to solve the NLP and the time required to refine the meshes.The program of these algorithms is expanded from open-source software GPOPS 43 .The transcribed NLP is solved by SNOPT 44, 45 , using a 2.4-GHz Core 2 Duo, 2G DDR RAM personal computer with MATLAB R2009b.For all examples, the following values are chosen for the parameters described in Section 4: r max 8 and γ 10.In hp 1 and hp 2 , the maximum number of collocation points on each mesh is 20, and the upper bound of iteration times is 20.

Table 1 :
Initial grid used in Example 6.1.

Table 2 :
Summary of results for Example 6.1 using various collocation strategies and accuracy tolerances.

Table 3 :
Initial grid used in Example 6.2.

Table 4 :
Summary of results for Example 6.2 using various collocation strategies and accuracy tolerances.