An Integrated Multicriterion hp-Adaptive 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 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 ℎ-method, which uses many low-degree 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 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 𝑝-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 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 [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 LS-SEM 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 multiple-interval 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. 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𝑔(𝐱(𝑡),𝐮(𝑡),𝑡)𝑑𝑡 (2.1) subject to the dynamic constraints 𝑑𝐱(𝐱(𝑡)(𝑡)),𝑑𝑡=𝐟,𝐮,𝑡 (2.2) the inequality path constraints 𝐂(𝐱(𝑡),𝐮(𝑡),𝑡)≤0, (2.3) and the boundary conditions 𝜓𝐱𝑡₀,𝑡₀,𝐱𝑡_𝑓,𝑡_𝑓=0, (2.4) 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 𝑡_𝑘−1 and end at the mesh point 𝑡_𝑘. The time domain in each element 𝑘, 𝑡∈[𝑡_𝑘−1,𝑡_𝑘] is transformed to 𝜏∈[−1,1] by the formula 𝑡𝜏=2𝑡−_𝑘−1+𝑡_𝑘𝑡_𝑘−𝑡_𝑘−1,𝑡_𝑘−1<𝑡_𝑘 (2.5)𝜏₀≡−1, and 𝜏_𝑓≡1. Furthermore, the inverse transformation is given as 𝑡𝑡=^𝑘−𝑡_𝑘−1𝜏+𝑡_𝑘−1+𝑡_𝑘2. (2.6)

The nonlinear continuous-time 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 𝐽=𝜙𝐱⁽¹⁾𝜏₀,𝑡₀,𝐱^(𝐾)𝜏_𝑓,𝑡_𝐾+^𝐾_𝑘=1𝑡_𝑘−𝑡_𝑘−12^𝜏𝑓_𝜏0𝑔𝐱^(𝑘)(𝜏),𝐮^(𝑘)(𝜏),𝜏;𝑡₀,𝑡_𝐾𝑑𝜏. (2.7) Next, the dynamics constraints, the boundary conditions, and the inequality path constraints are given, respectively, as 𝑑𝐱^(𝑘)(𝜏)=𝑡𝑑𝜏^𝑘−𝑡_𝑘−12𝐟𝐱^(𝑘)(𝜏),𝐮^(𝑘)(𝜏),𝜏;𝑡₀,𝑡_𝐾,()𝜓𝐱𝑘=1,2,…,𝐾₍₁₎𝜏₀,𝑡₀,𝐱^(𝐾)𝜏_𝑓,𝑡_𝐾𝑡=0,_𝑘−𝑡_𝑘−12𝐂𝐱^(𝑘)(𝜏),𝐮^(𝑘)(𝜏),𝜏;𝑡₀,𝑡_𝐾≤0,(𝑘=1,2,…,𝐾) (2.8)(2.9)(2.10) and continuity constraints at element interfaces 𝐱^(𝑘)𝜏_𝑓=𝐱^(𝑘+1)𝜏₀,(𝑘=1,2,…,𝐾−1). (2.11)

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.

Figure 1 

Flow chart of ℎ𝑝-adaptive pseudospectral method.

hp-adaptive 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 [𝑡_𝑘−1,𝑡_𝑘],𝑡₀≡𝑡_𝑘−1,𝑡_𝑓≡𝑡_𝑘, 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 𝐱(𝜏)≈𝐗(𝜏)=^𝑁_𝑖=0𝐗𝜏_𝑖𝐿_𝑖(𝜏),𝐮(𝜏)≈𝐔(𝜏)=^𝑁_𝑖=0𝐔𝜏_𝑖𝐿_𝑖(𝜏), (3.1)(3.2) where (𝜏₁,…𝜏_𝑁)  are collocation points, plus the start point 𝜏₀≡−1, 𝑁 is the number of collocation points. 𝐿_𝑖(𝜏),(𝑖=0,1,…,𝑁), is a basis of Lagrange polynomial with a variable degree. Then differential constraint of (2.8) can be transcribed into the algebraically form 𝐿̇^𝑁_𝑖=0𝑖𝜏_𝑖𝑥^𝑖−𝑡_𝑓−𝑡₀2𝐟𝐗_𝑖,𝐔_𝑖,𝜏;𝑡₀,𝑡_𝑓=0, (3.3) where 𝐗_𝑖=𝐗(𝜏_𝑖),𝐔_𝑖=𝐔(𝜏_𝑖),(𝑖=1,…,𝑁). The initial state 𝐗₀≡𝐗(𝜏₀), and terminal state 𝐗_𝑓≡𝐗(𝜏_𝑓) satisfies Gauss quadrature formula 𝐗_𝑓≡𝐗₀+𝑡_𝑓−𝑡₀2^𝑁_𝑖=1𝜔_𝑖𝐟𝐗_𝑖,𝐔_𝑖,𝜏_𝑖;𝑡₀,𝑡_𝑓=0. (3.4) Differentiating 𝐗(𝜏) in (3.1) with respect to 𝜏, we get 𝑑𝐗(𝜏)≡𝐗̇(𝜏)=𝑑𝜏^𝑁_𝑖=0𝐗𝜏_𝑖𝐿̇_𝑖(𝜏)=^𝑁_𝑖=0𝐗𝜏_𝑖𝑑𝐿_𝑖(𝜏).𝑑𝜏 (3.5) The continuous objective function of (2.7) can be approximated via LGR quadrature 𝐽≡𝜙𝐗₀,𝑡₀,𝐗_𝑓,𝑡_𝑓+𝑡_𝑓−𝑡₀2^𝑁_𝑖=0𝜔_𝑖𝑔𝐗_𝑖,𝐔_𝑖,𝜏;𝑡₀,𝑡_𝑓, (3.6) where 𝜔_𝑖 is the weight coefficients of Gauss quadrature. And the boundary conditions of (2.9) is 𝜓𝐗₀,𝑡₀,𝐗_𝑓,𝑡_𝑓=0. (3.7) Combining the state and control variances, the trajectory is constrained by 𝐂𝐗_𝑖,𝐔_𝑖,𝜏_𝑖;𝑡₀,𝑡_𝑓≤0. (3.8) 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 [𝑡_𝑘−1,𝑡_𝑘] in mesh, (𝑘=1,2,…,𝐾), we firstly transform the interval of variable 𝑡 into the interval of 𝜏∈[𝜏₀,𝜏_𝑓]. Then the approximation polynomials of terminal state are given as 𝐗^𝑘₀=𝐗₀^𝑘−1+^𝑁_𝑖=0𝐟𝐗_𝑖,𝐔_𝑖,𝜏_𝑖;𝑡_𝑘−1,𝑡_𝑘=𝐗₀^𝑘−1+^𝜏𝑓_𝜏0𝐗̇(𝜏)𝑑𝜏−^𝜏𝑓_𝜏0𝐗̇(𝜏)𝑑𝜏+^𝑁_𝑖=0𝐟𝐗_𝑖,𝐔_𝑖,𝜏_𝑖;𝑡_𝑘−1,𝑡_𝑘=𝐗^𝑘₀−^𝜏𝑓_𝜏0^𝑁_𝑖=0𝐗𝜏_𝑖𝐿̇_𝑖(𝜏)𝑑𝜏−^𝑁_𝑖=0𝐟𝐗_𝑖,𝐔_𝑖,𝜏_𝑖;𝑡_𝑘−1,𝑡_𝑘, (3.9) where 𝑁 is the number of collocation points in the 𝑘th element, (𝑖=0,1,…,𝑁). Hence the local error of the element is given as 𝐗^𝑘₀−𝐗^𝑘₀=^𝜏𝑓_𝜏0^𝑁_𝑖=0𝐗𝜏_𝑖𝐿̇_𝑖(𝜏)𝑑𝜏−^𝑁_𝑖=0𝐟𝐗_𝑖,𝐔_𝑖,𝜏_𝑖;𝑡_𝑘−1,𝑡_𝑘. (3.10) Furthermore, we can define the absolute error of the 𝑘th element interval as 𝜀_𝑘=1ℎ_𝑘^𝜏𝑓_𝜏0^𝑁_𝑖=0𝐗𝜏_𝑖𝐿̇_𝑖(𝜏)𝑑𝜏−^𝑁_𝑖=0𝐟𝐗_𝑖,𝐔_𝑖,𝜏_𝑖;𝑡_𝑘−1,𝑡_𝑘, (3.11) where ℎ_𝑘=𝑡_𝑘−𝑡_𝑘−1 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 𝜂^(𝑘)=max_𝑗𝜀_𝑘𝜔_𝑗, (3.12) where 𝜔_𝑗=1𝑁^𝑁_𝑖=0𝑥_𝑗,𝑖. (3.13)𝜔_𝑗 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 knowledge-based 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 𝑟≈1+𝛼, (4.1) 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 𝑝_𝑘+2≤𝑟_𝑘, then 𝑝-refinement is performed, otherwiseℎ-refinement is adopted. Furthermore, let 𝑝⁰_𝑘=𝑟⁰_𝑘−2 be an optimal choice for initial 𝑝-refinement.

Second, intermediate error can be adopted as a criterion for ℎ𝑝-adaptive strategy. Define an intermediate error 𝜂_𝐼=𝛾𝜂, (4.2) 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 𝑘 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 𝜅^(𝑘)(𝜏)=𝑋̈_𝑗^(𝑘)(𝜏)1+𝑋̇_𝑗^(𝑘)(𝜏)^23/2, (4.3) where 𝑋_𝑗^(𝑘)(𝜏) is the 𝑗th component of the state approximation in the 𝑘th element [17]. Assume 𝜅^(𝑘)_max and 𝜅^(𝑘) be the maximum and mean value of 𝜅^(𝑘)(𝜏). Furthermore, the ratio of the maximum curvature to the mean curvature is defined 𝑟_𝑘=𝜅^(𝑘)_max𝜅^(𝑘). (4.4)

For the element 𝑘 in mesh, if 𝑟_𝑘<𝑟_max (where 𝑟_max>0 is a user-defined parameter), then 𝑝-refinement is adopted and a larger degree polynomial is used to obtain a better approximation. If 𝑟_𝑘≥𝑟_max, 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.

Figure 2 

Flow chart of multicriterion ℎ𝑝-adaptive strategy.

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 𝜂=0.01 (1% error). Then the intermediate error 𝜂_𝐼 is specified, for example, 𝜂_𝐼=0.1. And 𝑟_max is specified, for example 𝑟_max=8.
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 Legendre-Gauss-Radau points. And for the initial 𝑝-refinement, the polynomial-order of 𝑘th element is determined via the assessment of local solution regularity, it is given as 𝑝⁰_𝑘=𝑟⁰_𝑘−2.

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.

Begin: For  𝑘=1,2,…,𝐾.

Step 3. If 𝜂^(𝑘)≤𝜂, then continue (proceed to next 𝑘+1 element). 𝜂^(𝑘) is the approximation error of the 𝑘th element in mesh.

Step 4 (adaptive h-refinement). If either 𝜂^(𝑘)>𝜂_𝐼, or 𝑟_𝑘≥𝑟_max, 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 𝑟_𝑘<𝑟_max. The new number of element intervals is computed by formula 𝑛_𝑘=log₁₀𝜂^(𝑘)𝜂+𝑐_ℎ×log_𝑒𝜂^(𝑘)𝜂_𝐼=log₁₀𝜂^(𝑘)𝜂+𝑐_ℎ×log_𝑒10×log₁₀𝜂^(𝑘)𝜂_𝐼=log₁₀𝜂^(𝑘)𝜂+log₁₀𝜂^(𝑘)𝜂_𝐼^{𝑐ℎ×log𝑒10}=log₁₀𝜂^(𝑘)𝜂×𝜂^(𝑘)𝜂_𝐼^{𝑐ℎ×log𝑒10}, (5.1) 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, 𝑐_ℎ=2. 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 byℎ-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 ℎ-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 h-refinement and 𝑝-refinement.
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 𝜌(𝜏)=𝑐𝜅(𝜏)^1/3, (5.2) where 𝑐 is a constant to satisfy ⁺¹₋₁𝜌(𝜁)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 heavy-tailed 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 𝐹(𝜏) be the cumulative distribution function given by 𝐹(𝜏)=^𝜏₋₁𝜌(𝜁)d𝜁. (5.4) The 𝑛_𝑘 new mesh points are chosen by 𝐹(𝜏)=𝑖−1𝑛_𝑘,1≤𝑖≤𝑛_𝑘+1. (5.5) Finally, it is noted that if 𝑛_𝑘=1, then no subintervals are created. Therefore, the minimum value for 𝑛_𝑘 is 2.

Step 5 (adaptive𝑝-refinement). 𝑝-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 𝑝_𝑘=𝑝_𝑘+log₁₀𝜂^(𝑘)𝜂+log_𝑒𝜂^(𝑘)𝜂_𝐼=𝑝_𝑘+log₁₀𝜂^(𝑘)𝜂+log_𝑒10×log₁₀𝜂^(𝑘)𝜂_𝐼=𝑝_𝑘+log₁₀𝜂^(𝑘)𝜂+log₁₀𝜂^(𝑘)𝜂_𝐼^log𝑒10=𝑝_𝑘+log₁₀𝜂^(𝑘)𝜂×𝜂^(𝑘)𝜂_𝐼^log𝑒10, (5.6) 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  𝑘=1,2,…,𝐾.

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 heavy-tailed probability distribution function, a slowly decayed autocorrelation function (ACF), and a power spectrum function (PSD) of 1/𝑓 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 ℎ-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 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 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 “take-off” and “landing” segments. And the second example is a reusable launch vehicle re-entry 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 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: 𝑟_max=8 and 𝛾=10. 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 1𝐽=2^𝑡𝑓₀𝑥²+𝑢²𝑑𝑡, (6.1) subject to the dynamic constraint 𝑥̇=−𝑥³+𝑢, (6.2) and the boundary conditions 𝑥(0)=1.5,𝑥𝑡_𝑓=1, (6.3) where 𝑡_𝑓 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, 𝑡_𝑓=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 shows the state and control curve of solution on the final mesh for ℎ𝑝⁽²⁾-3 method, at the accuracy tolerance grade 𝜂=1𝑒−3. The rings denote the location of collocation points. From analysis of these records, some results are illustrated.

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.

(a) State 𝑥 versus time 𝑡(s)

(b) Control 𝑢 versus time 𝑡(s)

Figure 3 

State and control curve on the final mesh of Example 6.1 for ℎ𝑝(2)-3 and 𝜂=1𝑒−3.

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 1𝑒−3 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 “take-off” 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 re-entry trajectory optimization is an OCP of maximizing the cross range during the atmospheric entry [1]. Minimize the cost function 𝐽=−𝜙𝑡_𝑓, (6.4) subject to the differential dynamic constraint ̇̇𝐷𝑟=𝑣sin𝛾,𝑣=𝑚̇−𝑔sin𝛾,𝜃=𝑣cos𝛾sin𝜓𝑟,𝛾̇=𝐿cos𝜎−𝑔𝑚𝑣𝑣−𝑣𝑟̇cos𝛾,𝜙=𝑣cos𝛾cos𝜓𝑟,𝜓̇=𝐿sin𝜎+𝑚𝑣cos𝛾𝑣cos𝛾sin𝜓tan𝜙𝑟, (6.5) and the boundary conditions 𝑟(0)=79248+𝑅_𝑒m,𝑟𝑡_𝑓=24384+𝑅_𝑒𝜃(0)𝑡m,=0deg,𝜃_𝑓𝜙(0)𝑡=Free,=0deg,𝜙_𝑓𝑣(0)𝑡=Free,=7803m/s,𝑣_𝑓𝛾(0)𝑡=762m/s,=−1deg,𝛾_𝑓𝜓(0)𝑡=−5deg,=90deg,𝜓_𝑓=Free, (6.6) 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 𝑅_𝑒=6371203.92m.

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 𝜂=1𝑒−4. The state and control curves are relatively smooth for this example. From analysis of these records, some results are illustrated.

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.

(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)

Figure 4 

State and control curve on the final mesh of Example 6.2 for ℎ𝑝(2)-6 and 𝜂=1𝑒−4.

First, ℎ𝑝-adaptive methods ℎ𝑝⁽¹⁾ and ℎ𝑝⁽²⁾ are similar or better than 𝑝-method in this example. 𝑝-method can achieve 1𝑒−5 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 andℎ𝑝-adaptive 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 low-accuracy even with much iteration. When a high-accuracy 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 (𝛾, 𝑟_max, 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 power-laws 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 [−1,1]. 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 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 𝑡 can be expressed as a function of 𝜏 via 𝑡𝑡=^𝑓−𝑡₀𝜏+𝑡_𝑓+𝑡₀2,𝜏=2𝑡𝑡_𝑓−𝑡₀−𝑡_𝑓+𝑡₀𝑡_𝑓−𝑡₀, (A.1) where 𝑡₀ and 𝑡_𝑓 stand for the initial and final time, respectively.

The state is approximated using a basis of 𝑁+1 Lagrange interpolating polynomials, 𝐿_𝑖(𝜏)(𝑖=0,…,𝑁), as follows: 𝐱(𝜏)≈𝐗(𝜏)=^𝑁_𝑖=0𝐗𝜏_𝑖𝐿_𝑖(𝜏), (A.2) where 𝐿_𝑖(𝜏)=^𝑁_{𝑗=0,𝑗≠𝑖}𝜏−𝜏_𝑗𝜏_𝑖−𝜏_𝑗, (A.3) where 𝜏_𝑖,(𝑖=0,…,𝑁) are the LG points belonging to the set.

An approximation to the control, 𝐔(𝜏), is formed with a basis of 𝑁 Lagrange interpolating polynomials, 𝐿_𝑖(𝜏)(𝑖=1,…,𝑁), as follows: 𝐮(𝜏)≈𝐔(𝜏)=^𝑁_𝑖=1𝐔𝜏_𝑖𝐿_𝑖(𝜏). (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: 𝐱̇𝜏_𝑘≈𝐗̇𝜏_𝑘=^𝑁_𝑖=0𝐗𝜏_𝑖𝐿̇_𝑖𝜏_𝑘=^𝑁_𝑖=0𝐷_𝑘𝑖𝐗𝜏_𝑖,(𝑘=1,…,𝑁). (A.5) The differentiation matrix, 𝐷∈ℝ^{𝑁×(𝑁+1)}, can be determined offline from dynamic equations: 𝐷_𝑘𝑖=𝐿̇_𝑖𝜏_𝑘=1+𝜏_𝑘𝑃̇_𝑁𝜏_𝑘+𝑃_𝑁𝜏_𝑘𝜏_𝑘−𝜏_𝑖1+𝜏_𝑖𝑃̇_𝑁𝜏_𝑖+𝑃_𝑁𝜏_𝑖,𝑖≠𝑘1+𝜏_𝑖𝑃̈_𝑁𝜏_𝑖+2𝑃̇_𝑁𝜏_𝑖21+𝜏_𝑖𝑃̇_𝑁𝜏_𝑖+𝑃_𝑁𝜏_𝑖,𝑖=𝑘. (A.6)

The 𝑁 collocation equations require (A.5) to be equal to the right-hand side of the dynamic equations at the collocation points: ^𝑁_𝑖=0𝐷_𝑘𝑖𝐗𝜏_𝑖−𝑡_𝑓−𝑡₀2𝑓𝑋𝜏_𝑘,𝑈𝜏_𝑘,𝜏_𝑘;𝑡₀,𝑡_𝑓=0,(𝑘=1,…,𝑁). (A.7)

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: 𝐂𝐗_𝑘,𝐔_𝑘,𝜏_𝑘;𝑡₀,𝑡_𝑓≤0. (A.8)

Lastly, the integral term in the cost functional can be approximated with a Gauss quadrature as follows: 𝐽=Φ𝑋₀,𝑡₀,𝑋_𝑓,𝑡_𝑓+𝑡_𝑓−𝑡₀2^𝑁_𝑘=1𝜔_𝑘𝑔𝑋_𝑘,𝑈_𝑘,𝜏_𝑘;𝑡₀,𝑡_𝑓, (A.9) 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: min𝐹(𝐲)s.t.𝑔_𝑗(𝐲)ℎ≤0,𝑗=1,2,…,𝑝_𝑘(𝐲)=0,𝑘=1,2,…,𝑞. (A.10)

The transcribed NLP can be solved by well-developed 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).