Research Article  Open Access
Hongfu Liu, Shaofei Chen, Lincheng Shen, Jing Chen, "An Integrated Multicriterion hpAdaptive Pseudospectral Method for Direct Optimal Control Problems Solving", Mathematical Problems in Engineering, vol. 2012, Article ID 760890, 22 pages, 2012. https://doi.org/10.1155/2012/760890
An Integrated Multicriterion hpAdaptive Pseudospectral Method for Direct Optimal Control Problems Solving
Abstract
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 adaptive PMs consider only one heuristic criterion, which cannot produce satisfactory performance in many cases. In this paper, we propose a novel method which integrates multicriterion to adaptive PM, in order to further improve the performance. For this purpose, we first devise an OCP solving framework of adaptive PM. We then design a multicriterion adaptive 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.
1. 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 finitedimensional 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 method, which uses many lowdegree approximating subintervals. The terminology “” denotes the method in which the degree of the elements is fixed and convergence is achieved by refining the mesh size . The name of second one is method, which utilizes a few fixed numbers of intervals with a variable degree polynomial in each interval. The terminology “” denotes the method fixes the mesh and achieves convergence by increasing the polynomial degree of the elements. And the third approach is classified as adaptive method, which allows for a variable number of approximating intervals with a variable degree approximation within each interval. The terminology “” denotes the method which allows for refinement in both the element size , and the polynomial degree simultaneously. The three approaches derive from the finite element method (FEM) in mechanics and fluid dynamics [10, 11], which are FEM, FEM, and adaptive FEM, respectively.
Gauss PM (GPM) [12], Radau PM (RPM) [13], and Lobatto PM (LPM) [14] are three most welldeveloped 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 methods and they achieve convergence by increasing the degree of basis function polynomial. 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 methods and leads to an extremely large NLP in methods.
It is natural to combine the strongpoint of methods and methods to form a class of adaptive methods to improve the applicability of PM. One typical 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 refinement. Otherwise, the accuracy is improved by increasing the degree of approximating polynomial within an element; the process is named refinement. Another adaptive strategy decides refinement or refinement according to the convergence rate in each element [18]. The refinement is adopted only when exponential convergence is lost in an existing element. Based on the observation, single criterion is used by the aforementioned 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.
adaptive methods attract increasing amount of attentions in the research of FEM [11, 19–25]. Babuška discussed the mathematical properties of  and adaptive methods [10]. Galvão et al. addressed a adaptive leastsquares spectral element method (LSSEM) 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 atomictocontinuum coupling method constructed from the Arlequin framework [27]. In [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 adaptive LSSEM for the population balance equation [29]. An 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 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 refinement, and adaptive 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 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 adaptive PM for direct OCPs solving.
In this paper, we intend to integrate multicriterion to adaptive PM, for the purpose of further enhancing the performance of computation and approximation. We first devise an OCP solving framework of 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 adaptive algorithm. We then design a multicriterion 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 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 multipleinterval nonlinear continuous OCP in Bolza form. Section 3 describes the OCP solving framework of adaptive PM. Section 4 designs a novel multicriterion adaptive strategy. Section 5 details our integrated multicriterion adaptive PM. Section 6 provides two examples for illustrating our method. Finally, Section 7 concludes the paper.
2. MultipleInterval Nonlinear Continuous Optimal Control Problem
Consider the following multipleinterval nonlinear continuous optimal control problem in Bolza form. Minimize the cost function subject to the dynamic constraints the inequality path constraints and the boundary conditions where is the state, is the control, is the control and state constraints, is time, is the start time, and is the terminal time.
Next, consider dividing the nonlinear continuous optimal control problem defined above into element intervals. Let , where is the terminal time. An elements is defined to begin at the mesh point and end at the mesh point . The time domain in each element , is transformed to by the formula , and . Furthermore, the inverse transformation is given as
The nonlinear continuoustime optimal control problem on element , which is defined by (2.1)–(2.4), can be expressed as a 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
3. Framework of Direct Optimal Control Problems Solving Using Adaptive Pseudospectral Method
3.1. Flow Chart of Adaptive Pseudospectral Method
Pseudospectral method transcribes a continuous OCP into a nonlinear programming problem (NLP). The flow chart of 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 adaptive algorithm.
hpadaptive algorithm refines mesh and basis function adaptively, which is the key component of adaptive PM. adaptive algorithm involves adaptive strategy, adaptive refinement, and adaptive refinement. During an iterative procedure, the algorithm assesses the approximation error of solution and analyses the curvature of state curve; and then adaptive strategy decides which refinement (refinement or refinement) is more suitable. Adaptive refinement improves the accuracy of solution by refining the distribution of meshes and collocation points. And the goal of adaptive refinement is to achieve an exponential convergence rate by adjusting the degree of basis functions adaptively. In general, adaptive algorithm integrates the strongpoint of both methods and methods to optimize the mesh and the polynomial degree of each element, which can simultaneously improve accuracy and computational efficiency.
3.2. 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 , 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 loworder Lagrange interpolating polynomials, which are expressed as where are collocation points, plus the start point , is the number of collocation points. , is a basis of Lagrange polynomial with a variable degree. Then differential constraint of (2.8) can be transcribed into the algebraically form where . The initial state , and terminal state satisfies Gauss quadrature formula Differentiating in (3.1) with respect to , we get The continuous objective function of (2.7) can be approximated via LGR quadrature where 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.
3.3. 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 th element interval in mesh, , we firstly transform the interval of variable into the interval of . Then the approximation polynomials of terminal state are given as where is the number of collocation points in the th element, . Hence the local error of the element is given as Furthermore, we can define the absolute error of the th element interval as where is the length of the th element interval. And the maximum relative error over all state variable components of the differential equations is given by where is the even value of the th state variable on all collocation points. is the maximum relative error of all state variables on the th 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 th element is less than , then the approximation of the th element is considered to satisfy the accuracy tolerance.
4. Multicriterion 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 refinement or by a refinement. It is no longer sufficient to guide the adaptive refinement. A method for making that determination is called adaptive strategy. We learn the lesson of adaptive strategy in FEM. Based on the observation, we propose three classes of useful criterion for adaptive strategy: prior knowledge, intermediate error, and curvature.
First, the priori knowledge of solution regularity is useful for adaptive strategy. As one criterion, it guides a refinement on the smooth segments and an 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 knowledgebased criterion requires an estimator, which evaluates the regularity value of intermediate approximate solution. Here we define the regularity value as . 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 adaptive strategy, the regularity value of approximate solution on the th element is readily available by initialization. Let and , respectively, denote the initial degree and current degree of basis polynomial in the element. If , then refinement is performed, otherwiserefinement is adopted. Furthermore, let be an optimal choice for initial refinement.
Second, intermediate error can be adopted as a criterion for adaptive strategy. Define an intermediate error where is a parameter generally ranging from 5 to 20, is a specified error tolerance, for example let (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 in mesh, if the estimated error , it needs detailed refinement. 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 refinement [11]. Hence, an optimal adaptive mesh distribution should simultaneously reduce the global error and distribute the local error of each element equably.
Third, curvature is another criterion for adaptive strategy. The curvature of state curve describes the regularity of intermediate solution. The curvature of the th component of the state in element is given by where is the th component of the state approximation in the th element [17]. Assume and be the maximum and mean value of . Furthermore, the ratio of the maximum curvature to the mean curvature is defined
For the element in mesh, if (where is a userdefined parameter), then refinement is adopted and a larger degree polynomial is used to obtain a better approximation. If , then refinement is adopted in this element.
Furthermore, we design an adaptive strategy via combining prior knowledge, intermediate error and curvature as heuristic criterions. The flowchart of the integrated multicriterion adaptive strategy is described in Figure 2. The priori knowledge of singularities or flections is adopted as heuristic criterion for initial refinement. The assessment of solution regularity is used to determine the initial degree of refinement. In the iterative procedure, the criterions of intermediate error and curvature are complementary for adaptive strategy in various OCPs solving.
5. Integrated Multicriterion Adaptive Algorithm Procedure
In this section, we depict an integrated multicriterion adaptive algorithm for a nonlinear continuous OCP with intervals. The algorithm mainly includes six procedures.
Step 1 (initialization). Firstly, the error tolerance is specified, for example (1% error). Then the intermediate error is specified, for example, . And is specified, for example .
Consider the initial refinement for each state equation, if the equation is a piecewise function, the initial 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 refinement is made with uniform mesh element size . The discrete approximation in each element is defined by the set of LegendreGaussRadau points. And for the initial refinement, the polynomialorder of th 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 differentialalgebraic constraints on the collocation points.
Begin: For .
Step 3. If , then continue (proceed to next element). is the approximation error of the th element in mesh.
Step 4 (adaptive hrefinement). If either , or , then refine the th element. refinement keeps fixed and adaptively constructs elements with varying size. It maintains refinement until the solution achieves the intermediate error and curvature ratio . The new number of element intervals is computed by formula
where is the ratio between current estimative maximum error and the specified error tolerance, is the ratio between current estimative maximum error and the intermediate error, and is the operator that rounds to the next integer. is an integer constant as a parameter for adjusting iterative speed, for example, . If there is some prior knowledge about state curve, can be adjusted for flections.
Formula (5.1) determines a dissatisfactory element interval should be subdivided into how many subelements byrefinement. 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 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. controls the growth in the number of element intervals. If 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 is small, the mesh will increase slowly, but the algorithm may require many iterations. An appropriate value of should make tradeoff between hrefinement and refinement.
The locations of the new requiredelements 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 is a constant to satisfy
The density function expressed by (5.2) is a probability density function (pdf) with fractional powerlaws. Here we adopt 1/3 power in the function. In fact, from the point of view of powerlaws in stochastic processes and/or fractal time series, 1/3 power is a typical heavytailed case. By heavytailed we mean that the density decays slowly. And heavytailed pdfs imply that is in wild randomness due to infinite or very large variance [36].
Let be the cumulative distribution function given by
The new mesh points are chosen by
Finally, it is noted that if , then no subintervals are created. Therefore, the minimum value for is 2.
Step 5 (adaptiverefinement). 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 on the th element is determined by where 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.
End: For .
Step 6 (return to Step 2). Our 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 adaptive algorithm procedure, adaptive 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 heavytailed probability distribution function, a slowly decayed autocorrelation function (ACF), and a power spectrum function (PSD) of 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, heavytailed 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 powerlawtype pdf. The key problem of 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 time series and powerlaws 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. Onedimensional random function with long memory is introduced to model the sea level fluctuation [40]. Topics in powerlaws 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 cyberphysical networking systems in [36], the contribution of simplicial approach to fractal structures by Cattani et al. in [42], and the work of Li et al. in [38].
6. 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 “takeoff” and “landing” segments. And the second example is a reusable launch vehicle reentry trajectory optimization. The experimentations compare three algorithms: GPM [12], adaptive PM based on curvature [17], and our integrated multicriterion adaptive PM. For conciseness, the three algorithms are denoted by , , , 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 opensource software GPOPS [43]. The transcribed NLP is solved by SNOPT [44, 45], using a 2.4GHz 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: and . In and , the maximum number of collocation points on each mesh is 20, and the upper bound of iteration times is 20.
Example 6.1. Consider the following hypersensitive OCP adapted from [46]. It is a Hamiltonian boundary value OCP in long time interval. Minimize the cost function subject to the dynamic constraint and the boundary conditions where is fixed. This problem has a “takeoff”, “cruise”, and “landing” structure where all of the interesting behaviors occur in the “takeoff” and “landing” segments. In this example, .
The specification of initial grid is listed in Table 1. The efficiency and accuracy of various methods are listed in Table 2. Figure 3 shows the state and control curve of solution on the final mesh for 3 method, at the accuracy tolerance grade . The rings denote the location of collocation points. From analysis of these records, some results are illustrated.


(a) State versus time
(b) Control versus time
First, adaptive methods and outperform method in computational efficiency and approximation error. The results are better at all three accuracy tolerance grades. method even cannot achieve accuracy tolerance forever. Because 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 method is better than method in computational efficiency. At all three accuracy tolerance grades, the count of collocation points for 3, 5 is less than the count of 3, 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 3 is less than the count of 5, but 3 method needs more iterations. From the comparison of experimental records, 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 adaptive algorithm. Meanwhile, the distribution of collocation points on the smooth segments is sparse. Compared with method and method, our method improves computational efficiency and reduces approximation error for this kind of nonsmooth OCP solving.
Example 6.2. Reusable launch vehicle reentry trajectory optimization is an OCP of maximizing the cross range during the atmospheric entry [1]. Minimize the cost function subject to the differential dynamic constraint and the boundary conditions where is the geocentric radius, is the longitude, is the latitude, is the speed, is the flight path angle, is the azimuth angle, and is the bank angle. The radius of the earth .
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 6 method, at the accuracy tolerance grade . The state and control curves are relatively smooth for this example. From analysis of these records, some results are illustrated.


(a) Height (km) versus time (s)
(b) Longitude (deg) versus time (s)
(c) Latitude (deg) versus time (s)
(d) Velocity (km/s) versus time (s)
(e) Flight path angle (deg) versus time (s)
(f) Azimuth angle (deg) versus time (s)
First, adaptive methods and are similar or better than method in this example. method can achieve accuracy tolerance and has the fewest number of collocation points, but it needs much iteration. The computational time of method is longer than adaptive methods. The low computational efficiency of method is due to the fact that the transcribed NLP is dense.
Second, our method outperforms method in computational efficiency. The count of collocation points for is less than the count of at three accuracy tolerance grades. The comparison illustrates that the distribution of collocation points via is better than the distribution via .
Third, the count of collocation points for 6 is less than the count of 4. Based on the comparison of experimental records, 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 methods andadaptive methods. methods are fit for solving problems which their solutions are smooth. adaptive methods perform similarly or better to methods for smooth problems. For nonsmooth problems, methods cannot accurately capture the discontinuities and rapid changes in the trajectory. They only achieve lowaccuracy even with much iteration. When a highaccuracy tolerance is required to satisfy, 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 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.
7. 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 adaptive PM, in order to further improve the performance of computation and approximation. For this purpose, we first devised an OCP solving framework of adaptive PM. We then designed a multicriterion 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 adaptive PM based on curvature and effectively improves the convergence rate of computation and the accuracy of solution. Meanwhile, it enhances the applicability of 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 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 (, , and ) 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 powerlaws to guide the choice of density function for various problems. This point itself is an issue worth discussing in our future work.
Appendix
Basis of the Pseudospectral Method
Much of the details of pseudospectral methods (PMs) are extensively described elsewhere [2, 15, 16, 34]; here, we briefly summarize the basis of the PM. The basic principle of a PM is to solve the continuous optimal control problem (OCP) by discretizing it to the NLP problem. The whole process mainly includes three steps. First, the horizon of a continuous OCP is mapped to a finite spectral interval . Second, a PM is applied to discretize the problem. The discretized problem generates a sequence of approximate problems parameterized by the size of the grid. In the third step, each discretized problem in this sequence is solved by a globally convergent sequential quadratic programming (SQP) technique. In the PMs, the functions are approximated by Lagrange interpolating polynomials of order in which interpolation occurs at Gaussian quadrature points. There are three different families of Gauss quadrature points known as LegendreGauss, LegendreGaussRadau, and LegendreGaussLobatto. Gauss PM (GPM) [12], Radau PM (RPM) [13], and Lobatto PM (LPM) [14] are the three most welldeveloped 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 LegendreGauss (LG) points across the interval.
The standard interval considered here is denoted as . By using a linear transformation, the actual time can be expressed as a function of via where and stand for the initial and final time, respectively.
The state is approximated using a basis of Lagrange interpolating polynomials, , as follows: where where are the LG points belonging to the set.
An approximation to the control, , is formed with a basis of Lagrange interpolating polynomials, , as follows:
The differential dynamic constraints should be posed to ensure that the entire discrete state satisfies the dynamic equations. The lefthand side of the dynamic equations is approximated by differentiating the state approximation of (A.1) at the LG points as follows: The differentiation matrix, , can be determined offline from dynamic equations:
The collocation equations require (A.5) to be equal to the righthand 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: where 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:
The transcribed NLP can be solved by welldeveloped algorithms, such as SQP, interior point method, particle swarm optimization and so on.
Acknowledgments
Great thanks go to the reviewers for valuable comments on our revision of the paper. This work was supported by the National Natural Science Foundation of China (61005077), the National Basic Research Program of China (6138101001), and the National Defense Foundation of China (403060103).
References
 J. T. Betts, Practical Methods for Optimal Control and Estimation Using Nonlinear Programming, Advances in Design and Control, SIAM Press, Philadelphia, Pa, USA, 2nd edition, 2009. View at: Publisher Site
 D. Garg, M. Patterson, W. W. Hager, A. V. Rao, D. A. Benson, and G. T. Huntington, “A unified framework for the numerical solution of optimal control problems using pseudospectral methods,” Automatica, vol. 46, no. 11, pp. 1843–1851, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 K. W and B. N, “Pseudospectral optimal control theory makes debut flightsaves NASA $1 M in under 3 hrs,” SIAM News, vol. 40, no. 7, pp. 1–2, 2007. View at: Google Scholar
 M. Li, “An optimal controller of an irregular wave maker,” Applied Mathematical Modelling, vol. 29, no. 1, pp. 55–63, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. Castellucci, M. Carlini, M. Guerrieri, and T. Honorati, “Stability and control for energy production parametric dependence,” Mathematical Problems in Engineering, vol. 2010, Article ID 842380, 21 pages, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Soler, A. Olivares, and E. Staffetti, “Hybrid optimal control approach to commercial aircraft trajectory planning,” Journal of Guidance, Control, and Dynamics, vol. 33, no. 3, pp. 985–991, 2010. View at: Publisher Site  Google Scholar
 L. Xie, P. Yang, T. Yang et al., “DualEKFbased realtime celestial navigation for lunar rover,” Mathematical Problems in Engineering, vol. 2012, Article ID 578719, 16 pages, 2012. View at: Publisher Site  Google Scholar
 S. H. Pourtakdoust, M. Kiani, and A. Hassanpour, “Optimal trajectory planning for flight through microburst wind shears,” Aerospace Science and Technology, vol. 15, no. 7, pp. 567–576, 2011. View at: Publisher Site  Google Scholar
 M. Li, “An iteration method to adjusting random loading for a laboratory fatigue test,” International Journal of Fatigue, vol. 27, no. 7, pp. 783–789, 2005. View at: Publisher Site  Google Scholar
 I. Babuška and M. Suri, “The $p$ and $h$$p$ versions of the finite element method, basic principles and properties,” SIAM Review, vol. 36, no. 4, pp. 578–632, 1994. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 W. F. Mitchell and M. A. McClain, “A survey of hpadaptive strategies for elliptic partial differential equations,” in Recent Advances in Computational and Applied Mathematics, pp. 227–258, Springer, Amsterdam, The Netherlands, 2011. View at: Google Scholar
 D. A. Benson, G. T. Huntington, T. P. Thorvaldsen, and A. V. Rao, “Direct trajectory optimization and costate estimation via an orthogonal collocation method,” Journal of Guidance, Control, and Dynamics, vol. 29, no. 6, pp. 1435–1440, 2006. View at: Publisher Site  Google Scholar
 D. Garg, M. A. Patterson, C. L. Darby et al., “Direct trajectory optimization and costate estimation of finitehorizon and infinitehorizon optimal control problems using a Radau pseudospectral method,” Computational Optimization and Applications, vol. 49, no. 2, pp. 335–358, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 G. Elnagar, M. A. Kazemi, and M. Razzaghi, “The pseudospectral Legendre method for discretizing optimal control problems,” IEEE Transactions on Automatic Control, vol. 40, no. 10, pp. 1793–1796, 1995. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, New York, NY, USA, 1998. View at: Publisher Site
 Q. Gong, F. Fahroo, and I. M. Ross, “Spectral algorithm for pseudospectral methods in optimal control,” Journal of Guidance, Control, and Dynamics, vol. 31, no. 3, pp. 460–471, 2008. View at: Publisher Site  Google Scholar
 C. L. Darby, W. W. Hager, and A. V. Rao, “Direct trajectory optimization using a variable loworder adaptive pseudospectral method,” Journal of Spacecraft and Rockets, vol. 48, no. 3, pp. 433–445, 2011. View at: Publisher Site  Google Scholar
 C. L. Darby, W. W. Hager, and A. V. Rao, “An $hp$adaptive pseudospectral method for solving optimal control problems,” Optimal Control Applications & Methods, vol. 32, no. 4, pp. 476–502, 2011. View at: Publisher Site  Google Scholar
 A. K. Patra, A. Laszloffy, and J. Long, “Data structures and load balancing for parallel adaptive hp finiteelement methods,” Computers and Mathematics with Applications, vol. 46, no. 1, pp. 105–123, 2003. View at: Google Scholar
 L. Demkowicz, W. Rachowicz, and Ph. Devloo, “A fully automatic $hp$adaptivity,” Journal of Scientific Computing, vol. 17, no. 1–4, pp. 117–142, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 L. Demkowicz, “Computing with hpadaptive finite elements,” in One and Two Dimensional Elliptic and Maxwell Problems, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2007. View at: Google Scholar  Zentralblatt MATH
 P. Solín, J. Červený, and I. Doležel, “Arbitrarylevel hanging nodes and automatic adaptivity in the $hp$FEM,” Mathematics and Computers in Simulation, vol. 77, no. 1, pp. 117–132, 2008. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 R. Tews and W. Rachowicz, “Application of an automatic $hp$ adaptive finite element method for thinwalled structures,” Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 2126, pp. 1967–1984, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 P. Solín and M. Kuraz, “Solving the nonstationary Richards equation with adaptive hpFEM,” Advances in Water Resources, vol. 34, no. 9, pp. 1062–1081, 2011. View at: Publisher Site  Google Scholar
 A. Szymczak, A. Paszýnska, M. Paszýnski et al., “Anisotropic 2D mesh adaptation in hpadaptive FEM,” Procedia Computer Science, vol. 4, pp. 1818–1827, 2011. View at: Publisher Site  Google Scholar
 L. Galvão, M. Gerritsma, and B. D. Maerschalck, “$hp$adaptive least squares spectral element method for hyperbolic partial differential equations,” Journal of Computational and Applied Mathematics, vol. 215, no. 2, pp. 409–418, 2008. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 H. Ben Dhia, L. Chamoin, J. T. Oden, and S. Prudhomme, “A new adaptive modeling strategy based on optimal control for atomictocontinuum coupling simulations,” Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 37–40, pp. 2675–2696, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. T. Oden and S. Prudhomme, “Control of modeling error in calibration and validation processes for predictive stochastic models,” International Journal for Numerical Methods in Engineering, vol. 87, no. 1–5, pp. 262–272, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 C. A. Dorao and H. A. Jakobsen, “$hp$adaptive least squares spectral element method for population balance equations,” Applied Numerical Mathematics, vol. 58, no. 5, pp. 563–576, 2008. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 C. A. Dorao, M. Fernandino, H. A. Jakobsen, and H. F. Svendsen, “hpAdaptive spectral element solver for reactor modeling,” Chemical Engineering Science, vol. 64, no. 5, pp. 904–911, 2009. View at: Publisher Site  Google Scholar
 Z. Lu, “Adaptive mixed finite element methods for parabolic optimal control problems,” Mathematical Problems in Engineering, vol. 2011, Article ID 217493, 21 pages, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. T. Oden and A. Patra, “A parallel adaptive strategy for $hp$ finite element computations,” Computer Methods in Applied Mechanics and Engineering, vol. 121, no. 1–4, pp. 449–470, 1995. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Ainsworth and B. Senior, “An adaptive refinement strategy for hpfinite element computations,” Applied Numerical Mathematics, vol. 26, no. 12, pp. 165–178, 1998. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 G. T. Huntington, Advancement and Analysis of a Gauss Pseudospectral Transcription for Optimal Control Problems, Massachusetts Institute of Technology, Cambridge, Mass, USA, 2007.
 Y. Zhao and P. Tsiotras, “Density functions for mesh refinement in numerical optimal control,” Journal of Guidance, Control, and Dynamics, vol. 34, no. 1, pp. 271–277, 2011. View at: Publisher Site  Google Scholar
 M. Li and W. Zhao, “Visiting power laws in cyberphysical networking systems,” Mathematical Problems in Engineering, vol. 2012, Article ID 302786, 13 pages, 2012. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Li, “Fractal time series—a tutorial review,” Mathematical Problems in Engineering, vol. 2010, Article ID 157264, 26 pages, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Li, M. Scalia, and C. Toma, “Nonlinear time series: computations and applications,” Mathematical Problems in Engineering, vol. 2010, Article ID 101523, 5 pages, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Li and W. Zhao, “Quantitatively investigating the locally weak stationarity of modified multifractional Gaussian noise,” Physica A, vol. 391, no. 24, pp. 6268–6278, 2012. View at: Publisher Site  Google Scholar
 M. Li, C. Cattani, and S. Y. Chen, “Viewing sea level by a onedimensional random function with long memory,” Mathematical Problems in Engineering, vol. 2011, Article ID 654284, 13 pages, 2011. View at: Publisher Site  Google Scholar
 M. Li, “Approximating ideal filters by systems of fractional order,” Computational and Mathematical Methods in Medicine, vol. 2012, Article ID 365054, 6 pages, 2012. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 C. Cattani, E. Laserra, and I. Bochicchio, “Simplicial approach to fractal structures,” Mathematical Problems in Engineering, vol. 2012, Article ID 958101, 21 pages, 2012. View at: Publisher Site  Google Scholar
 A. V. Rao, D. A. Benson, and C. L. Darby, “User's manual for GPOPS: a Matlab software for solving optimal control problems using pseudospectral methods,” 2011, http://www.gpops.org,. View at: Google Scholar
 P. E. Gill, W. Murray, and M. A. Saunders, “SNOPT: an SQP algorithm for largescale constrained optimization,” SIAM Review, vol. 47, no. 1, pp. 99–131, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 K. Holmström, A. O. Göran, and M. M. Edvall, Users Guide For TOMLAB /SNOPT, Mälardalen University, Department of Mathematics and Physics, Västerås, Sweden, 2006.
 A. V. Rao, “Application of a dichotomic basis method to performance optimization of supersonic aircraft,” Journal of Guidance, Control, and Dynamics, vol. 23, no. 3, pp. 570–573, 2000. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2012 Hongfu Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.