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

Authors: Humin Lei; Tao Liu; Deng Li; Jikun Ye; Lei Shao
Journal: Mathematical Problems in Engineering
 (2017)

---

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


---

## Body

## 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:(1)J=Mx-1,t0,x+1,tf+tf-t02∫-1+1Lxτ,uτ,tτ,t0,tfdτ,subject to the constraints(2)dxdτ=tf-t02fxτ,uτ,tτ,t0,tf,Cxτ,uτ,tτ,t0,tf≤0,Bx-1,t0,x-1,tf≤0,where the term x(τ)∈Rnx denotes the state and the term u(τ)∈Rnu denotes the control. In (1)-(2), the time domain τ∈[-1,+1] is transformed from the time domain t∈[t0,tf] by the following affine transformation [12]:(3)τ=2ttf-t0-tf+t0tf-t0,t=tτ=12τtf-t0+t0+tf,where the terms t0 and tf 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 [-1,1], so it will be useful to transform the time domain.


## 3. Adaptive Mesh Iteration Methodology

### 3.1. Numerical Discretization

The domainτ∈[-1,+1] is divided into K mesh subintervals Sk when using mesh refinement. Then we have (4)⋃k=1KSk=-1,+1,Sk=Tk-1,Tk.The mesh points have the property-1=T0<T1<⋯<Tk=+1. The state in the subintervals Sk is approximated by the Hermite interpolating polynomial withnth order [11] (5)xτ=a0+a1τ+⋯+anτn,τ∈-1,+1.Then the cost function is approximated by Gauss-Lobatto quadrature rule as follows:(6)J≈Mx11,t0,xNK,tf+tf-t02∑k=1Khi2∑j=1NkωjkLxjk,ujk,tτik,where the term ωj(k) is the same as ωj 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 Sk,   k=1,…,K, with Nk HLGL points has been solved. The ensuing mesh with Mk=Nk+1 HLGL points τ^1k,…,τ^Mkk, where τ^1k=τ1k=Tk-1 and τ^Mkk=Tk. Assume further that xτ^1k,…,xτ^Mkk are the values of the state approximation at τ^1k,…,τ^Mkk. We then have (7)x^kτ^jk=x^kτk-1+tf-t02∑l=1Mkhl2fx^kτ^jk,ukτ^jk,tτ^ik,t0,tf,j=1,…,Mk.The absolute error and the relative error approximations atτ^1(k),…,τ^Mkk of the state are defined, respectively, as(8)Eikτ^lk=x^ikτ^lk-xikτ^lk,eikτ^lk=Eikτ^lk1+maxj∈1,…,Nk+1,k∈1,…,K⁡xikτlk,l=1,…,Mk,i=1,…,nx.The maximum relative error inSk is then defined as (9)emaxk=maxi∈1,…,nx,l∈1,…,Mk+1eikτ^jk.


### 3.3. Nonsmooth Solution Location

If a mesh interval has met the accuracy tolerance, that is,emax(k)≤ε, 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 κkτ be the curvature of the ith component of the state in mesh interval k, as (10)κkτ=X¨ikτ1+X˙ikτ23/2.Letκ-k and κmaxk be the mean and maximum value of κkτ, respectively. Then define rk as the ratio of the maximum to the mean curvature:(11)rk=κmaxkκ-k.


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

Whenemaxk>ε and if rk≤rmax, where rmax 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 k. Let NkM and NkM+1 denote the number of collocation points in interval k at mesh M and M+1, respectively, where M is the mesh refinement iteration. The number of points NkM+1 at mesh M+1 is calculated by the formula (12)NkM+1=NkM+Pk,Pk=ceil⁡log10⁡emaxkε.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 Nmax 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 Nk(M+1)>Nmax (i.e., NkM+1 exceeds the maximum allowable polynomial degree), then the mesh interval Sk must be divided into equally spaced subintervals.


### 3.5. Generation of New Mesh Segment

