Mathematical Problems in Engineering

Volume 2017, Article ID 2184658, 7 pages

https://doi.org/10.1155/2017/2184658

## Adaptive Mesh Iteration Method for Trajectory Optimization Based on Hermite-Pseudospectral Direct Transcription

^{1}Air and Missile Defense College, Air Force Engineering University, Xi’an 710051, China^{2}College of Education, Shaanxi Normal University, Xi’an 710062, China^{3}Nanyang No. 2 High School, Nanyang 473000, China

Correspondence should be addressed to Tao Liu; moc.361@redevoltl

Received 20 February 2017; Accepted 3 May 2017; Published 18 June 2017

Academic Editor: Vladimir Turetsky

Copyright © 2017 Humin Lei 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.

#### Abstract

An adaptive mesh iteration method based on Hermite-Pseudospectral is described for trajectory optimization. The method uses the Legendre-Gauss-Lobatto points as interpolation points; then the state equations are approximated by Hermite interpolating polynomials. The method allows for changes in both number of mesh points and the number of mesh intervals and produces significantly smaller mesh sizes with a higher accuracy tolerance solution. The derived relative error estimate is then used to trade the number of mesh points with the number of mesh intervals. The adaptive mesh iteration method is applied successfully to the examples of trajectory optimization of Maneuverable Reentry Research Vehicle, and the simulation experiment results show that the adaptive mesh iteration method has many advantages.

#### 1. Introduction

There are general two kinds of methods to solve the optimal control problems, indirect method and direct method [1–3]. The indirect method uses the Pontryagin minimum principle and the first-order optimal necessary conditions to convert the optimal problem to two or more points boundary value problem. This method could obtain the exact optimal solution. The advantages of the indirect method are obtaining the high precision solution, satisfying the optimality of the first-order necessary conditions. The disadvantages are the small convergence radius and highly sensitive costate which is difficult to estimate. In recently years, direct methods have been widely applied for the numerical solution of nonlinear optimal control problems [4, 5]. The state and control are discretized at a series of suitable points in a direct method, then the continuous-time optimal control is converted into a finite dimensional nonlinear programming problem (NLP), and the next procedure is using the NLP solver software to solve the NLP [6].

Among the popular direct methods are the pseudospectral method and Hermite-Simpson method [7, 8]. The pseudospectral method has the advantage of high rate of convergence [9–11], so this paper develops the theory about the Legendre pseudospectral method. There are also some considerable interest in developing theory related to Hermite-Simpson method due to its reasonable accuracy with highly sparse Hessian and constraint Jacobians matrix. Williams provides a framework for arbitrary order and arbitrary number of intervals for implementation on digital computers [12]. That method allows the trading of the mesh points within each interval with the number of intervals; thus, better accuracy can be achieved by increasing mesh points for smooth problems, while increasing the number of intervals to achieve better accuracy for nonsmooth problems. However, it is a fact that smooth regions and nonsmooth regions together exist in one solution of the problem, so it is difficult to trade the number of mesh points with the number of mesh intervals when solving a complicated problem. Motived by the desire to trade the number of mesh points with the number of mesh intervals, both the number of mesh points within each mesh interval and the number of mesh intervals are allowed to vary in the method described in this paper. Furthermore, the method also can improve computational efficiency by reducing the size of the mesh.

#### 2. Trajectory Optimization Problems

Without loss of generality, consider the following optimal control problems with inequality path constraints:

subject to the constraintswhere the term denotes the state and the term denotes the control. In (1)-(2), the time domain is transformed from the time domain by the following affine transformation [12]:where the terms and represent initial time and terminal time, respectively. The basic idea of the approach in this paper is based on interpolating functions for state and costate on Legendre-Gauss-Lobatto (LGL) quadrature nodes [4, 7]. As the LGL nodes points are distributed over the interval , so it will be useful to transform the time domain.

#### 3. Adaptive Mesh Iteration Methodology

##### 3.1. Numerical Discretization

The domain is divided into mesh subintervals when using mesh refinement. Then we have

The mesh points have the property . The state in the subintervals is approximated by the Hermite interpolating polynomial with* n*th order [11]

