# An Integrated Multicriterionhp-Adaptive Pseudospectral Method for Direct Optimal Control Problems Solving

Authors: Hongfu Liu; Shaofei Chen; Lincheng Shen; Jing Chen
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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. Existinghp-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 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 hp-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.


---

## Body

## 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 namedh-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 ofh-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.h
p-adaptive methods attract increasing amount of attentions in the research of FEM [11, 19–25]. Babuška discussed the mathematical properties of p- and 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 [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 tohp-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. Section2 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.


## 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(2.1)J=ϕ(x(t0),t0,x(tf),tf)+∫t0tfg(x(t),u(t),t)dt
subject to the dynamic constraints
(2.2)dxdt=f(x(t),u(t),t),
the inequality path constraints
(2.3)C(x(t),u(t),t)≤0,
and the boundary conditions
(2.4)ψ(x(t0),t0,x(tf),tf)=0,
where x(t)∈ℝm is the state, u(t)∈ℝn is the control, C(x(t),u(t),t)∈ℝp is the control and state constraints, t is time, t0 is the start time, and tf is the terminal time.Next, consider dividing the nonlinear continuous optimal control problem defined above intoK element intervals. Let t0<t1<t2<⋯<tK, where tK=tf is the terminal time.An elements k is defined to begin at the mesh point tk-1 and end at the mesh point tk. The time domain in each element k, t∈[tk-1,tk] is transformed to τ∈[-1,1] by the formula
(2.5)τ=2t-(tk-1+tk)tk-tk-1,(tk-1<tk)τ0≡-1, and τf≡1. Furthermore, the inverse transformation is given as
(2.6)t=[(tk-tk-1)τ+(tk-1+tk)]2.The nonlinear continuous-time optimal control problem on element[t0,tf], 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
(2.7)J=ϕ(x(1)(τ0),t0,x(K)(τf),tK)+∑k=1K[tk-tk-12∫τ0τfg(x(k)(τ),u(k)(τ),τ;t0,tK)dτ].
Next, the dynamics constraints, the boundary conditions, and the inequality path constraints are given, respectively, as
(2.8)dx(k)(τ)dτ=tk-tk-12f(x(k)(τ),u(k)(τ),τ;t0,tK),(k=1,2,…,K)(2.9)ψ(x(1)(τ0),t0,x(K)(τf),tK)=0,(2.10)tk-tk-12C(x(k)(τ),u(k)(τ),τ;t0,tK)≤0,(k=1,2,…,K)
and continuity constraints at element interfaces
(2.11)x(k)(τf)=x(k+1)(τ0),(k=1,2,…,K-1).


## 3. Framework of Direct Optimal Control Problems Solving Usinghp-Adaptive Pseudospectral Method

### 3.1. Flow Chart ofhp-Adaptive Pseudospectral Method

Pseudospectral method transcribes a continuous OCP into a nonlinear programming problem (NLP). The flow chart ofhp-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 hp-adaptive algorithm.Figure 1
Flow chart ofhp-adaptive pseudospectral method.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.


### 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 [tk-1,tk],t0≡tk-1,tf≡tk, 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
(3.1)x(τ)≈X(τ)=∑i=0NX(τi)Li(τ),(3.2)u(τ)≈U(τ)=∑i=0NU(τi)Li(τ),
where (τ1,…τN)  are collocation points, plus the start point τ0≡-1, N is the number of collocation points. Li(τ),(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
(3.3)∑i=0NL˙i(τi)xi-tf-t02f(Xi,Ui,τ;t0,tf)=0,
where Xi=X(τi),Ui=U(τi),(i=1,…,N). The initial state X0≡X(τ0), and terminal state Xf≡X(τf) satisfies Gauss quadrature formula
(3.4)Xf≡X0+tf-t02∑i=1Nωif(Xi,Ui,τi;t0,tf)=0.
Differentiating X(τ) in (3.1) with respect to τ, we get
(3.5)dX(τ)dτ≡X˙(τ)=∑i=0NX(τi)L˙i(τ)=∑i=0NX(τi)dLi(τ)dτ.
The continuous objective function of (2.7) can be approximated via LGR quadrature
(3.6)J≡ϕ(X0,t0,Xf,tf)+tf-t02∑i=0Nωig(Xi,Ui,τ;t0,tf),
where ωi is the weight coefficients of Gauss quadrature. And the boundary conditions of (2.9) is
(3.7)ψ(X0,t0,Xf,tf)=0.
Combining the state and control variances, the trajectory is constrained by
(3.8)C(Xi,Ui,τi;t0,tf)≤0.
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 thekth element interval [tk-1,tk] 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
(3.9)X^0k=X0k-1+∑i=0Nf(Xi,Ui,τi;tk-1,tk)=X0k-1+∫τ0τfX˙(τ)dτ-∫τ0τfX˙(τ)dτ+∑i=0Nf(Xi,Ui,τi;tk-1,tk)=X0k-[∫τ0τf(∑i=0NX(τi)L˙i(τ))dτ-∑i=0Nf(Xi,Ui,τi;tk-1,tk)],
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
(3.10)X0k-X^0k=∫τ0τf(∑i=0NX(τi)L˙i(τ))dτ-∑i=0Nf(Xi,Ui,τi;tk-1,tk).
Furthermore, we can define the absolute error of the kth element interval as
(3.11)εk=1hk|∫τ0τf(∑i=0NX(τi)L˙i(τ))dτ-∑i=0Nf(Xi,Ui,τi;tk-1,tk)|,
where hk=tk-tk-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
(3.12)η(k)=maxjεkωj,
where
(3.13)ωj=1N∑i=0N|xj,i|.ω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.


## 3.1. Flow Chart ofhp-Adaptive Pseudospectral Method

Pseudospectral method transcribes a continuous OCP into a nonlinear programming problem (NLP). The flow chart ofhp-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 hp-adaptive algorithm.Figure 1
Flow chart ofhp-adaptive pseudospectral method.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.


## 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 [tk-1,tk],t0≡tk-1,tf≡tk, 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
(3.1)x(τ)≈X(τ)=∑i=0NX(τi)Li(τ),(3.2)u(τ)≈U(τ)=∑i=0NU(τi)Li(τ),
where (τ1,…τN)  are collocation points, plus the start point τ0≡-1, N is the number of collocation points. Li(τ),(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
(3.3)∑i=0NL˙i(τi)xi-tf-t02f(Xi,Ui,τ;t0,tf)=0,
where Xi=X(τi),Ui=U(τi),(i=1,…,N). The initial state X0≡X(τ0), and terminal state Xf≡X(τf) satisfies Gauss quadrature formula
(3.4)Xf≡X0+tf-t02∑i=1Nωif(Xi,Ui,τi;t0,tf)=0.
Differentiating X(τ) in (3.1) with respect to τ, we get
(3.5)dX(τ)dτ≡X˙(τ)=∑i=0NX(τi)L˙i(τ)=∑i=0NX(τi)dLi(τ)dτ.
The continuous objective function of (2.7) can be approximated via LGR quadrature
(3.6)J≡ϕ(X0,t0,Xf,tf)+tf-t02∑i=0Nωig(Xi,Ui,τ;t0,tf),
where ωi is the weight coefficients of Gauss quadrature. And the boundary conditions of (2.9) is
(3.7)ψ(X0,t0,Xf,tf)=0.
Combining the state and control variances, the trajectory is constrained by
(3.8)C(Xi,Ui,τi;t0,tf)≤0.
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 thekth element interval [tk-1,tk] 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
(3.9)X^0k=X0k-1+∑i=0Nf(Xi,Ui,τi;tk-1,tk)=X0k-1+∫τ0τfX˙(τ)dτ-∫τ0τfX˙(τ)dτ+∑i=0Nf(Xi,Ui,τi;tk-1,tk)=X0k-[∫τ0τf(∑i=0NX(τi)L˙i(τ))dτ-∑i=0Nf(Xi,Ui,τi;tk-1,tk)],
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
(3.10)X0k-X^0k=∫τ0τf(∑i=0NX(τi)L˙i(τ))dτ-∑i=0Nf(Xi,Ui,τi;tk-1,tk).
Furthermore, we can define the absolute error of the kth element interval as
(3.11)εk=1hk|∫τ0τf(∑i=0NX(τi)L˙i(τ))dτ-∑i=0Nf(Xi,Ui,τi;tk-1,tk)|,
where hk=tk-tk-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
(3.12)η(k)=maxjεkωj,
where
(3.13)ωj=1N∑i=0N|xj,i|.ω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.


## 4. Multicriterionhp-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 anh-refinement or by a p-refinement. 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 forhp-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 asr. A method for estimating the local regularity of a solution on a given mesh was proposed in [33], it indicated that
(4.1)r≈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 ofhp-adaptive strategy, the regularity value rk of approximate solution on the kth element is readily available by initialization. Let pk0 and pk, respectively, denote the initial degree and current degree of basis polynomial in the k element. If pk+2≤rk, then p-refinement is performed, otherwiseh-refinement is adopted. Furthermore, let pk0=rk0-2 be an optimal choice for initial p-refinement.Second, intermediate error can be adopted as a criterion forhp-adaptive strategy. Define an intermediate error
(4.2)ηI=γη~,
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 elementk in mesh, if the estimated error η(k)>ηI, it needs detailed h-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 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.Third, curvature is another criterion forhp-adaptive strategy. The curvature of state curve describes the regularity of intermediate solution. The curvature of the jth component of the state in element k is given by
(4.3)κ(k)(τ)=|X¨j(k)(τ)||[1+X˙j(k)(τ)2]3/2|,
where Xj(k)(τ) is the jth component of the state approximation in the kth element [17]. Assume κmax(k) 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)rk=κmax(k)κ-(k).For the elementk in mesh, if rk<rmax (where rmax>0 is a user-defined parameter), then p-refinement is adopted and a larger degree polynomial is used to obtain a better approximation. If rk≥rmax, then h-refinement is adopted in this element.Furthermore, we design anhp-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.Figure 2
Flow chart of multicriterionhp-adaptive strategy.


## 5. Integrated Multicriterionhp-Adaptive Algorithm Procedure

In this section, we depict an integrated multicriterionhp-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 rmax is specified, for example rmax=8.
Consider the initialh-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 h0. The discrete approximation in each element is defined by the set of Legendre-Gauss-Radau points. And for the initial p-refinement, the polynomial-order of kth element is determined via the assessment of local solution regularity, it is given as pk0=rk0-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 Section3.3. The accuracy of the approximation in each element is assessed by calculating the differential-algebraic constraints on the collocation points.Begin: For  k=1,2,…,K.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 (adaptiveh-refinement).
If eitherη(k)>ηI, or rk≥rmax, then h-refine the kth element. h-refinement keeps p=p0 fixed and adaptively constructs elements with varying size. It maintains refinement until the solution achieves the intermediate error ηI and curvature ratio rk<rmax. The new number of element intervals is computed by formula
(5.1)nk=⌈log10(η(k)η~)+ch×loge(η(k)ηI)⌉=⌈log10(η(k)η~)+ch×loge10×log10(η(k)ηI)⌉=⌈log10(η(k)η~)+log10(η(k)ηI)ch×loge10⌉=⌈log10((η(k)η~)×(η(k)ηI)ch×loge10)⌉,
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. ch is an integer constant as a parameter for adjusting iterative speed, for example, ch=2. If there is some prior knowledge about state curve, ch can be adjusted for flections.
Formula (5.1) determines a dissatisfactory element interval should be subdivided into how many subelements byh-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. ch controls the growth in the number of element intervals. If ch 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 ch is small, the mesh will increase slowly, but the algorithm may require many iterations. An appropriate value of ch should make tradeoff between h-refinement and p-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
(5.2)ρ(τ)=cκ(τ)1/3,
where c is a constant to satisfy
(5.3)∫-1+1ρ(ζ)dζ=1.
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].
LetF(τ) be the cumulative distribution function given by
(5.4)F(τ)=∫-1τρ(ζ)dζ.
The nk new mesh points are chosen by
(5.5)F(τ)=i-1nk,1≤i≤nk+1.
Finally, it is noted that if nk=1, then no subintervals are created. Therefore, the minimum value for nk is 2.Step 5 (adaptivep-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
(5.6)p^k=pk+⌈log10(η(k)η~)+loge(η(k)ηI)⌉=pk+⌈log10(η(k)η~)+loge10×log10(η(k)ηI)⌉=pk+⌈log10(η(k)η~)+log10(η(k)ηI)loge10⌉=pk+⌈log10((η(k)η~)×(η(k)ηI)loge10)⌉,
where pk 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  k=1,2,…,K.Step 6 (return to Step2).
Ourhp-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 ourhp-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 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], 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: rmax=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.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
(6.1)J=12∫0tf(x2+u2)dt,
subject to the dynamic constraint
(6.2)x˙=-x3+u,
and the boundary conditions
(6.3)x(0)=1.5,x(tf)=1,
where tf 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, tf=1000.The specification of initial grid is listed in Table1. 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 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.Table 1
Initial grid used in Example6.1.
Strategy
Number of collocation points
Number of intervals
p
60
1
h
p
(
1
)-5, hp(2)-5
10
5
h
p
(
1
)-3, hp(2)-3
10
3Table 2
Summary of results for Example6.1 using various collocation strategies and accuracy tolerances.
η
~
Strategy
Approximate error
Collocation points
Elements
Mesh iterations
CPU time (s)
p
9
.
23
e
-
2
203
1
50
197.36
h
p
(
1
)-5
6.51
e
-
2
366
72
8
66.90
1
e
-
2
h
p
(
2
)-5
2.92
e
-
2
269
52
8
25.49
h
p
(
1
)-3
5.03
e
-
2
337
102
10
49.67
h
p
(
2
)-3
1.89
e
-
2
244
68
9
21.08
p
—
—
—
—
—
h
p
(
1
)-5
8.63
e
-
3
478
93
8
91.85
1
e
-
3
h
p
(
2
)-5
4.10
e
-
3
347
66
7
42.93
h
p
(
1
)-3
6.45
e
-
3
436
125
11
82.72
h
p
(
2
)-3
2.17
e
-
3
329
87
9
39.80
p
—
—
—
—
—
h
p
(
1
)-5
7.14
e
-
4
614
119
9
117.22
1
e
-
4
h
p
(
2
)-5
6.59
e
-
4
426
75
8
66.08
h
p
(
1
)-3
6.97
e
-
4
565
152
11
107.79
h
p
(
2
)-3
3.86
e
-
4
407
108
10
56.16State and control curve on the final mesh of Example6.1 for hp(2)-3 and η~=1e-3.
(a)
Statex versus time t(s)
(b)
Controlu versus time t(s)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 proposedhp(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 forhp-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 Figure3, the state and control curves rapid change on the “take-off” 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
(6.4)J=-ϕ(tf),
subject to the differential dynamic constraint
(6.5)r˙=vsinγ,v˙=Dm-gsinγ,θ˙=vcosγsinψr,γ˙=Lcosσmv-(gv-vr)cosγ,ϕ˙=vcosγcosψr,ψ˙=Lsinσmvcosγ+vcosγsinψtanϕr,
and the boundary conditions
(6.6)r(0)=79248+Rem,r(tf)=24384+Rem,θ(0)=0deg,θ(tf)=Free,ϕ(0)=0deg,ϕ(tf)=Free,v(0)=7803m/s,v(tf)=762m/s,γ(0)=-1deg,γ(tf)=-5deg,ψ(0)=90deg,ψ(tf)=Free,
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 Re=6371203.92m.The specification of initial grid is listed in Table3. 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.Table 3
Initial grid used in Example6.2.
Strategy
Number of collocation points
Number of intervals
p
50
1
h
p
(
1
)-6, hp(2)-6
10
6
h
p
(
1
)-4, hp(2)-4
10
4Table 4
Summary of results for Example6.2 using various collocation strategies and accuracy tolerances.
η
~
Strategy
Approximate error
Collocation points
Elements
Mesh iterations
CPU time (s)
p
7
.
61
e
-
3
58
1
4
70.88
h
p
(
1
)-6
5.12
e
-
3
90
18
2
19.55
1
e
-
3
h
p
(
2
)-6
3.55
e
-
3
70
12
3
18.05
h
p
(
1
)-4
6.81
e
-
3
113
32
4
32.55
h
p
(
2
)-4
4.19
e
-
3
83
18
4
19.66
p
9.17
e
-
4
75
1
12
88.94
h
p
(
1
)-6
4.08
e
-
4
151
29
6
33.76
1
e
-
4
h
p
(
2
)-6
3.23
e
-
4
122
21
4
24.89
h
p
(
1
)-4
7.99
e
-
4
202
56
6
63.71
h
p
(
2
)-4
5.83
e
-
4
153
34
4
40.24
p
9.74
e
-
5
90
1
16
99.50
h
p
(
1
)-6
8.02
e
-
5
214
40
6
63.64
1
e
-
5
h
p
(
2
)-6
7.50
e
-
5
183
32
5
53.44
h
p
(
1
)-4
8.93
e
-
5
293
70
5
89.46
h
p
(
2
)-4
7.69
e
-
5
251
51
5
84.88State and control curve on the final mesh of Example6.2 for hp(2)-6 and η~=1e-4.
(a)
Heighth (km) versus time t (s)
(b)
Longitudeθ (deg) versus time t (s)
(c)
Latitudeϕ (deg) versus time t (s)
(d)
Velocityv (km/s) versus time t (s)
(e)
Flight path angleγ (deg) versus time t (s)
(f)
Azimuth angleψ (deg) versus time t (s)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, ourhp(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 forhp-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 ofp-methods andhp-adaptive methods. p-methods are fit for solving problems which their solutions are smooth. hp-adaptive 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.


## 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 tohp-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 forhp-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 (γ, rmax, and ch) 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.



---

*Source: 760890-2012-11-12.xml*