Assumeemax(k)>ε and rk>rmax, and then the kth 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 Nmin collocation points. Third, the new number of mesh intervals, Bk, is given by the formula(13)Bk=ceil⁡Bulog10⁡emaxkε,where Bu 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 Bk, can be rewritten as (14)Bk=max2,ceil⁡Bulog10⁡emaxkε.


### 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, andrk<rmax, and then the number of collocation points should be decreased. The number of points NkM+1 at mesh M+1 is calculated by the formula (15)NkM+1=NkM-Pk′,Pk′=ceil⁡log10⁡emaxkε.


### 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 Section3.5 and then generally estimate the number of mesh interval points. If Nk+1≠Nk, the mesh interval Sk+1=Tk,Tk+1 and mesh interval Sk=Tk-1,Tk 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(16)J~≈J+ρ∑k=1K∑j=1Nkmax⁡0,CXjk,Ujk,tτik,where the term ρ represents penalty factor and Cj(k)=CXjk,Ujk,tτik. 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].(17)pjkε,ρ=0,Cjk<-ερn,ρn2εCjk+Cjk+ε2ρn,-ερn≤Cjk<0,Cjk+ε2ρn,Cjk>0.Obviously, the functionp is first-order differentiable with respect to C and has a good approximation to the max operation. When the smoothing factor ε→0, then we have pjkε,ρ=max⁡0,Cjk. The optimal index is converted into (18)J~≈J+ρ∑k=1K∑j=1Nkpjkε,ρ.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:(19)Hx,u,λ=L+λTf+pε,ρ.


### 3.9. Mesh Refinement Method

The schematic of adaptive mesh iteration method is shown in Figure1. The adaptive mesh iteration method is summarized as follows.Figure 1
Schematic of adaptive mesh iteration method.Step 1.
SetM=0 and supply initial mesh, S=⋃k=1KSk=-1,+1, where ⋂k=1KSk=⌀.Step 2.
Solve NLP on current meshS.Step 3.
Compute maximum relative erroremax(k) in Sk, k=1,…,K, if emax(k)≤ε for all k=1,…,K or M>Mmax, and then quit. Otherwise, proceed to Step 4.Step 4.
Ifemax(k)>ε, k=1,…,K, proceed to Step 5; otherwise, proceed to Step 6.Step 5.
Compute the ratio between the maximum and the mean curvaturerk in Sk, and if rk≤rmax, set the number of collocation points according to (15); else divide the interval Sk into Bk subintervals, where Bk is given by (14). Then proceed to Step 7.Step 6.
For the single mesh interval, reduce the number of collocation points. Merge the adjacent mesh intervals if they satisfy the merger conditions.Step 7.
SetM=→M+1, and return to Step 2.


## 3.1. Numerical Discretization

The domainτ∈[-1,+1] is divided into K mesh subintervals Sk when using mesh refinement. Then we have (4)⋃k=1KSk=-1,+1,Sk=Tk-1,Tk.The mesh points have the property-1=T0<T1<⋯<Tk=+1. The state in the subintervals Sk is approximated by the Hermite interpolating polynomial withnth order [11] (5)xτ=a0+a1τ+⋯+anτn,τ∈-1,+1.Then the cost function is approximated by Gauss-Lobatto quadrature rule as follows:(6)J≈Mx11,t0,xNK,tf+tf-t02∑k=1Khi2∑j=1NkωjkLxjk,ujk,tτik,where the term ωj(k) is the same as ωj 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 Sk,   k=1,…,K, with Nk HLGL points has been solved. The ensuing mesh with Mk=Nk+1 HLGL points τ^1k,…,τ^Mkk, where τ^1k=τ1k=Tk-1 and τ^Mkk=Tk. Assume further that xτ^1k,…,xτ^Mkk are the values of the state approximation at τ^1k,…,τ^Mkk. We then have (7)x^kτ^jk=x^kτk-1+tf-t02∑l=1Mkhl2fx^kτ^jk,ukτ^jk,tτ^ik,t0,tf,j=1,…,Mk.The absolute error and the relative error approximations atτ^1(k),…,τ^Mkk of the state are defined, respectively, as(8)Eikτ^lk=x^ikτ^lk-xikτ^lk,eikτ^lk=Eikτ^lk1+maxj∈1,…,Nk+1,k∈1,…,K⁡xikτlk,l=1,…,Mk,i=1,…,nx.The maximum relative error inSk is then defined as (9)emaxk=maxi∈1,…,nx,l∈1,…,Mk+1eikτ^jk.