Then the cost function is approximated by Gauss-Lobatto quadrature rule as follows:where the term is the same as in [12].

##### 3.2. Approximation of Solution Error

The estimate method for the relative error is similar to the error estimate obtained for numerically solving a different equation through using the modified Euler Runge-Kutta scheme [13, 14]. Suppose the NLP on mesh , , with HLGL points has been solved. The ensuing mesh with HLGL points , where and . Assume further that are the values of the state approximation at . We then have

The absolute error and the relative error approximations at of the state are defined, respectively, as

The maximum relative error in is then defined as

##### 3.3. Nonsmooth Solution Location

If a mesh interval has met the accuracy tolerance, that is, , where is the desired relative tolerance, then mesh size is reduced by decrease of the number of collocation points or merging adjacent mesh intervals; otherwise the mesh size needs to be modified by increasing points or dividing that mesh interval into several subintervals. Let be the curvature of the th component of the state in mesh interval , as

Let and be the mean and maximum value of , respectively. Then define as the ratio of the maximum to the mean curvature:

##### 3.4. Increasing the Number of Mesh Points within a Mesh Interval

When and if , where is a user-defined parameter, the curvature is considered uniform in this mesh interval and then the number of collocation points should be increased in interval . Let and denote the number of collocation points in interval at mesh and , respectively, where is the mesh refinement iteration. The number of points at mesh is calculated by the formula

It is noted in (12) that the ratio of the maximum to the error tolerance has a direct effect on the increases of the polynomial degree in mesh interval. An upper limit is set for the maximum allowable polynomial degree to make sure that the number of collocation points does not grow to an unreasonably large value. If (i.e., exceeds the maximum allowable polynomial degree), then the mesh interval must be divided into equally spaced subintervals.

##### 3.5. Generation of New Mesh Segment

Assume and , and then the th mesh interval should be refined. The following procedure is the strategy for mesh interval division. Firstly, the predicted polynomial determines the number of all the collocation points in the new subinterval. Secondly, the number of collocation points should be no fewer than the minimum allowable number. In other words, whenever dividing a mesh interval, each interval will contain at least collocation points. Third, the new number of mesh intervals, , is given by the formulawhere is a user-defined positive integer. In this process, it is ensured that the number of new intervals should be at least two. Thus, the number of new subintervals, denoted as , can be rewritten as

##### 3.6. Reducing the Number of Collocation Points in a Mesh Interval

The relative error of the mesh interval is less than the desired relative tolerance, and , and then the number of collocation points should be decreased. The number of points at mesh is calculated by the formula

##### 3.7. Merging Adjacent Mesh Subintervals

Before the adjacent mesh subintervals merging, it is necessary to decrease the number of each interval according to the method in Section 3.5 and then generally estimate the number of mesh interval points. If , the mesh interval and mesh interval cannot be merged because the highest polynomial orders of the two adjacent mesh intervals are not equal. All the matching points of the original two mesh intervals are combined, and the conditions for the merging of the two mesh subintervals are mainly three:(1)The two mesh subintervals must be adjacent.(2)The relative error estimates of the two grid intervals are not more than .(3)The relative error of the new mesh interval after the merger is not larger than .

##### 3.8. Penalty Function Method for Solving NLP

For convenience, the cost function and the path constraints are adjoined together, so that we obtain the augmented cost function where the term represents penalty factor and . When the penalty factor is greater than a certain threshold, the solution of problem is equal to the solution of original problem. The max operator has a great difference on the time spent of solving the problem, so it is time-consuming to solve the discretized NLP by using the existing mature gradient-based optimization algorithm. This paper develops the smoothing approximation function method proposed in [15].

Obviously, the function is first-order differentiable with respect to and has a good approximation to the max operation. When the smoothing factor , then we have . The optimal index is converted into

The gradient-based NLP algorithm also needs to compute the gradient information of the objective function for the optimal variable. Define the Hamiltonian function of optimal control problem:

##### 3.9. Mesh Refinement Method

The schematic of adaptive mesh iteration method is shown in Figure 1. The adaptive mesh iteration method is summarized as follows.