Abstract
To solve the multiobjective optimization problem on hypersonic glider vehicle trajectory design subjected to complex constraints, this paper proposes a multiobjective trajectory optimization method that combines the boundary intersection method and pseudospectral method. The multiobjective trajectory optimization problem (MTOP) is established based on the analysis of the feature of hypersonic glider vehicle trajectory. The MTOP is translated into a set of general optimization subproblems by using the boundary intersection method and pseudospectral method. The subproblems are solved by nonlinear programming algorithm. In this method, the solution that has been solved is employed as the initial guess for the next subproblem so that the time consumption of the entire multiobjective trajectory optimization problem shortens. The maximal range and minimal peak heat problem is solved by the proposed method. The numerical results demonstrate that the proposed method can obtain the Pareto front of the optimal trajectory, which can provide the reference for the trajectory design of hypersonic glider vehicle.
1. Introduction
Boost-glide vehicle has become a research hotspot because of its unique advantages, such as increasing range and improving ability of penetration. Trajectory optimization technology is one of the key aircraft design technologies. Reentry trajectory optimization usually needs to consider multiple performance indicators. Some of these objectives are conflicting, such as velocity and heating value [1, 2]. To solve these conflicts of multiobjective trajectory optimization problem, the following methods are usually adopted: the most important objective is chosen as single objective and other objectives are converted to constraints; multiobjective optimization problem is converted to single objective optimization problem using the weight-added sum method. But these methods can only get single solution instead of the optimal Pareto front. This paper mainly proposes a research on multiobjective trajectory optimization problem of hypersonic glide vehicle.
In recent years, there have been many studies on trajectory planning, such as pseudospectral method [3, 4] and support vector machine-based method [5, 6]. Because of the wide application of the pseudospectral method in solving hypersonic glider vehicle (HGV) trajectory optimization problem, this paper is based on pseudospectral method. With the wide application of multiobjective optimization, many new multiobjective optimization methods based on weighted method have appeared. One of the most representative methods is boundary intersection (BI) method. Das and Dennis [7] proposed the normal boundary intersection (NBI) method to solve the multiobjective optimization problem. After that, Messac et al. [8] proposed the normalized normal constraint (NNC) method. These two methods belong to boundary intersection method [9]. The basic theory of BI is as follows: the Pareto front of consequent multiobjective optimization problem (MOP) is part of the lower left boundary of objective set. The basic idea of the BI method is to look for the intersection points of the boundary and a set of straight lines and to approximate the Pareto front with intersection points. If this group of straight lines is evenly distributed, a uniform approximation of Pareto front will be gotten. The BI method converts MOP to a group of single objective optimization problems.
In this paper, the pseudospectral method is combined with boundary intersection method so that the multiobjective trajectory optimization problem can be converted to a set of single objective optimization subproblems, and then the subproblems are solved with pseudospectral method individually. In order to shorten the computing time, the solution of the subproblem is used as the initial value of the next subproblem. In this paper, the solving method of HGV multiobjective trajectory optimization problem has the following advantages:(1)Compared with the weight-added sum method, this method can avoid repeating design when engineering practice experience is lacking.(2)Compared with multiobjective trajectory optimization method based on evolutionary algorithm, this method can get sufficient accurate Pareto optimal solution with relatively less calculating quantity.
Content is arranged as follows: first, mathematical model of HGV trajectory optimization, typical objectives, and constraints are analyzed; then, multiobjective trajectory optimization algorithm based on pseudospectral method and the NBI method is proposed; finally, the multiobjective trajectory optimization problem on HGV maximal range-minimal peak heat flux is solved with the method in this paper and the Pareto front result is analyzed.
2. Glide Vehicle Trajectory Optimization Model
2.1. Reentry Dynamics Model
HGV is a lifting body with flat shape and it adopts bank-to-turn (BTT) swerve technique generally; therefore, flight sideslip angle can be set to zero. In this paper, the Earth rotation is ignored for simplification and three-degree-of-freedom motion model is established as (1). Location parameters of HGV are described by geocentric distance , longitude , and latitude . Velocity parameters of HGV are described by velocity magnitude , heading angle , and flight-path angle . Heading angle is the angle between velocity vector and the local horizontal. is positive when the velocity vector points to the water level upwards. Flight-path angle is the angle between velocity local level projection of velocity vector and north direction. is positive when the angle from north direction to the velocity vector is clockwise. These parameters are shown in Figure 1. The three-degree-of-freedom dynamics equation is described as follows [10]:
In the equations above, is the bank angle, that is, the angle between the lateral axis of an aircraft in flight and the horizontal. is aerodynamic drag force and is aerodynamic lift force. is the mass of HGV. is acceleration of gravity at current location. Hence,where is atmospheric density. is vehicle reference area. is drag coefficient and is lift coefficient. and are generally calculated as functions of angle of attack and Mach number Ma.
From the perspective of the optimal control, the state variables of vehicle reentry dynamics system are location parameters and speed parameters; control variables are the angle of attack and angle of bank . The purpose of trajectory optimization design is to find the optimal angle of attack and angle of bank to make the performance indicators minimum.
2.2. Trajectory Constraints
Constraints of flight process, terminal parameters, and control variable must be considered in trajectory design. Flight process constraints are the constraints where trajectory parameters must be satisfied in the flight process. HGV reentry is a complex flight process, so heat flow, dynamic pressure, overload, and maneuvering capability of the aircraft must be considered. These elements should not exceed the affordability of aircraft. Terminal parameter constraints are the conditions that should be satisfied in the end of the aircraft trajectory. Control constraints refer to the limit of attack angle and bank angle.
2.2.1. Heating Rate Constraint
Reentry trajectory design must consider the influence of aerodynamic thermal protection system (TPS). Material of TPS decides the limit of aircraft surface temperature and heat flux. To ensure the safety of aircraft to return, generally the stagnation heat flow restrictions along the trajectory of the flight are set aswhere is heating rate (its unit is ) and depends on the heat transfer coefficient of aircraft head shape. The stagnation heat flow is influenced by material of TPS. and are constants; because of the hypersonic reentry problems, they can be set as or and [10].
2.2.2. Aerodynamic Load Constraint
The impact of aerodynamic load on the internal structure should be considered with the aircraft reentry. The horizontal aerodynamic load constraint of HGV reentry is not serious, so reentry trajectory design mainly considers vertical aerodynamic load constraint:
2.2.3. Dynamic Pressure Constraint
Considering dynamic pressure effect on the spacecraft attitude control system and stability, dynamic pressure is an important parameter. Aerodynamic torque and pneumatic torque are directly related to dynamic pressure. In addition, the dynamic pressure affects the stability and execution efficiency of the aircraft aerodynamic control surfaces. Postural stability (especially the lateral stability) requires that the dynamic pressure does not exceed a certain range, so that the requirement on stability of flight control system will be satisfied:
2.2.4. Equilibrium Glide Constraint
Reentry vehicles should have sufficient capacity to meet the requirements of motorized guidance and control systems. When the trajectory height is too large, due to the thin air, the available lift is insufficient to balance gravity; once the aircraft is disturbed, vehicles cannot track the aircraft predetermined trajectory. Thus, it will affect the performance of guidance and even make the mission unable to be completed. To ensure controllable trajectory, the aircraft should keep the vertical upward force available to balance other forces:
2.2.5. Control Variables Constraints
During the flight, due to hardware limitations, the control variable such as attack angle and bank angle should not exceed the constraint value:
2.2.6. Terminal Constraints
Terminal constraint is determined by the mission. For example, when the requirement is to hit the target, the trajectory terminal location parameters should be consistent with the target. And when the requirement is security reentry, the trajectory terminal velocity will be requested. According to practical requirements, the reentry vehicle usually requires landing speed and landing trajectory inclination angle:
2.3. Performance Indicators
According to the design requirements, trajectory optimization chooses different performance indicators. When analyzing the aircraft performance, the maximum range, maximum speed, and other parameters are usually taken as performance indicators. For example, the maximum terminal trajectory range as optimization indicator is taken as follows:
And when the end position of trajectory is determined, optimization indicators such as the shortest path, the shortest time cost, or the minimum total amount of heat may need to be considered. For example, the minimum total amount of heat as optimization indicator is chosen as follows:
3. Multiobjective Trajectory Optimization Method
3.1. Multiobjective Trajectory Optimization Problem
The difference between multiobjective optimization problem and general trajectory optimization problem is the number of optimization objectives. These constraints are basically the same. With respect to single objective trajectory optimization, multiobjective trajectory optimization problem can be stated as follows [3].
The control function (and times ) is selected to make multiple performance functions minimize:
For state equation constraints,
For boundary condition constraints,
For path constraints,
is the state vector and boundary conditions are , . Single objective function is as follows:
Equations (11)–(15) give the mathematical model of multiobjective optimization trajectory. An effective method for solving this kind of problem is transforming it into a general multiobjective parametric optimization problem and then using gradient-based algorithms (such as sequential quadratic programming) to solve. In this paper, NBI and pseudospectral methods are used to solve the multiobjective HGV trajectory optimization problem.
3.2. Normal Boundary Intersection Method
To overcome the defect of weight-added sum method, Das and Dennis [7] proposed normal boundary intersection method. The basic idea of NBI is to find the intersection of the border and a group of straight lines to approximate Pareto front. In general, Pareto front of continuous MOP is a part of the lower left reachable border set. If this set of lines in some sense is uniformly distributed, it is available to get an approximate of Pareto front.
Calculation process of NBI method is as follows: first, find the single objective advantage; Then, create hyperplanes in the objective space through single objective optimal point (defined as Convex Hull of Individual Minima, CHIM) and build a group normal lines of this superplane. At last, solve lower left border intersection points of normal lines and reachable objectives set by single objective optimization methods to approximate Pareto front. Single objective optimization problem is converted by NBI method as
In (16), is weight vector, is configuration matrix, the th column element is , and is the quasi normal from CHIM to ideal spot. Constraint ensures that the objective function corresponding to solution is located in quasi-normal line. Optimization objective is the distance from the objective vector to CHIM.
Sometimes, for convenience, a set of lines from the ideal point can be used. In this case, the polymerization objectives are converted to find out the minimum distance spots from line. The spots are on a straight line to ideal spots in reachable space as follows:
The advantage of the NBI method is that it is not sensitive to the shape of Pareto front, and the resulting Pareto optimum is evenly distributed. But it needs to increase the number of equality constraints (consistent with the objective number), and more than one dominating solution may be found; then, the solution set should be filtered. Meanwhile, for some optimization problems with more than two goals, the solutions obtained by this method may not cover the entire Pareto frontier [11]. Figure 2 shows a schematic view of NBI method. Figure 2(a) is using the normal line of different points on CHIM as constraints, and Figure 2(b) is using the quasi-normal line from the ideal point as constraints.
(a) Translation
(b) Rotation
3.3. Pseudospectral Method
In recent years, the pseudospectral method has become a common tool for solving complex optimal control problems because of its high accuracy solution. The basic idea is to use the roots of orthogonal polynomial as discrete points and then separate the state variables and control variables of continuous optimal control problem and then use Lagrange interpolation of full region to approximate the state and control variables. Thus, the optimal control problem is converted to the nonlinear programming problem to be solved [4].
Currently, adaptive mesh method is widely used to select discrete points and automatically adjust mesh. adaptive pseudospectral method will change the mesh in accordance with the accuracy requirements in the iterative process [12]. relates to the number of segments and relates to the order of interpolation polynomial. When the error magnitude is similar on a segment, the number of points will be increased in this segment. When error magnitude of some isolated points is anomalous to other points, this segment will be refined to several segments. It is found that the adaptive mesh method leads to higher accuracy solutions with less computational effort and memory than what is required in a global pseudospectral method.
In this paper, hypersonic glide vehicle multiobjective trajectory optimization problem is converted to general multiobjective parameter optimization problem using adaptive pseudospectral method. The KKT (Karush-Kuhn-Tucker) condition of nonlinear programming problems transformed with pseudospectral method is consistent with the first-order optimality condition of discrete Hamilton boundary value problems under certain conditions [13]. Thereby, the covariates of the original optimal control problem will be obtained to make up the defect that the direct method cannot obtain covariate information. In order to apply the pseudospectral method and boundary intersection method for solving multiobjective trajectory optimization problem, pseudospectral method is modified on the basis of GPOPS Toolbox as follows:(1)In the objective function, Mayer and Lagrange type indexes are outputted, respectively, corresponding to each objective.(2)The sparse matrix template is modified.(3)The normal boundary intersection method is added.(4)The subproblem solutions obtained are taken as the initial value of next subproblem solving process.
Normal boundary intersection method obtains different optimization subproblems by changing the direction of the vector. When these problems are solved, Pareto optimum solutions are obtained. In general, the subproblems of the adjacent vector direction are relatively close to the objective function, so it is considered that the optimal parameters are also close relatively. Therefore, optimal solution of adjacent subproblem is used as initial value of next subproblem in order to improve efficiency of solving the optimization problem.
3.4. Process of Multiobjective Trajectory Optimization
The process of hypersonic glide vehicle multiobjective trajectory optimization based on normal boundary intersection method and pseudospectral method is shown as flow chart in Figure 3. Hence, we have the following:(1)Multiobjective optimization problem is constructed. Optimization objectives are chosen according to the task. Differential equations of reentry dynamics are established. The constraints of state variables and control variables needed to be satisfied are considered. Then, the mathematical model of multiobjective optimization problem is constructed. It should be noted that suboptimization objectives selected should be conflicting. For example, in trajectory optimization problem, minimize heat flux and maximize range or the minimum total absorption heat and maximum range, and so forth.(2)Multiobjective optimization problem is dispersed. Multiobjective mathematical model is dispersed by normal boundary intersection method. Weight vector is generated automatically and weight vector needs to meet and . Then, a series of optimal control subproblems are obtained and the solutions of these subproblems correspond to Pareto optimal solutions.(3)Single objective optimization subproblem is solved with pseudospectral method. The single objective optimization subproblems are solved with pseudospectral method and sequential quadratic programming method. Then, the subproblem solutions obtained are taken as the initial value of next subproblem solving process so that the computing time is shortened. Multiobjective Pareto front and corresponding solutions of optimal controls and states are obtained.
4. Maximum Range-Minimum Peak Heat Multiobjective Trajectory Optimization
4.1. Parameter Setting
No rotating spherical Earth model or exponent atmospheric model is adopted. Specific parameters are as follows.
4.1.1. Aircraft Parameters
Lift coefficient , drag coefficient , stagnation heat flux , , and . Other parameters in the above model are calculated as follows: sea-level atmospheric density kg/m3, reference height km, vehicle reference area m2, aircraft quality kg, and coefficients , , , , , , , , and .
4.1.2. Optimal Conditions
The initial states (reentry points parameters) are taken as follows: geocentric distance km ( is Earth radius), longitude , latitude , velocity m/s, heading angle , and flight-path angle .
Terminal constraint does not limit flight time but considers the height, speed, and trajectory angle requirements. Set parameters km, m/s, and .
Process constraints such as heat constraint, overload constraint, and dynamic pressure constraint are taken into account. Attack angle constraint and bank angle constraint are set as and .
4.1.3. Optimization Objectives
Maximum terminal range and minimum peak heat flux are taken as optimization objectives.
The first cost function is terminal range:
In the equation above, is initial longitude, is terminal longitude, is initial latitude, is terminal latitude, and is Earth radius.
Because , , and , in this paper, cost function is simplified as the projection of range on initial velocity direction:
The second cost function is aerodynamic heating on the vehicle wing leading edge [14]:
In the equations above, , , , , and .
For optimization objectives,
4.2. Analysis of Result
Pareto front of multiobjective trajectory problem is solved as in Figure 4. Minimum peak heat flux optimization solution and maximum terminal range optimization solution are shown in Table 1. Because computing resources of pseudospectral method closely depend on the number of iterations and the mesh grid nodes, the computing time of each subproblem has no linear relationship. Total computing time is adopted as calculation efficiency indicator. Calculation is carried out on Matlab platform on a computer with quad-core processor with processor core frequency of 3.30 GHz. Mesh tolerance in pseudospectral method is set as 10−3. Total computing time with iteration method in this paper is 391.67 minutes. Total computing time of all the subproblems solved severally by pseudospectral method is 463.50 minutes. Iteration method in this paper shortens the computing time.
From the results, Pareto front points are evenly distributed. In maximum terminal range solution, the value of terminal range is 6102 km, and its peak heat flux is 128.3 W/cm2. In minimum peak heat solution, the value of terminal range is 5597 km, and its peak heat flux is 67.92 W/cm2. The remaining objective points are between these two limits; this feature means that these two objectives appear to be in strong conflict. In addition, as can be seen from Pareto front in Figure 4, upper-right corner gradient of peak heat to range curve is larger. The curve of range to peak heat is similar. This phenomenon shows that the cost of edge improvement in these two optimization objectives should be big (the other cost function is sharply bad) and needs to be focused on in vehicle trajectory design.
Solution curves of minimum peak heat flow and maximum terminal range are shown in Figure 5. There are curves of height-time, velocity-time, latitude-time, peak heat flow-time, attack angle-time, and bank angle-time in Figure 5.
(a) Altitude-time curve
(b) Speed-time curve
(c) Range-time curve
(d) Peak heat-time curve
(e) Angle of attack-time curve
(f) Angle of bank-time curve
The simulation results show that the method based on boundary intersection method and pseudospectral method in this paper can effectively solve the multiobjective glide vehicle trajectory optimization problem, and this method can provide a useful reference to trajectory design.
5. Conclusion
In hypersonic glide vehicle trajectory design, multiple conflicting optimization objectives usually need to be considered using multiobjective optimization method. To solve this problem, this paper presents multiobjective hypersonic glide vehicle trajectory optimization method based on normal boundary intersection method and pseudospectral method. Multiobjective method optimization problem is converted to multiple single objective optimization problems with normal boundary intersection method, and then the optimal control problems are converted to parameter optimization problems with pseudospectral method, so that it can be solved with nonlinear programming algorithm. Hypersonic glide vehicle trajectory multiobjective optimization problem about maximum range and minimum peak heat is numerically simulated, and Pareto front solution is evenly distributed relatively. The varying range of trajectory conflicting objectives can be obtained from Pareto front, and it provides series of candidate solutions. This paper has reference value for hypersonic glide vehicle trajectory design.
Competing Interests
The authors declare that they have no competing interests.