## 3.3. Nonsmooth Solution Location

If a mesh interval has met the accuracy tolerance, that is,emax(k)≤ε, 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 κkτ be the curvature of the ith component of the state in mesh interval k, as (10)κkτ=X¨ikτ1+X˙ikτ23/2.Letκ-k and κmaxk be the mean and maximum value of κkτ, respectively. Then define rk as the ratio of the maximum to the mean curvature:(11)rk=κmaxkκ-k.


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

Whenemaxk>ε and if rk≤rmax, where rmax 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 k. Let NkM and NkM+1 denote the number of collocation points in interval k at mesh M and M+1, respectively, where M is the mesh refinement iteration. The number of points NkM+1 at mesh M+1 is calculated by the formula (12)NkM+1=NkM+Pk,Pk=ceil⁡log10⁡emaxkε.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 Nmax 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 Nk(M+1)>Nmax (i.e., NkM+1 exceeds the maximum allowable polynomial degree), then the mesh interval Sk must be divided into equally spaced subintervals.


## 3.5. Generation of New Mesh Segment

Assumeemax(k)>ε and rk>rmax, and then the kth 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 Nmin collocation points. Third, the new number of mesh intervals, Bk, is given by the formula(13)Bk=ceil⁡Bulog10⁡emaxkε,where Bu 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 Bk, can be rewritten as (14)Bk=max2,ceil⁡Bulog10⁡emaxkε.


## 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, andrk<rmax, and then the number of collocation points should be decreased. The number of points NkM+1 at mesh M+1 is calculated by the formula (15)NkM+1=NkM-Pk′,Pk′=ceil⁡log10⁡emaxkε.


## 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 Section3.5 and then generally estimate the number of mesh interval points. If Nk+1≠Nk, the mesh interval Sk+1=Tk,Tk+1 and mesh interval Sk=Tk-1,Tk 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(16)J~≈J+ρ∑k=1K∑j=1Nkmax⁡0,CXjk,Ujk,tτik,where the term ρ represents penalty factor and Cj(k)=CXjk,Ujk,tτik. 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].(17)pjkε,ρ=0,Cjk<-ερn,ρn2εCjk+Cjk+ε2ρn,-ερn≤Cjk<0,Cjk+ε2ρn,Cjk>0.Obviously, the functionp is first-order differentiable with respect to C and has a good approximation to the max operation. When the smoothing factor ε→0, then we have pjkε,ρ=max⁡0,Cjk. The optimal index is converted into (18)J~≈J+ρ∑k=1K∑j=1Nkpjkε,ρ.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:(19)Hx,u,λ=L+λTf+pε,ρ.


## 3.9. Mesh Refinement Method

The schematic of adaptive mesh iteration method is shown in Figure1. The adaptive mesh iteration method is summarized as follows.Figure 1
Schematic of adaptive mesh iteration method.Step 1.
SetM=0 and supply initial mesh, S=⋃k=1KSk=-1,+1, where ⋂k=1KSk=⌀.Step 2.
Solve NLP on current meshS.Step 3.
Compute maximum relative erroremax(k) in Sk, k=1,…,K, if emax(k)≤ε for all k=1,…,K or M>Mmax, and then quit. Otherwise, proceed to Step 4.Step 4.
Ifemax(k)>ε, k=1,…,K, proceed to Step 5; otherwise, proceed to Step 6.Step 5.
Compute the ratio between the maximum and the mean curvaturerk in Sk, and if rk≤rmax, set the number of collocation points according to (15); else divide the interval Sk into Bk subintervals, where Bk is given by (14). Then proceed to Step 7.Step 6.
For the single mesh interval, reduce the number of collocation points. Merge the adjacent mesh intervals if they satisfy the merger conditions.Step 7.
SetM=→M+1, and return to Step 2.


