An improved adaptive hp mesh refinement method in solving optimal control problems

In order to improve the efficiency of the adaptive hp pseudospectral method in solving the optimal control problems, an improved adaptive hp pseudospectral method is proposed in this paper. This method mainly improves the mesh refinement process of the traditional hp pseudospectral method to reduce the number of mesh iterations, so as to achieve the purpose of shortening the solving time. The focus of this method is twofold. Firstly, by judging the position of the sudden change points of control curve, the mesh is refined in advance around them. Secondly, the non‐uniform segmentation iteration is carried out by using the curvature of the system state curve as the criterion. These two points together improve the solving efficiency of the proposed method. Through the convergence analysis and simulation verification, it is shown that the proposed method can effectively solve the optimal control problem, and it is also proved that the proposed method has higher solving efficiency than the traditional adaptive hp pseudospectral method.

In solving the optimal control problem, many scholars have actively studied the related theory and application of pseudospectral method. For three different types of collocation points (LG, LGR and LGL), Garg et al. [4][5][6] study the theoretical foundation and algorithm framework of pseudospectral method, which lays an important foundation for the research of this algorithm. In reference 4, a unified framework of pseudospectral method for three collocation strategies is proposed, and the characteristics of the algorithm under each strategy are deeply analyzed. For the trajectory optimization problem in finite-horizon and infinite-horizon, the accuracy of solution of pseudospectral method based on LGR collocation points is studied in reference 5. For the infinite-horizon optimization control problem, Garg et al. 6 study the influence of several different time domain transformation methods on the solution results of pseudospectral method based on the two collocation points of LG and LGR. In order to improve the solving efficiency of NLP, Patterson et al. 7 propose a method to calculate the first and second derivatives of the NLP function. For the optimal control problem with terminal constraints, Fahroo et al. 10 reveal the superiority of the pseudospectral method based on Legendre polynomial in dealing with this type of problem. Using the covector mapping principle, Gong et al. 11,12 prove that the discrete solution obtained by the pseudospectral method can converge to a continuous optimized solution.
However, if only a fixed number of collocation points and a single mesh interval are used, the pseudospectral method will be difficult to obtain a satisfactory solution for complex optimal control problems. In order to solve this problem, many scholars have put forward two ideas. One is the p method, 13 which improves the solving precision by increasing the collocation points. The other is the h method, [14][15][16] which increases the number of mesh intervals by continuously dividing the time domain to ensure the solving precision. According to these two ideas, many scholars have made meaningful exploration. If there is a smooth analytical solution to the optimal control problem, the p method can get good results. Hussaini et al. [17][18][19] has shown that in this case, by increasing the collocation points, a good approximation of the analytical solution and the exponential convergence rate can be achieved. However, if the solution of the optimal control problem is non-smooth, such as there are discontinuities or sudden changes, in this case, even if the number of collocation points is increased to increase the order of the interpolation polynomial, it is difficult to obtain a satisfactory solution. In response to this difficulty, h method is proposed, where original time domain is divided and analytic solution is approximated by establishing different mesh intervals. Relevant research has proved that this method can effectively make up for the shortcomings of p method. [14][15][16] Considering that the h method and the p method have their own advantages for different types of questions, some scholars have appropriately combined the two to further improve the solving precision and efficiency. [20][21][22][23] Patterson et al. 20 propose the ph method, which is actually a mesh refinement strategy using p method first and then h method. In other words, when using the ph method for mesh refinement, the order of the interpolation polynomial will first be increased to improve the solving precision. Then, if the order reaches the specified maximum limit, the mesh will be segmented instead until a satisfactory solution is obtained. Unlike ph method, a mesh refinement method called hp method is proposed in references 21-23, which adaptively decide whether to use h method or p method to refine the mesh through certain criteria. Darby et al. 21 take the consistency of approximation error as the criterion, that is, if the approximation error of each collocation point is relatively consistent, the p method is used to refine mesh, otherwise, the h method is used. Different from this, Liu et al. 22,23 take the smoothness of the solution as the criterion, that is, if the obtained solution is smooth, the h method is used for refinement, otherwise, the p method is used. Due to the efficient solving ability and the flexibility of the hp algorithm, its related theories and applications have been widely studied by many scholars. [24][25][26][27][28] Li et al. 29,30 apply the adaptive hp method to solve trajectory optimization problem of missile, which obtain satisfactory solving precision. Considering obstacles avoidance and overload minimization, Xia et al. 31 propose an hp-pseudospectral sequential convex programming method to achieve trajectory optimization of unmanned aerial vehicle. Based on adaptive hp method, Xu et al. 32 complete the optimal motion planning of the space manipulator. Liu et al. 33 use hp-adaptive method realize path-tracking of vehicle. Yazdani et al. 34 analyze the feasibility of hp-adaptive Radau pseudospectral method. By improving the mesh refinement process, some effective variants of the hp method are proposed in references 35-39. In addition, several numerical methods [40][41][42][43] are proposed, which also provide important ideas for solving various types of optimal control problem.
Although the aforementioned works 21-39 on the adaptive hp method has made great progress, there still are several issues that need further consideration.
1. When there are sudden change points in the optimal control curve (such as bang-bang control, 44,45 switch control 46 ), in order to obtain a satisfactory solution, some existing adaptive hp methods often configure a large number of collocation points near the sudden change point, or segments the mesh frequently, which largely increases the dimensions of the NLP, therefore, affects the computational efficiency. 2. Most of the existing hp methods divide the mesh evenly according to the criterion, which is inefficient and may lead to frequent segmentation in the process of mesh refinement. This process will inevitably increase the number of iterations and prolong the computing time. 3. When hp method is used to solve optimal control problem with Bolza form, most existing works do not give the strict proof of convergence of this algorithm. Although the convergence of hp method is analyzed in reference 52, it is only discussed for solving unconstrained optimal control problems.
In view of the aforementioned consideration, this paper proposes an improved adaptive hp mesh refinement method, which mainly has three contributions.
1. By accurately locating the sudden change points of the control curve, the mesh points are added in the vicinity of these points. With the targeted mesh refinement in advance, the number of iterations of the grid is reduced and the efficiency of the algorithm is improved. 2. In order to further improve the efficiency of mesh refinement, a mesh segmentation method based on the cumulative sum of state curvature values is designed. This non-uniform segmentation method (unlike the uniform mesh segmentation method in references [21][22][23][24] can ensure that the locations that need more refinement can be quickly allocated to more mesh points. 3. For the Bolza optimization control problem with equality and inequality constraints, we extend the results in reference 52 and further prove the convergence of the improved adaptive hp algorithm proposed in this paper.
The remaining part of the paper is organized as follows. The Section 2 is the problem formulation, which describes the optimal control problem to be solved in this paper. In the Section 3, based on LGR collocation method, the original optimal control problem is transformed into a nonlinear programming problem. The Section 4 is the description of the improved adaptive hp mesh refinement method, and the core algorithm of this paper is introduced in detail. In Section 5, we prove the convergence of the proposed algorithm. The simulation example in the Section 6 verifies the effectiveness and superiority of the method in this paper. In the Section 7, the work of this paper is summarized.

PROBLEM FORMULATION
Without loss of generality, we define an optimal control problem with Bolza form as in (1), that is, to find the state-control function (x, u) ∈ R n x × R n u such that Where, x ∈ R n x is the state variable, u ∈ R n u is the control variable, and x and u are defined as row vectors. t 0 and t f represent initial time and terminal time, respectively, and . E(⋅) and F(⋅) are scalar functions, which represent the terminal type and integral type performance indicators, respectively. f(⋅) ∈ R n f represents dynamic constraints, e(⋅) ∈ R n e represents boundary conditions, h(⋅) ∈ R n h represents inequality path constraints. Define f, e and h as row vectors. All these nonlinear functions (E, F, f, e, h) are continuously differentiable in their domains, and their gradients are Lipschitz continuous.
According to Pontryagin minimum principle, 47 we can reduce the first-order necessary condition of problem B to the following boundary-value problem. 48,49 Where H(t) represents the Lagrangian of the Hamiltonian.
H is the Hamiltonian function.
E represents the Lagrangian form of the endpoint constraint.
In addition, the Lagrangian multipliers , , and in Equations (3)-(5) are row vectors of n h dimension, n f dimension and n e dimension, respectively. In order to solve the problem B, here we adopt the idea of mesh division in the hp pseudospectral method, select K + 1 mesh points . Let x (k) (t) and u (k) (t) denote the state and the control on the mesh interval S k , respectively. Therefore, problem B Equation (1) can be equivalent to Since the state must be continuous in time interval , the following conditions should be satisfied at each internal mesh point.

LEGENDRE-GAUSS-RADAU COLLOCATION
In order to transform problem B into nonlinear programming problem (NLP), here we use the LGR collocation method [4][5][6] to discretize Equation (6). First of all, we need to introduce time-domain variation to transform time t of each grid interval . Introduce a new time variable and make the transformation as shown below.
Let h k denote the length of each mesh interval S k , that is, h k = T k − T k−1 . Combined with Equations (6) ∼ (8), problem B can be rewritten as Then, select N k LGR collocation points (k) is the root of polynomial P N k −1 + P N k (P N k is Legendre polynomial of order N k ). In order to discretize Equation (9), we use LGR collocation points (k) 1 , … , (k) N k and noncollocated point (k) N k +1 = +1 as sampling points on each mesh interval S k to achieve Lagrange interpolation on x (k) ( ). The interpolation polynomial based on these N k + 1 points is expressed as follows Thus, the continuous state x (k) ( ) on mesh interval S k can be approximated as Here, . Taking the derivative of Equation (11) at the collocation . Then, by interpolating , its integral can be approximated as Where quadrature weight (k) j is Combining Equations (10)- (14), we can transform the optimal control problem (Equation (9)) into an NLP, which is given For problem B N , we can use the NLP solver SNOPT 9 to solve it. However, the optimal solution solved directly may not be able to meet the accuracy requirements of the original problem B. Therefore, we need to continuously adjust the position of mesh points [T 0 , T 1 , … , T K ] (h method) and the number of LGR collocation points in each mesh interval (P method), and then iteratively solve problem B N until the optimal solution meets the accuracy of the original problem B. In this paper, we propose an improved adaptive hp mesh refinement method.

IMPROVED ADAPTIVE HP MESH REFINEMENT METHOD
For some optimal control problems, the sudden changes of the control and the excessive curvature of the state curve will bring difficulties to the mesh refinement process of the adaptive hp method. However, in most of the existing studies, [20][21][22][23][24] the traditional h method is used to segment the mesh evenly, which leads to frequent mesh refinement, increasing the number of iterations and reducing the solving efficiency. In response to these problems, this paper has improved the traditional adaptive hp mesh refinement method and formed a new algorithm with two core points. One is to find out the interval containing control sudden changes by evaluating the derivativeU (k) i at different allocation points, and then refine the interval emphatically. Second, different from the uniform segmentation method adopted by the traditional adaptive hp method, this paper determines the location of the mesh points according to the curvature of the state curve, that is, where the curvature is large, the mesh points are dense, otherwise they are sparse.
This section mainly introduces the improved adaptive hp mesh refinement method proposed in this paper. Section 4.1 is the evaluation of the solution error. In Section 4.2, we give a detailed description of the core algorithm of mesh refinement proposed in this paper.

Evaluation of solution error
After solving the solution of problem B N using the NLP solver, it is also necessary to judge whether the solution can achieve satisfactory accuracy, that is, whether the error between X (k) ( ) and x (k) ( ) in Equation (11) can be within a tolerance range . Since we cannot know the analytical solution x (k) ( ) of problem B in advance, here we take the following method to evaluate the solution error in each mesh interval. First, add an LGR collocation point in each grid interval. For the mesh interval S k (k = 1, 2, … , K), denote the new Γ k collocation points as Then, based on the numerical solution and the corresponding sampling points , the state ) and the control U (k) can be approximately expressed as can be approximately expressed as Finally, combining the X (k) and obtained above, integrate again to obtain the stateX using the integral form of the discrete state equation.
For the new LGR collocation points , construct a Γ k × (Γ k + 1) differentiation matrixD (k) , and the Similar to the derivation of the discrete state equation in Equation (15), the discrete state equation at collocation points can be expressed as follows , At the same time, an Γ k × Γ k integration matrixÎ is the i th column element of the differentiation matrixD (k) . Therefore, combined with Equation (19), we can obtain the integral form of the discrete state equation as followŝ .
So far, we can use to approximate the error between the numerical solution X (k) ( ) and the analytical solution x (k) ( ) of problem B, which is expressed as follows , respectively.

