Adaptive mesh refinement method for optimal control based on Hermite-Legendre-Gauss-Lobatto direct transcription

3.1. Numerical discretization and approximation

The domain $\tau \in [-1,+1]$ is divided into $K$ mesh subintervals $S_k$ when using mesh refinement method. Then we have:

The mesh points have the property$-1=T_0<T_1<\cdots <T_k=+1$. The state in the subintervals $S_k$ is approximated by the Hermite interpolating polynomial with nth order [9]:

For ensure the integration accuracy and interpolation accuracy, the collocation points and nodes are defined as LGL points ${\upsilon }_k(k=1,...,n)$ within each interval. Note that there is no distinction between collocation points ${\tau }_j$ and nodes ${\varsigma }_j$ in some PS method [2], whereas the collocation points are used to formulate the residual equations for the NLP, and the nodes are used to form the interpolating polynomial in this paper. Thus, the collocations points and nodes defined according to:

Note that the values of the states and states derivatives at the points ${\tau }_j$ determine the Hermite interpolating coefficients ${\mathbf{a}}_k(k=1,...,n)$:

where the interval length is defined by $h_i=(t_{i+1}-t_i)/2$.

It is noted that the left side of Eq. (9) determined by the number and location of the nodes. Let $\mathbf{a}=[{\mathbf{a}}_0,{\mathbf{a}}_1,...,{\mathbf{a}}_n{]}^T$, then $\mathbf{a}={\mathbf{A}}^{-1}\mathbf{b}$, where $\mathbf{A}$ is the matrix on the right side of Eq. (9), and $\mathbf{b}$ is the vector of the left side of Eq. (9). Then we have:

Differentiating $\mathbf{x}\left({\varsigma }_j\right)$ in Eq. (10) with respect to $\tau$, we obtain:

where:

The system constraints in per interval can be expressed as:

where:

The cost function is approximated by Gauss-Lobatto quadrature rule as:

where the term ${\omega }_j^{\left(k\right)}$ is the same as the ${\omega }_j$ in Ref. [16].

3.2. Approximation of solution error

The estimate method for the relative error is similar to the error estimate obtained for numerically solving a differential equation through using the modified Euler Runge-Kutta scheme. Suppose the NLP of Eqs. (1-4) on a mesh $S_k$, $k=1,...,K$ with $N_k$ HLGL points has been solved. The ensuing mesh with $M_k=N_k+1$ HLGL points $\left({\widehat{\tau }}_1^{\left(k\right)},...,{\widehat{\tau }}_{M_k}^{\left(k\right)}\right)$.

where:

Assume further that $\left(x\left({\widehat{\tau }}_1^{\left(k\right)}\right),...,x\left({\widehat{\tau }}_{M_k}^{\left(k\right)}\right)\right)$ are the values of the state approximation at $\left({\widehat{\tau }}_1^{\left(k\right)},...,{\widehat{\tau }}_{M_k}^{\left(k\right)}\right)$. We then have:

The absolute error and the relative error approximations at $\left({\widehat{\tau }}_1^{\left(k\right)},...,{\widehat{\tau }}_{M_k}^{\left(k\right)}\right)$ of the state are defined, respectively, as:

The maximum relative error in $S_k$ is then defined as:

3.3. Refining the mesh

If a mesh interval has met the accuracy tolerance, that is $e_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}\le \epsilon$, where $\epsilon$ is the desired relative tolerance, then mesh size is reduced by decreasing mesh polynomial degree or merging adjacent mesh interval, otherwise the mesh size need to be modified by increasing points or dividing the mesh interval into several subintervals. Let ${\kappa }^{\left(k\right)}\left(\tau \right)$ be the curvature of the $i$th component of the state in mesh interval $k$, as:

Let ${\stackrel{-}{\kappa }}^{\left(k\right)}$ and ${\kappa }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}$ be the mean and maximum value of ${\kappa }^{\left(k\right)}\left(\tau \right)$, respectively. Then define $r_k$ as the ratio of the maximum to the mean curvature:

When $e_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}>\epsilon$, and if $r_k<r_{\mathrm{m}\mathrm{a}\mathrm{x}}$, where $r_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is a user-defined parameter, the curvature is considered uniform in this interval mesh then the number of collocation points should be increased in interval $k$. Let $N_k^{\left(M\right)}$ and $N_k^{(M+1)}$ denote the number of collocation points in interval $k$ at mesh $M$ and $M+1$ respectively, where $M$ is the mesh refinement iteration number. The number of points $N_k^{(M+1)}$ at mesh $M+1$ is calculated by the equation:

It is noted in Eq. (21) that the ratio of the maximum to the error tolerance have a direct effect on polynomial degree in mesh interval. An upper limit $N_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is set for the maximum allowable polynomial degree to make sure that the number of collocation points does not grow an unreasonably large value. If $N_k^{(M+1)}>N_{\mathrm{m}\mathrm{a}\mathrm{x}}$ (i.e.$N_k^{(M+1)}$exceeds the maximum allowable polynomial degree), then the mesh interval $S_k$ must be divided into equally spaced subintervals.

3.4. Generation of new mesh segment

Assume $e_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}>\epsilon$ and $r_k>r_{\mathrm{m}\mathrm{a}\mathrm{x}}$, then the $k$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 $N_{\mathrm{m}\mathrm{i}\mathrm{n}}$ collocation points. Third, the new number of mesh intervals $B_k$, is given by the equation:

where $B_u$ 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 $B_k$, can be rewritten as:

3.5. 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 if $r_k<r_{\mathrm{m}\mathrm{a}\mathrm{x}}$, then the number of collocation points should be decreased. The number of points $N_k^{(M+1)}$ at mesh $M+1$ is calculated by the equation:

3.6. 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, then generally estimate the number of mesh interval points. If $N_{k+1}\ne N_k$, the mesh intervals $S_{k+1}=[T_k,T_{k+1}]$ and $S_k=[T_{k-1},T_k]$ cannot be merged because highest polynomial order 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 estimations of the two grid intervals are not more than $\epsilon$.

(3) The relative error of the new mesh interval after the merger is not larger than $\epsilon$.

3.7. Mesh refinement method

The schematic of adaptive mesh refinement method is shown in Fig. 1. The adaptive mesh refinement method is summarized as follows.

Step 1: Set $M=$ 0 and supply initial mesh, $S=\bigcup _{k=1}^K{S}_k=[-1,+1]$, where $\bigcap _{k=1}^K{S}_k=\mathrm{\varnothing }$.

Step 2: Solve NLP on current mesh $S$.

Step 3: Compute maximum relative error $e_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}$ in $S_k$, $k=1,...,K$, if $e_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}\le \epsilon$ for all $k=1,...,K$ or $M>M_{\mathrm{m}\mathrm{a}\mathrm{x}}$, then quit. Otherwise, proceed to Step 4.

Step 4: If $e_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}>\epsilon$, $k=1,...,K$, proceed to Step 5, otherwise, proceed to Step 6.

Step 5: Compute the ratio between the maximum and the mean curvature $r_k$ in $S_k$, if $r_k\le r_{\mathrm{m}\mathrm{a}\mathrm{x}}$, set the number of collocation points increase by $P_k$, else divide the interval $S_k$ into $B_k$ subintervals, where $B_k$ is given by Eq. (23). Then proceed to Step 7.

Step 6: For the single mesh interval, decrease the number of collocation points. Merge the adjacent mesh interval if they satisfy the conditions of the merger.

Step 7: Set $M\overrightarrow{=}M+1$, and return to Step 2.

4.1. Example 1

Consider the following Bang-Bang optimal control problem from Ref. [18] to illustrate the effectiveness of the method in this paper. Firstly, the method is able to accurately solve the optimal control problem. Secondly, the mesh refinement with elimination of unnecessary mesh points is able to improve the algorithm performance.

Minimize the cost function:

The analytic solution for this optimal control problem is given as:

Fig. 2(a) and Fig. 2(b) and Fig. 2(c) show the exact solutions of the optimal control problem.

Fig. 2a) x(t) vs. t, b) u(t) vs. t, c) λ(t) vs. t

a)

b)

c)

Fig. 3(a) and Fig. 3(b) show the collocation point’s distribution of the solutions obtained using FOI method and VOI method respectively. According to the exact solutions of Eq. (29), we can know that near time $t=tc$, the state variables and control variables are rather changeable. Fig. 4(a) and Fig. 4(b) show the collocation points are mainly located near $t=tc$ instead of located at both ends of the solution, which is because more number of collocation points are needed to capture the changes at $t=tc$ of the solution. When $M\ge \text{2}$, the mesh point in Fig. 3(b) does not increase with each mesh iteration, but decreases. This is because the VOI method has the properties of reducing the interval number and the number of mesh points, and the FOI method does not have this kind of property, and mesh points in Fig. 3(a) is fixed. The iteration program is terminated at $M=$ 5 when the accuracy tolerance is satisfied, and the number of mesh points using the VOI method is 60.

Fig. 3a) FOI mesh point history, b) VOI mesh point history