## 4. Simulation Analysis

The order and intervals of the method in [12] are fixed in each simulation, while those of the method described in this paper are variable. For the convenience of narration, the mesh refinement method in [12] is called FMRM (fixed mesh refinement method), and the method in Section 3 is called AMIM (adaptive mesh iteration method). The term M denotes the mesh refinement iteration, and M=0 means the mesh initialization, and the term N and term K denote the total collocation points and interval number, respectively. The number of collocation points within each interval of the two methods is at least 2. The maximum of all mesh interval allowable error values is ε, where ε=10-6. When the mesh is initialization, the whole mesh is divided into 10 intervals, and each interval has a number of 2 collocation points. The value of term rmax is 1.2, and the maximum number of collocation points with each interval is 12. The penalty factor ρ is 100 [15]. The simulation results of this paper were performed on a 3.4 GHz Intel Core i7 CPU computer and MATLAB Version R2013.Consider the following trajectory optimization problem taken from [16], and the state equations for a three-dimensional point mass vehicle which is commonly used in midcourse guidance systems of a Maneuverable Research Reentry Vehicle (MaRRV) are listed as follows:(20)r˙=vsin⁡γ,ϕ˙=vcos⁡γsin⁡χrcos⁡θ,θ˙=Vcos⁡γcos⁡χr,v˙=-Dm-gsin⁡γ,γ˙=1vLcos⁡β-g-v2rcos⁡γ,χ˙=Lsin⁡βmvcos⁡γ+vrcos⁡γsin⁡χtan⁡θ,where r is the radial position and the terms ϕ and θ are the latitude and longitude, respectively. The term v is the total velocity, and the terms γ and χ denote the flight path angle and azimuth angle, respectively. The symbol m represents the vehicle mass and g is the gravity (μ/r2 with μ being the gravitational parameter). The atmospheric density with altitude is assumed to be exponential model as(21)ρr=ρ0e-r-Re/H,where ρ0 and H are constants in Table 1.Table 1
Hypersonic boost-glide vehicle date and physical constants.
Quantity
Numerical values
Mass (kg)
1200
W
/
S (kg/m2)
420
R
e (m)
6.378
×
1
0
6
ρ
0 (kg/m3)
1.225
μ (m2/s2)
3.986
×
1
0
14
β (m)
7200The lift and drag are defined as(22)L=qdSCLα,D=qdSCDα,where (23)CLα=-0.034+0.93α,CDα=0.037-0.01α+0.736α2+0.936α3.The objective is to maximize the range, so the cost function is given by(24)J=min⁡-θf.Figures2(a), 2(b), and 2(c) show the final solutions of the optimal control problem. It can be seen from Figures 2(a), 2(b), and 2(c) that the solution to the trajectory is relatively smooth, though there seems to be a rapid change in the attack angle in Figure 2(c) near t=1000 s. It is noted that the control variable is limited to -5deg,15deg and that the solutions showed in Figure 2(c) satisfy that constraint. The results demonstrate the ability of approaches described in this paper for solving optimal control problems with inequality path constraints. In particular, it can be seen that the attack angle remains of 15 deg near t=1000 s. One might suppose that if the attack angle is not constrained, its value would likely be larger than 15 deg.Figure 2
(a) Altitude versus range. (b) Velocity versus time. (c) Attack angle versus time.
(a)
(b)
(c)Figure3 shows the collocation points distribution of the solutions obtained using AMIM. When M=0, the initial mesh is constructed. From M=0 to M=2, the mesh points are added largely, while the increment is relatively small from M=2 to the iteration process terminated. As Figures 4(a) and 4(b) show about the trajectories iteration process and attack angle iteration process, respectively, it can be seen that there is a great difference from mesh initialization to Mesh Iteration 1 and Mesh Iteration 2 of the trajectories; moreover, the gap between Mesh Iteration 2 and Mesh Iterations 3–4 is very small, and the situation of attack angle iteration process shown in Figure 4(b) is the same as the trajectories process, which means that the solutions are gradually converged on the actual solution with each mesh refinement iteration. It is noted that a higher number of mesh points are needed to ensure the convergence of the solutions, and that explains why the increments in Figure 3 from mesh initialization to Mesh Iteration 2 and from Mesh Iteration 2 to Mesh Iteration 4 are different.Figure 3
Mesh point history.Figure 4
(a) Trajectories iteration. (b) Attack angle iteration.
(a)
(b)A comparison of the implementation of the AMRM and FMRM with different higher-order and intervals solutions are given in Table2. Comparisons are made in terms of computation time, the number of total mesh points N and the number of mesh intervals K and the number of mesh refinement iterations M, and the cost function values. In each case, the initial guesses are randomized controls, with randomized state. A total of 100 samples are used to produce the results in this paper. The terminology AMRM (2, 11) refers to the AMRM where the number of mesh points within each interval can vary between 2 and 11; furthermore, the number of mesh intervals in AMRM can vary as well. All the simulations parameters are shown in black roman font in Table 2, and all the results are shown in black italics. As can be seen from the results listed in Table 2, the AMRM result in the smallest overall times compared with other cases. The reason is that the AMRM has the properties of reducing the unnecessary points and intervals, while FMRM (other cases) does not have this kind of property but only keeps the number of mesh points and intervals fixed until the simulation terminated. It is a fact that computation times mostly depend on the number of mesh refinement iterations and mesh size, while the growth of computation time for cases 4, 7, and 10 is due mostly to the increment number of mesh points and mesh intervals. Interestingly, the number of mesh intervals in case 1 is not set parameter but the result from the simulation, where the numbers of mesh intervals in cases 2–10 are set parameters. Case7 using 7 mesh points with 55 intervals gives a terminal cost function of −4163.4 km, which is the most optimality one in cases 2–9. The results show, as expected, that using larger mesh size results in improvements in accuracy at the expense of more runtime, due to the denser Jacobians. For FMRM, 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 [12]. 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 and the number of mesh intervals when solving a complicated problem. The simulations show that the AMEM can trade the number of mesh points within intervals with the number of mesh intervals and obtain an accurate solution with a relatively small mesh size.Table 2
Compare results using VOIand FOI methods.
Case
Methods
n
K
Mean times/s
N
M
Cost function
1
AMRM
(2, 11)
65
6.3
214
4
−4165.3
2
FMRM
5
50
19.6
201
8
−
4050.2
3
FMRM
5
60
33.8
241
7
−
4103.7
4
FMRM
5
65
41.3
261
6
−
4128.5
5
FMRM
7
40
25.2
241
7
−
4085.2
6
FMRM
7
45
38.1
271
6
−
4114.6
7
FMRM
7
55
47.8
331
5
−4163.4
8
FMRM
9
25
23.0
201
8
−
4016.1
9
FMRM
9
30
36.4
241
7
−
4073.9
10
FMRM
9
35
49.7
281
5
−
4122.8


## 5. Conclusions

This paper develops an adaptive mesh iteration method for solving optimal control problems. The method has the ability to trade the number of mesh points with number of mesh intervals compared to other mesh refinement methods. The mesh size is guided by a previously derived maximum relative error and the ratio of the maximum to the mean curvature of the mesh intervals. The number of mesh intervals is increased in regions where solution is nonsmooth, while the mesh points increased in regions where the solution is smooth. Furthermore, the mesh size can be decreased either by reducing the mesh points or by combining adjacent mesh intervals which share the same number of mesh points. The method applied successfully trajectory optimization of MaRRV from the open literature and compared with other mesh refinement methods. The simulation results demonstrate the superiority of the AMRM.



---

*Source: 2184658-2017-06-18.xml*