Refining the mesh
we use e (k) max calculated in Equation (21) to represent the error of numerical solution . Set is the accuracy tolerance range. If e (k) max ≤ holds for all mesh interval S k , k = 1, … , K, then the numerical solution obtained is considered to be satisfactory and the iterative solving is stopped. On the contrary, it is necessary to refine the original mesh and continue the iterative solving. The detailed mesh refinement method is as follows.

4.2.1
The mesh refinement in advance based on control sudden changes For many optimal control problems, there are some sudden changes on the control curve, such as bang-bang control, 44,45 switching control 46 and so on. When hp pseudospectral method is used to solve problem B N , the sudden changes of control F I G U R E 1 Optimization results of system (22) curve will greatly increase the number of iterations and affect the efficiency of solution. For example, consider the time optimal control problem of the following second-order integrator model.
In order to find out the optimal control u * (t) to minimize the cost function J, the adaptive hp method in reference 21 is used to solve Equation (22). The results are as Figure 1.
According to the optimal control theory, this example is bang-bang control, and its analytical solution is By comparing Equation (23) with Figure 1A, when the adaptive hp method 21 is used to solve this example, the optimization result is in good agreement with the analytical solution. On the other hand, in the process of solving the example, in order to achieve the desired approximation accuracy, dense mesh points and collocations are needed at the sudden change point (t = 2.22s) of the control curve, which greatly increases the number of iterations and reduces the efficiency of the solving (see Figure 1B). In order to effectively solve the above problems, in the improved adaptive hp method, we refine the mesh in advance based on the sudden change points of control. The basic idea is to accurately locate the sudden change points, and then refine the mesh around the points first to improve the efficiency.
Firstly, in order to locate the sudden change points, take time derivative along u (k) Where t i = ( i + 1) h k ∕2 + T k−1 , Equation (24) can be further written aṡ Then, according to the derivative of the control and the set threshold, the interval of the control changing too fast is selected. Suppose that the derivative of control at each collocation point in interval S k is as shown in the Figure 2.
In the Figure 2, is the set threshold. When the derivative of the control exceeds (such as t and Figure 2), we determine that there may be control sudden change points in this area, so we need to refine the mesh in advance.
Finally, mesh points are added in the area where the control changes sharply, so as to refine the grid in this area first. Taking Figure 2 as an example, we will add new mesh points t −1 , t +1 , t −2 , t +2 to further refine the grid of . In order to describe the subsequent mesh refinement process, the new mesh points are denoted as [ , and the new mesh intervals are denoted as

4.2.2
Criterion of mesh classification refinement based on state curvature When using hp method to solve the optimal control problem, the mesh interval should be classified first and then refined. [20][21][22][23] The state in some mesh intervals change dramatically, so the segmented method is suitable for this kind of interval, that is, h method. While for those mesh intervals where the state change smoothly, it is suitable to use the method of adding collocation points to refine, that is, p method. In the improved adaptive hp method, we take the curvature of state as the criterion of classification. The set of collocation points in the new mesh interval S ′ k is denoted as Ω t k , and the state at each collocation point is expressed as Ω x(t k ) . Based on Ω t k and Ω x(t k ) , N points with equal time intervals are selected in the interval Lagrange interpolation, and the coordinates of the interpolation points are recorded as . The curvature of each point can be calculated by using the interpolation point coordinates Combined with the curvature, we set the criterion of classification refinement as If the curvature threshold parameter is set to and Cri (k) > is true, it can be considered that the state in the mesh interval changes too fast and need to be refined by segmentation. On the contrary, if Cri (k) ≤ , it is necessary to refine the mesh by adding collocation points.

Increasing collocation points based on maximum approximation error
If it is determined in section 4.2.2 that the number of collocations needs to be further increase in mesh interval S ′ k , here we will increase the number of collocations according to the maximum approximation error e (k) max in Equation (21) in the following way Where, N k is the original allocation number of interval S ′ k , N ′ k is the increased allocation number, and is the set accuracy tolerance range.
[⋅] represents the least integer function.
The rationality of this method is that if the solution error of the mesh interval is large, more points need to be added to improve the accuracy of the solution. If the error is small, a small number of collocations is added to ensure the solving speed.

4.2.4
The segmentation method based on cumulative sum of curvature When the traditional adaptive hp method performs grid segmentation, the mesh interval that needs to be segmented is usually equally divided into several sub-intervals. This approach does not maximize the efficiency of segmentation, that is, it does not allocate more mesh points to those locations that need more refinement. In order to further improve the efficiency of mesh refinement, we design a segmentation method based on the cumulative sum of curvature, that is, assigning more mesh points to the position where the curvature of the state curve is larger. Firstly, the number of subintervals to be added is estimated. The principle is that if the approximation error of the solution in the mesh interval is large, the interval needs to be divided into more subintervals. On the contrary, the number of subintervals is less. According to this principle, the number of sub-intervals divided in mesh interval S ′ k is determined by the following formula Where max(•, * ) is the maximum of • and * . is the set product factor. Then, according to the curvature of the state curve in mesh interval S ′ k , the position of the new mesh point is determined. It is known from Equation (29) that interval S ′ k needs to be divided into M (k) subintervals, that is, M (k) − 1 mesh points need to be added between T ′ k−1 and T ′ k . In order to determine the location of these mesh points, we need to refer to the curvature (k) i at each interpolation point in Equation (26) to ensure that more dense mesh points are configured at the position where the state changes rapidly. Firstly, the continuous curvature curve is taken as an example to analyze, assuming that the curvature changes with time as shown in the figure below. According to the above mesh segmentation principle, if the original interval needs to be divided into four sub-intervals (as shown in Figure 3 respectively, and S represents the total area bounded by the curvature curve and the time axis in the whole interval . Referring to the above analysis, for the discrete curvature, the mesh point configuration rule of the improved adaptive hp method in this paper is to select M (k) − 1 points from the N interpolation points (k) i in the interval S ′ k , and record them as m (k) i , i = 1, … , M (k) − 1, so that the sum of the curvatures at the interpolation points in each sub is approximately equal. According to this rule, the algorithm for selecting new mesh points is as follows

Algorithm procedure
So far, the key steps of the improved adaptive hp method proposed in this paper have been described, and the process of its core algorithm is summarized as follows 1.
Step 1: Select the initial mesh points [T 0 , T 1 , … , T K ], and set the initial collocation points Step 2: Based on the current mesh points and collocation points, use the NLP solver to solve the problem B N , and get the numerical solution Step 3: Use the method in section 4.1 to approximate the error between the numerical solution and the analytical solution, and compare the maximum error e (k) max in each mesh interval S k with the accuracy tolerance to determine whether the original grid needs further refinement. If e (k) max ≤ is established in all mesh intervals, the existing numerical solution can be considered as a satisfactory solution, and skip to Step 4. Otherwise, skip to Step 5. 4.
Step 4: Complete the iteration of the solution, and take the current numerical solution as the final output. 5.
Step 5: The existing numerical solutions do not meet the accuracy requirements, so the mesh needs to be further refined: i Determine the location of the sudden change points based on the change rate of control, add new mesh points around them, and refine the original mesh points into [ , so as to continue to mesh classification Refinement. If Cri (k) in Equation (27) is greater than the threshold , go to step iii. Otherwise, go to step iv. iii According to the method in section 4.2.4, the mesh interval S ′ k is segmented based on curvature accumulation. iv Based on the maximum approximation error e (k) max in Equation (21), as described in 4.2.3, add a certain number of collocation points in the mesh interval S ′ k .

F I G U R E 4
Algorithm procedure 1.
Step 6: The refined mesh points and collocation points in step 5 are brought back to step 2 for a new round of solving.
The procedure of this method is shown in Figure 4.

CONVERGENCE ANALYSIS
In order to prove the convergence of the algorithm, the following assumptions are made for problem B (Equation (1)).

Assumption 1. For l ≥ 3, in the domain
there is an optimal solution ( x * , u * , * , * , * ) for problem B (Equation [1]) so that the Pontryagin minimum condition (Equation [2]) holds. At the same time, there is > 0 and open set Ω ⊂ R n x × R n u such that B (x * (t), u * (t)) ⊂ Ω holds for any t ∈ x * ( t f )) .

Assumption 3. For any
Next, analyze the convergence of the improved adaptive hp algorithm proposed in this paper, i.e. to prove that for the NLP B N (Equation [15]) there is a minimum point and Lagrangian multiplier vectors (k) i , (k) i and (k) such that converges to the optimal solution Theorem 1. If Assumption 1∼Assumption 3 hold, then for problem B N (Equation [15]), there is a minimum point ( In addition, k = 1, … , K. N k is the number of LGR points in the mesh interval S k , h k represents the interval length of S k , C k is a constant greater than 0, and the definition of l k is similar to the definition of l in Assumption 1. Proof. First, consider the first-order necessary conditions for the minimum point of problem B N (Equation [15]).
Introduce Lagrangian multiplier vectors (k) i ∈ R n f , m (k) i ∈ R n h and (k) ∈ R n e (all three are assumed to be row vectors), and construct the following Lagrangian function . At the same time (k) = ( . Calculate the gradient of Lagrangian function L in Equation (32) to respectively. The first-order necessary conditions for the minimum point are as follows (1) ∇ X (K) (2) ∇ X (1) Combining Equations (37)(38)(39)(40), rewrite the first-order necessary conditions (33) ∼ (36) into the following form (2) ∇ X (1) ) The proof of Theorem 1 requires the following Proposition. 51 ▪ Proposition 1. If 51 Assumption 1∼Assumption 3 are established. For * ∈ and r > 0, let the mapping Γ ∶  → be continuously Frechet differentiable within B r ( * ) . Where is a Banach space and is a normed linear space. If the following conditions are true.
Then there is unique ∈ B r ( * ) such that Γ( ) = 0, and at the same time there is . At the same time, combining Equations (15) and From reference 52, we know that the conditions C1 and C2 in Proposition 1 are established. Next, consider the size of According to Assumption 1, * satisfies Equations (1) and (2), so for Γ ( * ) defined above, there is Next, in order to continue to examine the size of , first introduce the following lemma. 52

Lemma 1. If the Lebesgue constant N is defined as
Here, the definition of Lagrange interpolation polynomial l j ( ) is shown in Equation (16). Then there is N = O(log N).

Lemma 2. If the Lebesgue constant N−1 is defined as
The definition of Lagrange interpolation polynomial l j ( ) is shown in Equation (17). Then there is, ) .

Lemma 3. For any y
is defined as shown in Equation (16), then Where, the definition of Lebesgue constant N is shown in Equation (56).

Lemma 5. According to the definition in Equations
That is, if p is a polynomial of degree N k − 1 and the ith element of vector P ∈ R N k is P i = p i = p ( i ), then there is The proof of Lemma 1 ∼ 5 can be found in reference 52.
In order to analyze the size of ‖ ‖ ‖ Γ 3 ( * ) ‖ ‖ ‖∞ , combining with Equation (2), we can know that in Equation (48), * (k) . Therefore, Equation (47) can be simplified to From Lemma 5, Where, the superscript N k − 1 denotes that * (k) is a polynomial of order N k − 1 obtained by Lagrange interpolation. Combining Equations (61) and (62), there is Let According to Lemmas 2, 3, and 4, we can get Since Equations (64) and (65) hold for any i and j, there is And for ‖ ‖ ‖ Γ 2 ( * ) ‖ ‖ ‖∞ , according to the similar analysis method, we can also get Finally, in order to analyze the size of is a differentiation matrix, Equation (51) can be written as Where, the superscript N k denotes that x * (k) is a polynomial of order N k obtained by Lagrange interpolation. Let Combined with Lemma 2, 3, and 4, we can get

SIMULATION ANALYSIS
We verify the effectiveness of the improved adaptive hp algorithm through several groups of simulation examples, [53][54][55] and illustrate the advantages of the proposed method by setting up the simulation comparison with the methods in the related research. In this part, all simulations rely on Lenovo desktop computers, whose processor is Intel(R) Core(TM) i5-4570, CPU clock speed is 3.20 GHz, memory is 4.00 G. And the simulation software platform is MATLAB 2014a. In addition, the algorithm proposed in this paper is developed based on the open source hp pseudospectral optimization software GPOPS, 56 and the NLP solver used is SNOPT. 9 Example 1. Soft lunar landing problem.
The optimal control problem for soft landing on the moon was proposed in reference 53, and it is restated here as follows: Find the optimal control u such that Where h is the height from the lunar surface and is the falling velocity. g represents the gravitational constant of the lunar surface, where g = 1.6. According to the optimal control theory, problem (74) is a bang-bang control problem, and the analytical solution of the optimal control u * (t) can be expressed as Combining with the initial conditions set by Equation (74), we can calculate that the switching time t * s = 1.3117 s and the terminal time t * f = 4.2394 s. At the same time, we use the improved adaptive hp algorithm proposed in this paper to solve the problem (74), and the relevant parameters in the algorithm are selected as follows: = 10 −6 , = 1, = 1.2, = 3. The initial mesh points are taken as [−1, 1], and the number of initial LGR collocation points is taken as 6. In addition, under the same tolerance range of mesh accuracy, the hp method in reference 23 and the ph-(N min , N max ) method in reference 20 are used to solve the problem and compare with the method in this paper. The notation ph-(N min , N max ) here represents the range of the number of collocation points of the ph method from N min to N max . The simulation results are shown in Figure 5.
As shown in Figure 5A-C, all three methods can accurately solve the problem under the same tolerance range of mesh accuracy, while the accuracy of the method in this paper is slightly higher than that of the other two methods. On the other hand, at the control sudden change point t = 1.3117 s, the hp method and the ph- (3,14) method need to greatly increase the number of collocation points and the number of mesh intervals to achieve an effective solving to the original problem, which greatly increases the number of mesh iterations and sacrifices the efficiency of solving (see Figure 5E,F). As shown in Figure 5D, the improved adaptive hp method proposed in this paper can accurately locate the position of the control sudden change point, so as to achieve a higher solution accuracy with a smaller number of collocation points and iterations. In addition, from the comparison results in Table 1, it can be seen that the proposed method needs less mesh iterations (only 2 iterations), fewer collocation points and mesh intervals, which makes the method more efficient (CPU time is 0.14 s, significantly less than the other two methods).

Example 2. Hyper-sensitive problem.
Consider the Hyper-sensitive optimal control problem mentioned in reference 54 and restate it as follows: The terminal time t f is fixed, which is taken as t f = 5000 here. The improved adaptive hp method proposed in this paper is used to solve the problem (76), and the initial mesh points are set to [−1, 1], and the initial number of LGR points is set to 6. The relevant parameters in the algorithm are set as: = 10 −6 , = 1, = 1.2, = 3. Similarly, the hp method 23 and the ph-(N min , N max ) method 20 are compared with the method in this paper. The simulation results are shown in Figure 6. It can be seen from the simulation results that the results obtained by the three methods are basically the same, and a large number of LGR points are allocated in the initial and final stages of x(t) curve (see Figure 6A,C,E). However, the mesh refinement process of these three methods is different: Although the ph- (3,14) method can continuously improve the solution accuracy by increasing the number of LGR points and mesh intervals, it requires a large number of LGR points and the dense grid to deal with this problem (see Figure 6F), which affects the efficiency of the solving. The reason is that when the method is performing mesh refinement, only when the number of collocation points increases to the upper limit N max will the mesh-segmented refinement begin, so that a large number of LGR points are inevitably needed at positions where the state changes rapidly. While hp method can determine whether the grid needs to be segmented refined by calculating the curvature of the state, which makes the mesh refinement process more targeted to a certain extent, reduces the number of LGR points required (see Figure 6D), and improves the efficiency of solving. The algorithm proposed in this paper improves the method of evenly dividing the grid in hp method, so that the position where the faster the state changes, the more collocations can be allocated, which further strengthens the pertinence of hp method for mesh refinement, and significantly reduces the number of mesh iterations and the total number of collocations required (see Figure 6B). Table 2 also shows the same conclusion. As shown in  Table 2, compared with the other two methods, the improved adaptive hp method proposed in this paper requires significantly fewer collocation points and mesh iterations when solving the problem, and the solving efficiency is also higher (CPU time is 0.14 s).

Example 3. Dynamic soaring problem.
Considering that the models of the optimal problems involved in the first two simulation examples are simple systems, in order to further test the performance of the algorithm proposed in this paper to solve the optimal problems of complex systems, the following control problems are introduced here. 55 Find the optimal control (C L , ) such that = (x(0), y(0), h(0), V(0), (0), Ψ(0) − 2 ) Here, lift L and drag D are expressed as The definitions and values of related variables and parameters involved in Equations (77) and (78) are detailed in reference 21.
In order to find out the optimal control (C L , ) satisfying Equation (77), the method proposed in this paper is used to solve the problem. The initial mesh point is set as [−1, 1], and the initial number of LGR points is set as 7. The relevant parameters of the algorithm in this paper are set as: = 10 −6 , = 1, = 1.2, = 3. At the same time, the hp method 23 and the ph-(N min , N max ) method 20 are used to solve the problem, which is compared with the method in this paper. The simulation results are shown in Figure 7. It can be seen from the simulation results that the method proposed in this paper can effectively solve the optimal control problem of a complex system of the form (77), and the other two methods have basically the same solution results (see Figure 7A,C). However, compared with the other two methods, the method in this paper has obvious advantages in mesh refinement, which requires fewer collocation points and mesh iterations (see Figure 7E,F). This is also reflected in the statistical results in Table 3. As shown in Table 3, when the method proposed in this paper is used to solve the problem, the grid only needs to be iterated 3 times, and the total number of collocation points required is significantly less than the other two methods, which has higher solving efficiency (CPU time is 4.24 s).

CONCLUSION
In order to improve the solving efficiency of the optimal problem, this paper proposes an improved adaptive hp mesh refinement method. In this method, the sudden change points of the control curve are accurately located, and then the mesh is refined first near these points, which greatly improves the efficiency of the algorithm. At the same time, a segmentation method based on the cumulative sum of curvature is designed to ensure that more mesh points are allocated at the position with larger curvature of the state curve. This method significantly reduces the number of mesh iterations and shortens the calculation time. Simulation results also show that the number of mesh iterations, the number of required collocation points, and the CPU time of the proposed algorithm are less than other methods. In fact, more mesh points also need to be configured in the position where the numerical solution approximation error is large to improve the solving precision. Therefore, in the future work, the insertion position of mesh points can be more reasonably determined by comprehensively considering the curvature value of state variables and the numerical solution approximation error, so as to further improve the solving efficiency of the algorithm.