a)

b)

Next, we analyze the approximation ability of solution obtained by the VOI method to the exact solution. Fig. 4(a) and Fig. 4(b) show the state on each mesh refinement iteration alongside the analytic solution using the VOI method. In addition, it is seen from Fig. 4(a) show that the resulting solution gradually converge to the exact solution with each mesh refinement iteration. Fig. 4(b) shows the states near $t=tc$ on each mesh refinement iteration, it is apparent that the difference between Mesh Iteration 1 and Analytic Solution is great, and the Mesh Iteration 2 is much closer to the analytic than Mesh Iteration 1, moreover, the gap between Mesh Iteration 2 and Mesh Iteration 3-5 is very small, which mean that the solutions are gradually converged on the analytic solution with each mesh refinement iteration. The number of mesh points near $t=tc$ is also increasing with the continuous refinement of the mesh, because there are larger state changes near $t=tc$, thus the precise requirement of the solution is satisfied with more mesh points.

Fig. 4a) x(t) vs. t, b) x(t) vs. t

a)

b)

A comparison of the implementation of the VOI method and FOI method with different higher-order and intervals solutions are given in Table 1. 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 iteration $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 VOI (2, 12) refers to the VOI method where the number of mesh points within each interval can vary between 2 and 12, furthermore, the number of mesh intervals in VOI method can vary as well. All the simulations parameters are shown in Table 1, and all the results are shown in blue. As it can be seen from the results listed in Table 1, the VOI method result in the smallest overall times compared with other cases. The reason is that the VOI method has the properties of reducing the unnecessary points and intervals, while FOI method (other cases) doesn’t have this kind of property but only keeps the number of mesh points and intervals fixed until the simulation terminated. It is a fact computation times are mostly depended on the number of mesh refinement iterations and mesh size, while the growth of computation time for case 4, 7, 10 are due mostly to the increase in 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 number of mesh intervals in case 2-10 are set parameters. The case 10 using 9 mesh points with 15 intervals gives a terminal cost function of 0.4572, which is the most optimality one in all cases. The results show, as expected, that using larger mesh size results in improvements in accuracy at the expense of increases in runtime, due to the denser Jacobians. For FOI method, better accuracy can be achieved by increasing mesh points for smooth problems, whereas increasing the number of intervals to achieve better accuracy for nonsmooth problems [16]. 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 VOI method 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.

Next, we analyze the convergence of mesh refinement. Table 2 shows the estimated maximum relative errors and exact relative errors for each mesh refinement iteration by using the VOI method (case 1) and FOI method (case 2, 7, 10). First, it can be seen from the Table 2 that the relative error on final mesh is quite small at $\approx$ 10^-7 for the state. The consistency in the exact relative error and the relative error approximation demonstrates the accuracy of the estimate derived in section 3.2.

4.2. Example 2

Consider the following trajectory optimization problem taken from Ref. [19] of maximizing the downrange of a Maneuverable Research Re-entry Vehicle (MaRRV). Minimize the cost function:

The state equations for hypersonic vehicle which is commonly used in midcourse guidance systems are listed as followings:

It is noted that the model (physical model and wing-body vehicle model) are taken from Ref. [19]. The initial conditions and terminal constraints are listed in Table 3. A typical solution of this problem is shown in Fig. 5(a)-(d) by using the VOI (2, 12) method.

It is seen that the solution to this example is relative smooth, especially the control variable (attack angle) is slowly changing in Fig. 5(d), meaning that it is easier to apply to engineering. As a result, it is possible to achieve an accurate solution with a relatively small number of mesh points when compared with FOI method. Table 3 shows that the terminal constraints are satisfied with error of less than 0.5 % in all parameters. A comparison of the implementation of the VOI method and FOI method solutions are given in Table 4.

Fig. 5a) Altitude vs. downrange, b) velocity vs. time, c) flight path angle vs. time, d) attack angle vs. time

a)

b)

c)

d)

It is seen that the VOI (2, 11) and FOI (5, 55) methods [shown in Table 4] are almost the most computationally efficient. FOI (9, 45) method takes the most time (with 185.7 seconds), while FOI (7, 55) produce the most optimal solution with the smallest cost function. As a result, for the vast majority of the solution, the largest decrease in error is achieved by using more mesh intervals and more mesh points in each mesh interval. The reason for more computation time needed to obtain the desired solution is that the number of mesh intervals and mesh points in each mesh interval are fixed when uses FOI methods. The simulations show that the proposed method can efficiently obtain a reliable, accurate solution for re-entry trajectory optimization problem of MaRRV.