A trajectory optimization method for the jump flight phase of a high-mach number aircraft

Abstract

The application discloses a high-Mach-number aircraft jump flight section trajectory optimization method and relates to the field of aircraft trajectory optimization, and comprises the following steps: S1, a dynamics equation and a kinematics equation are established, and a simplified model of a mass center motion equation set is obtained; S2, constraint conditions and an optimization target function are determined, and a trajectory optimization problem is described as a standard dynamic optimal control problem; S3, an AGA algorithm is used to globally optimize control variables, and a quasi-optimal solution of a globally optimal solution is obtained, so that the problem that a GPM-SQP algorithm is abnormally sensitive to an initial value when solving a nonlinear programming problem is solved; and S4, the quasi-optimal solution is used as an initial value of the GPM-SQP algorithm, and then the GPM-SQP algorithm is used to optimize the control variables in a guess time period, so that a globally optimal solution or a final solution close to the globally optimal solution is obtained based on the advantages of high calculation accuracy and fast convergence speed of the GPM-SQP algorithm. The method comprehensively utilizes the strong global search capability of the AGA algorithm and the fast convergence speed of the GPM-SQP algorithm, and thus an ideal trajectory meeting the constraint conditions and having the longest distance is easily obtained.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of aircraft trajectory optimization, in particular to a high Mach number aircraft jump flight segment trajectory optimization method. BACKGROUND

[0002] High Mach number aircraft has high strategic value such as high-altitude reconnaissance, high-speed penetration, long-range precision strike, etc., and has become a hot spot in the field of aerospace research. The "jump flight segment" is a trajectory of the high Mach number aircraft after reentry into the atmosphere, which slows down and descends until reaching the predetermined target. In order to achieve the precise strike target of the high Mach number aircraft and complete the flight mission, the jump flight segment trajectory optimization is the most important, which is based on the aircraft dynamics model, and a trajectory that meets various constraint conditions and has the optimal performance index is planned by solving the optimal time sequence of the control variables and the motion state.

[0003] There are the following difficulties in solving the jump flight segment trajectory optimization problem: first, the aircraft motion model is nonlinear and strongly coupled, and there is a certain difficulty in model establishment and trajectory optimization; second, the flight environment is complex and changeable, and the flight trajectory is disturbed by many factors such as model uncertainty and disturbance factors, which puts certain requirements on the robustness of the trajectory optimization algorithm; third, in the jump flight process, not only the constraints such as dynamic pressure, overload, heat flux rate need to be considered, but also the multiple no-fly zones formed by natural environment and other factors need to be avoided, and in addition, the strict restrictions on the terminal state, which poses a certain challenge to the study of trajectory optimization problem.

[0004] According to whether the performance index is directly optimized, the trajectory optimization method can be divided into indirect method and direct method. The essence of the indirect method is based on the classical variational method or Pontryagin minimum principle, which converts the original dynamic optimal control problem into a Hamiltonian two-point boundary value problem (HBVP) for solving, but for the high Mach number aircraft with complex dynamic characteristics, the convergence domain of the boundary value problem is small, and the initial value estimation of the adjoint variable is difficult. Unlike the indirect method, the direct method discretizes and parameterizes the variables in the original dynamic optimal control problem, converts the dynamic optimal control problem into a nonlinear programming problem (NLP), and further solves it by combining with the numerical optimization method, without the need to derive the first-order necessary condition, and is widely used. SUMMARY

[0005] The purpose of the present application is to design a high Mach number aircraft jump flight segment trajectory optimization method to solve the above problems.

[0006] The present application achieves the above-mentioned purpose by the following technical solutions:

[0007] A high-Mach-number aircraft jump flight section trajectory optimization method, comprising:

[0008] S1, ignoring the earth rotation and curvature influence, assuming the earth as a plane, respectively establishing dynamic equation and kinematic equation in the trajectory coordinate system and the ground coordinate system, obtaining 3 degrees of freedom high-Mach-number aircraft mass center motion equation group simplified model;

[0009] S2, determining the constraint condition suffered in the jump flight section flight process and the target function to be optimized, describing the trajectory optimization problem as a standard dynamic optimal control problem;

[0010] S3, based on the converted standard dynamic optimal control problem, first, the AGA algorithm is used to globally optimize the control quantity by taking advantage of the strong global search capability, and the quasi-optimal solution close to the global optimal solution is obtained;

[0011] S4, the quasi-optimal solution is used as the initial value of the GPM-SQP algorithm, and then the GPM-SQP algorithm is used to optimize the control quantity in the guess time period, and the global optimal solution or the final solution close to the global optimal solution is obtained.

[0012] The beneficial effects of the present application are that: the method comprehensively utilizes the strong global search capability of the adaptive genetic algorithm (AGA) and the fast convergence speed of the GPM-SQP algorithm, and the ideal trajectory meeting the constraint condition and the longest distance is easily obtained. Secondly, the trajectory optimization process is completed offline, the GPM-AGA-SQP hybrid algorithm proposed by the present application is based on the high-Mach-number aircraft as the research background, and is proposed for the deficiency that the GPM-SQP algorithm is sensitive to the initial guess value when solving the high-Mach-number aircraft jump flight section trajectory optimization problem under the condition of complex multiple constraints, but the algorithm also solves the problems of slow convergence speed and weak local optimization capability of the AGA algorithm when approaching the global optimal solution, and combines the advantages of strong robustness of the AGA algorithm, so that the GPM-AGA-SQP hybrid algorithm can still complete the trajectory optimization task under the condition of model uncertainty and disturbance interference. Finally, the GPM-AGA-SQP hybrid algorithm is essentially an optimization method, which is not only suitable for high-Mach-number aircraft jump flight section trajectory optimization, but also suitable for general flight trajectory optimization problems. BRIEF DESCRIPTION OF DRAWINGS

[0013] Figure 1 is the drag coefficient;

[0014] Figure 2 is the lift coefficient;

[0015] Figure 3 is the lift-drag ratio;

[0016] Figure 4 is the flowchart of the genetic algorithm of the present application;

[0017] Figure 5 is the "quasi-optimal solution" algorithm flowchart of the global optimal solution obtained by the AGA algorithm of the present application;

[0018] Figure 6 is the flowchart of the present application "quasi-optimal solution" as the initial solution of the GPM-SQP algorithm, quickly finding the global optimal solution;

[0019] Figure 7 is the best fitness of the adaptive genetic algorithm when simulating the trajectory optimization method of the high-Mach number aircraft in the jump flight section of the present application;

[0020] Figure 8 is the curve of the angle of attack angle changing with time after the adaptive genetic algorithm completes the optimization control quantity in the simulation;

[0021] Figure 9 is the curve of the speed tilt angle changing with time after the adaptive genetic algorithm completes the optimization control quantity in the simulation;

[0022] Figure 10 is the speed change curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0023] Figure 11 is the trajectory inclination angle change curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0024] Figure 12 is the trajectory deflection angle change curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0025] Figure 13 is the longitudinal range change curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0026] Figure 14 is the height change curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0027] Figure 15 is the lateral range change curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0028] Figure 16 is the angle of attack change curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0029] Figure 17 is the angle of attack change rate curve in the trajectory optimization problem with the farthest longitudinal range as the target in the simulation;

[0030] Figure 18Speed tilt angle variation curve simulated in trajectory optimization problem with the goal of maximum longitudinal range

[0031] Figure 19 Altitude-range variation curve simulated in trajectory optimization problem with the goal of maximum longitudinal range

[0032] Figure 20 Heat flow rate variation curve simulated in trajectory optimization problem with the goal of maximum longitudinal range

[0033] Figure 21 Dynamic pressure variation curve simulated in trajectory optimization problem with the goal of maximum longitudinal range

[0034] Figure 22 G-Force variation curve simulated in trajectory optimization problem with the goal of maximum longitudinal range

[0035] Figure 23 Flight path plot simulated in trajectory optimization problem with the goal of maximum longitudinal range

[0036] Figure 24 Three-dimensional plot simulated in trajectory optimization problem with the goal of maximum longitudinal range

[0037] Figure 25 Speed variation curve simulated in trajectory optimization problem with the goal of maximum lateral range

[0038] Figure 26 Trajectory inclination angle variation curve simulated in trajectory optimization problem with the goal of maximum lateral range

[0039] Figure 27 Trajectory bank angle variation curve simulated in trajectory optimization problem with the goal of maximum lateral range

[0040] Figure 28 Range variation curve simulated in trajectory optimization problem with the goal of maximum lateral range

[0041] Figure 29 Altitude variation curve simulated in trajectory optimization problem with the goal of maximum lateral range

[0042] Figure 30 Lateral range variation curve simulated in trajectory optimization problem with the goal of maximum lateral range

[0043] Figure 31 Angle of attack variation curve simulated in trajectory optimization problem with the goal of maximum lateral range

[0044] Figure 32The simulated attack angle rate curve in the trajectory optimization problem with the target of the maximum lateral range;

[0045] Figure 33 The simulated speed bank angle curve in the trajectory optimization problem with the target of the maximum lateral range;

[0046] Figure 34 The simulated altitude-range change curve in the trajectory optimization problem with the target of the maximum lateral range;

[0047] Figure 35 The simulated heat flow rate change curve in the trajectory optimization problem with the target of the maximum lateral range;

[0048] Figure 36 The simulated dynamic pressure change curve in the trajectory optimization problem with the target of the maximum lateral range;

[0049] Figure 37 The simulated overload change curve in the trajectory optimization problem with the target of the maximum lateral range;

[0050] Figure 38 The simulated flight path chart of avoiding the forbidden area in the trajectory optimization problem with the target of the maximum lateral range;

[0051] Figure 39 The simulated three-dimensional solid figure of avoiding the forbidden area in the trajectory optimization problem with the target of the maximum lateral range;

[0052] Figure 40 The performance index comparison curve of three algorithms (performance index: maximum range);

[0053] Figure 41 The performance index comparison curve of three algorithms (performance index: maximum lateral range);

[0054] Figure 42 The control variable attack angle comparison curve of three algorithms (performance index: maximum range);

[0055] Figure 43 The control variable speed bank angle comparison curve of three algorithms (performance index: maximum range);

[0056] Figure 44 The control variable attack angle comparison curve of three algorithms (performance index: maximum lateral range);

[0057] Figure 45 The control variable speed bank angle comparison curve of three algorithms (performance index: maximum lateral range);

[0058] Figure 46 The flight path chart of avoiding the forbidden area comparison curve of three algorithms (performance index: maximum range);

[0059] Figure 47 Three algorithm evasion no-fly zone overhead view comparison curve (performance index: the farthest range);

[0060] Figure 48 Monte Carlo simulation flight trajectory (performance index: the farthest longitudinal range);

[0061] Figure 49 Monte Carlo simulation flight trajectory (performance index: the farthest range). DETAILED DESCRIPTION

[0062] In order to make the objects, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are some but not all of the embodiments of the present application. The components of the embodiments of the present application described and shown in the drawings herein can be arranged and designed in various different configurations.

[0063] Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope of the claimed present application, but only represents selected embodiments of the present application. All other embodiments obtained by a person of ordinary skill in the art without creative labor based on the embodiments in the present application belong to the scope of protection of the present application.

[0064] The specific embodiments of the present application will be described in detail below with reference to the drawings.

[0065] A trajectory optimization method for a high-Mach number aircraft in a jump flight section, comprising:

[0066] S1, ignoring the earth rotation and curvature effect, assuming the earth as a plane, respectively establishing the dynamic equation and kinematic equation in the trajectory coordinate system and the ground coordinate system, obtaining a 3-DOF high-Mach number aircraft mass center motion equation group simplified model:

[0067]

[0068] In the formula: V(t) is the aircraft speed, θ(t) and ψ v (t) are the trajectory inclination angle and trajectory deflection angle respectively; y(t) is the height, x(t) and z(t) are the longitudinal range and lateral range respectively; m(t) is the aircraft mass, g is the gravity acceleration, α(t) is the attack angle, γ v (t) is the speed inclination angle. α B is the equilibrium attack angle, β B is the equilibrium side slip angle; Y B and Z B are α B , β BThe corresponding equilibrium lift and equilibrium lateral force, where X is the drag, are calculated using the following expressions:

[0069]

[0070] In the formula: C x C y For the drag coefficient and lift coefficient, q = 0.5 * ρ * V 2 , where ρ is the air density at the altitude of the high Mach number aircraft, and S is the reference area of ​​the aircraft.

[0071] It can be seen from equation (2) that the aerodynamic coefficient (including the drag coefficient C) x and lift coefficient C y The accuracy of the calculation directly affects the accuracy of aerodynamic forces (including drag and lift). This invention takes HTV-like aircraft as the research object and obtains aerodynamic coefficients under different flight conditions through numerical simulation calculations and evaluations. Figures 1-3 As shown, the reference length and reference area in the calculation of the lift coefficient and drag coefficient are taken as 1.0m and 1.0m, respectively. 2 .

[0072] S2. Determine the constraints and objective function to be optimized during the hop flight segment, describing the trajectory optimization problem as a standard dynamic optimal control problem. The constraints include self-control constraints, process constraints during the hop flight segment, and terminal constraints upon reaching the destination. The objective function is selected as the longest possible range (longest longitudinal range / longest lateral range), as follows:

[0073] Control constraints include angle of attack α and velocity tilt angle γ. v To ensure stable flight of the aircraft, the control variables are angle of attack α and velocity roll angle γ. v It needs to be limited to a certain range:

[0074]

[0075] Considering practical engineering applications, the rate of change of angle of attack is typically... Certain restrictions are also necessary, namely:

[0076]

[0077] At this point, the rate of change of angle of attack can be... It is treated as a virtual control variable, and the original control variable angle of attack α is treated as a state variable.

[0078] Process constraints include heat flow rate constraints, dynamic pressure constraints, overload constraints, and no-fly zone constraints. High Mach number aircraft experience significant aerodynamic heating during skip flight; therefore, heat flow rate constraints are necessary to ensure the structural safety of the aircraft.

[0079]

[0080] In the process of skip flight, the dynamic pressure is very large. In order to prevent the structural failure of the aircraft, the dynamic pressure needs to be constrained:

[0081] q = 0.5 * p * V 2 ≤ q dmax (6)

[0082] In order to ensure the normal operation of the related equipment of the aircraft, the overload of the aircraft needs to be limited within a certain range:

[0083]

[0084] The no-fly zone is a region formed due to natural environmental factors, which is not allowed to fly into. The no-fly zone of the present application selects a typical regular shape of a cylinder, the center point of which is (x i ,y i ,z i ), and the radius is r i . The specific expression is as follows:

[0085]

[0086] When the high Mach number aircraft skips, the trajectory thereof cannot enter the interior of the no-fly zone, that is, the shortest distance between the aircraft and the no-fly zone must be ensured to be positive. Let L(x(t), t) be the shortest distance from the aircraft to the no-fly zone, and let e n be an infinitesimal normal quantity. Then the no-fly zone constraint can be expressed as:

[0087] L(x(t), t) > e n

[0088] In the formula, K Q is a constant value related to the high Mach number aircraft, is the maximum allowable heat flow rate, q dmax is the maximum allowable dynamic pressure, and n Lmax is the maximum allowable overload.

[0089] The terminal constraints include the height, speed and trajectory inclination angle constraints at the terminal time. In order to meet the requirements of the mid- terminal guidance handover, the height, speed and trajectory inclination angle at the terminal time of the skip flight segment need to be set as constraints:

[0090]

[0091] In the formula, y f , V f and q f are the height, speed and trajectory inclination angle at the terminal time of the skip flight, respectively. y df , Vdf , θ df are the altitude, velocity and flight-path angle of the terminal time of the jump flight respectively; ε h , ε v , ε θ are the error upper bounds of the three constraints respectively;

[0092] The objective function refers to a function that makes a certain performance index optimal by adjusting the control variable on the basis of satisfying the above constraints. According to the specific task selection, the objective function can be generally divided into maximum range (maximum longitudinal range / maximum lateral range), shortest flight time, maximum velocity at the terminal time, etc. In the present application, the farthest range (the farthest longitudinal range / the farthest lateral range) is taken as the performance index, and the objective function is:

[0093] min J = -L df (10)

[0094] wherein L df refers to the range (longitudinal range / lateral range) of the flight vehicle.

[0095] Based on the above analysis, the trajectory optimization problem of the high-Mach number flight vehicle in the jump flight section can be described as: under the conditions of satisfying the mass center motion equation set (1), the control constraint equation (3)-(4), the process constraint equation (5)-(8), the terminal constraint equation (9), finding the optimal time sequence of the control variable attack angle α and velocity tilt angle γ v so that the objective function equation (10) is minimum. More generally, the optimal control problem can be described as follows: finding the control variable u(t) so that the Bolza type performance index (i.e. the combination of Mayer type and Lagrange type) is minimum

[0096]

[0097] and satisfying the mass center motion equation set constraint

[0098]

[0099] boundary condition constraint

[0100] φ(x(t0),t0,x(t f ),t f ) = 0 (13)

[0101] and inequality constraint

[0102] C(x(t),u(t),t)≤0(14)

[0103] wherein t0 is the starting time, t f is the terminal time.

[0104] S3, a quasi-optimal solution of the trajectory optimization problem (dynamic optimal control problem) is obtained by using an adaptive genetic algorithm (AGA) to solve the defect of sensitivity of the GPM-SQP algorithm to initial values (terminal time and initial and final time guess values of control variables). Genetic algorithm (GA) is a random global search optimization method. Its basic principle is to simulate the natural selection mechanism of "survival of the fittest" in the biological world. Through selection, crossover and mutation operations, the optimal traits are constantly evolved, and finally the chromosomes (individuals) that meet the optimization goal are generated. GA algorithm has the advantages of high parallelism, randomness and adaptive global optimization probability search, and has solved many complex engineering problems (such as multi-objective problem solving), and has been widely used in scientific research and engineering practice. The specific algorithm process is shown in Figure 4 The adaptive genetic algorithm (AGA) can be improved in several aspects such as fitness function, crossover probability and mutation probability based on the GA algorithm. Considering that the selection of the fitness function in the GA algorithm is directly related to the convergence speed and whether the optimal solution can be found, the AGA algorithm designed in this patent mainly improves the fitness function as follows:

[0105] 1) Adaptive adjustment of fitness function

[0106] In order to increase population diversity and further quickly and accurately approach the global optimal solution, the fitness function is designed to be self-adaptive adjustment with the increase of the number of population iterations. That is, in the initial stage of the algorithm, the fitness values of individuals in the population are quite different. In order to prevent the loss of individuals with poor fitness in the initial stage, the fitness function can be changed to retain them. When the algorithm approaches convergence, the fitness values of individuals in the population are small. In order to speed up the convergence speed, the fitness function can be adjusted to increase the difference between the fitness values of individuals.

[0107] 2) Logarithmic transformation of fitness function

[0108] Considering that the performance index and the terminal constraint condition have large difference in parameter magnitude, and cannot be directly used as the fitness function, it needs to be transformed, such as linear transformation, power law transformation and exponential transformation. This patent proposes an adaptive function based on the inverse function form of exponential transformation (i.e. logarithmic transformation), which can better reflect the advantages and disadvantages of individuals.

[0109] 3) Introducing multi-objective stratified planning method (Lexicographically Stratified Programming, LSP)

[0110] It is difficult to handle the fitness function that both satisfies the terminal constraints and optimizes the performance index at the terminal time. A multi-objective hierarchical planning method (LSP) can be introduced, i.e. multiple objectives are divided into multiple priority levels according to their importance. In the trajectory optimization of the skip flight segment, the performance index optimization is under the premise of satisfying the terminal constraints, i.e. the importance of satisfying the terminal constraints is greater than that of the performance index optimization.

[0111] The optimization process of the AGA algorithm is as follows: firstly, a control quantity dataset is constructed within the control constraint range, the control quantity (angle of attack and velocity tilt angle) segmentation time points are determined, and the control quantity-time history is smoothed by using a cubic spline interpolation; secondly, a penalty function is constructed to process the process constraints (dynamic pressure, overload, heat flow rate, and flight prohibited area, etc.) and terminal constraints, and a target function is established in combination with the trajectory optimization target (the farthest range / the farthest lateral range); finally, the control quantity is optimized by means of the adaptive genetic algorithm to obtain a trajectory (the farthest range / the farthest lateral range) that satisfies the constraint conditions and has the farthest range. It should be noted that in the process of unpowered skip flight, the maximum longitudinal range of skip flight is closely related to the maximum lift-drag ratio. Under the premise of not considering various constraints and other factors, generally speaking, the greater the lift-drag ratio, the greater the farthest range. Therefore, when determining the control quantity segmentation time points, the angle of attack corresponding to the maximum lift-drag ratio can be taken as a reference for segmentation. S3 is mainly to obtain a "quasi-optimal solution" close to the global optimal solution by running the AGA algorithm for a number of generations (which is much smaller than the number of generations when the adaptive genetic algorithm is used alone). The algorithm flow is shown in Figure 5 and specifically includes:

[0112] S31, population initialization coding: the feasible solution of the problem to be solved is expressed as a chromosome or individual in the genetic algorithm space by using real number coding;

[0113] S32, the fitness value of each individual in the population is calculated by the fitness function, the fitness function is an index for evaluating the chromosome, and is also the key to the "survival of the fittest" of the AGA algorithm. In the present patent, the Joines method is used to construct a comprehensive fitness function: the performance index at the terminal time and the terminal constraints are weighted to obtain the fitness function. The Joines method improves the static penalty function by setting the penalty coefficient to be dynamically self-adaptive, i.e. the penalty coefficient will be self-adaptively adjusted as the number of iterations increases. The specific form is as follows:

[0114]

[0115] In the formula: R1=(cd) η , f i (x) is the penalty function term for violating the i-th terminal constraint, and the specific form is as follows:

[0116]

[0117] wherein p i is the weight coefficient of the i-th terminal constraint penalty term.

[0118] Remaining constraint handling:

[0119] 1) Control constraints: Control constraints are automatically satisfied by constructing the control data set within the constraint limits.

[0120] 2) Process constraints (dynamic pressure, overload, heat flux rate and forbidden flight zone): The fitness function value of not satisfying the maximum dynamic pressure q dmax , maximum overload n Lmax , maximum heat flux rate and crossing the forbidden flight zone is set to infinity, i.e. F = inf at this time.

[0121] wherein c, η, are constants, and c = 0.5, is a feasible scheme; d is the number of iterations, J is the trajectory optimization objective function term, q is the corresponding weight coefficient, and K is a constant (to ensure that the fitness is positive).

[0122] S33, determine whether the maximum number of iterations is reached, if yes, obtain a quasi-optimal solution, and enter S4; otherwise, enter S34;

[0123] S34, perform genetic operations of selection, crossover and mutation on individuals according to the fitness value to generate a new generation of population, and enter S32;

[0124] Selection genetic operation: select excellent individuals from the population to form a new population to breed the next generation of individuals, and the higher the individual fitness value, the greater the probability of being selected. In this patent, the roulette method is selected, and the probability of individual i being selected is:

[0125]

[0126] Crossover genetic operation: Crossover operation is the most important operation of AGA algorithm, which is directly related to the global search ability of the algorithm. The core idea is to randomly select two individuals from the population, and combine the two individuals (chromosomes) to make the excellent characteristics of the parent excellent individuals inherit to the offspring, and produce new excellent individuals. The real number crossover method is selected in this patent, that is, the crossover operation of the jth gene of the kth chromosome a k and the jth gene of the lth chromosome a l is:

[0127]

[0128] Mutation genetic operation: In order to maintain the diversity of the population, by selecting a suitable mutation rate, to prevent premature convergence and avoid falling into local optimum. Mutation operation is the change of a gene value of an individual string in the population, the mutation operation method of the jth gene of the ith individual a ij is:

[0129]

[0130] Where, F i is the fitness value of individual i, Q is the number of population individuals, b is a random number in [0, 1], a max , a min are the upper and lower bounds of gene a ij , respectively; f(g) = r2(1-g0 / G max ) 2 , r2 is a random number, g0 is the current iteration number, G max is the maximum evolution number, and r is a random number in [0, 1].

[0131] S4, the quasi-optimal solution is used as the initial value of the GPM-SQP algorithm (including the terminal time and the initial and final time values of the control variable), and then the GPM-SQP algorithm is used to optimize the control variable in the guessed time period to obtain the global optimal solution or the final solution close to the global optimal solution.

[0132] The optimization process of GPM-SQP algorithm is as follows:

[0133] Firstly, the original optimal control problem is converted into NLP by using GPM algorithm, and the differential equation is converted into algebraic equation; then the SQP algorithm is used to divide the nonlinear constraint optimization problem into a series of quadratic programming sub-problems for solving. S4 mainly uses the terminal time and the initial and final time values of the control variable obtained by the AGA algorithm as the initial value of the GPM-SQP algorithm, so as to find the global optimal solution or the final solution close to the global optimal solution. The algorithm flow is shown in Figure 6 . Gauss-SQP algorithm includes two steps of variable discretization and parameterization and numerical optimization method, which are as follows:

[0134] S41, variable discretization and parameterization - Gauss pseudospectral method (GPM)

[0135] The core idea of the Gauss pseudospectral method is to discretize the state variable and the control variable on a series of Legendre-Gauss (LG) nodes at the same time, and to approximate the state variable and the control variable by using global polynomial interpolation, so that the system differential constraints are converted into algebraic constraints, and then the state variable and the control variable at the nodes are taken as the optimization design variables, so that the trajectory optimization problem (dynamic optimal control problem) is converted into NLP for solving, the specific steps are as follows:

[0136] ①, time domain conversion: let the starting time and the terminal time of the original dynamic optimal control problem be t0 and t f , and the GPM discrete points are selected in [-1, 1], so the following variable substitution is needed: x(τ) = (x(τ) - x0) / (x1 - x0) = x(τ) - x0

[0137]

[0138] ②, discrete state variable (state variable approximation): the collocation points τ k (k = 1, 2,...,N) are N LG points, and the collocation points are distributed in the interval (-1, 1), the roots (τ1, τ2,..., τ k ) of the N-order Legendre polynomial and τ0 = -1 are taken as the collocation nodes of the system, and N+1 Lagrange polynomials are constructed to approximate the state variable, which is expressed as:

[0139]

[0140] Among them, the discrete state value at the collocation node is exactly equal to the actual state value of the original system, that is, x(τ i ) = X(τ i ).

[0141] ③, discrete control variable (control variable approximation): the roots (τ1, τ2,..., τ k ) of the N-order Legendre polynomial are taken as the collocation nodes of the system, and N Lagrange polynomials are constructed to approximate the control variable, which is expressed as:

[0142]

[0143] ④, convert the differential equation at the collocation point into an algebraic equation, which is specifically:

[0144] For the approximate state variable, take the derivative with respect to τ

[0145]

[0146] The derivative of each Lagrange polynomial at the LG points can be represented by a differential approximation matrix:

[0147]

[0148] where:

[0149]

[0150] Thus, the differential equation constraints corresponding to the aircraft motion equations at the collocation points can be converted into algebraic constraints. That is:

[0151]

[0152] 5. Discretize the constraints;

[0153] Control constraints discretization: In GPM, the control variables only need to satisfy the constraints at the LG points, that is:

[0154]

[0155] Process constraints discretization: For process constraints, it can be considered that as long as the constraints are satisfied at the discrete nodes, the entire process satisfies the process constraints, that is:

[0156] C(X(τ k ),U(τ k ),τ k ; t0,t f )≤0,k=1,2,...,N (28)

[0157] Terminal state constraints discretization: The boundary conditions (initial state and terminal state) constraints can be described as follows:

[0158] φ(X(t0),t0,X(t f ),t f )=0 (29)

[0159] The initial state is fixed, and x(τ0)≡X(τ0)≡X(-1) is known. Since the terminal time τ f is not considered when approximating the state variable in equation (21), the terminal state needs to be supplemented. According to equation (26), we have:

[0160]

[0161] In addition, the state variable X(τ f ) at the terminal time of the system can be obtained by Gaussian integration through X k (k=0,1,2,...,N) and U k (k=1,2,...,N): ​

[0162]

[0163] where is a Gaussian weight, where then: τ k is the LG point.

[0164] 6. Discretize the farthest range performance index function.

[0165] The Bolza-type performance index function is:

[0166]

[0167] The integral term in the above equation is approximated by a Gaussian integral, and the performance index function in GPM is:

[0168]

[0169] Therefore, by:

[0170]

[0171] The original dynamic optimal control problem can be converted into a nonlinear programming problem (NLP), and the converted nonlinear programming problem has a consistent solution with the original control problem.

[0172] In summary, the trajectory optimization problem can be described as: determining the state variable X(τ k )(k=0, 1, 2,..., N), the control variable U(τ k )(k=1, 2,..., N), the initial time t0, the terminal time t f , within the time [t0, t f ] so that the performance index (33) is minimized under the above constraints (the converted aircraft motion algebraic equation constraint (26), the control constraint (27), the process constraint (28), and the terminal state constraint (31)).

[0173] S42, numerical optimization method—sequential quadratic programming (SQP)

[0174] For the discretized nonlinear programming problem (NLP), sequential quadratic programming is a method for directly processing constrained optimization problems. The basic idea is to determine a descending direction by solving a quadratic programming sub-problem in each iteration to reduce the value function to obtain a step size. Repeat these steps until the solution to the original problem is obtained.

[0175] Conclusion

[0176] Aiming at the trajectory optimization problem of hypersonic vehicles in the skip flight phase under complex constraint conditions including multiple forbidden zones, a solution based on GPM-AGA-SQP hybrid algorithm is proposed.

[0177] 1) In order to further quickly and accurately approach the "quasi-optimal solution" of the global optimal solution, three improvements are made to the fitness function in the genetic algorithm. A comprehensive fitness function in logarithmic form is constructed, which can be self-adapted with the increase of the number of population iterations. A new adaptive genetic algorithm is established.

[0178] 2) In order to solve the problem that GPM-SQP algorithm relies too much on initial guess value, the adaptive genetic algorithm is used to optimize the control variables globally, and the approximate global optimal solution is used as the initial guess value of GPM-SQP algorithm, so as to reduce the sensitivity of GPM-SQP algorithm to initial guess value.

[0179] 3) The GPM-AGA-SQP hybrid algorithm solves the problem of slow convergence speed and weak local optimization ability of AGA algorithm when approaching the global optimal solution, and combines the advantage of strong robustness of AGA algorithm, so that the hybrid algorithm can still complete the trajectory optimization task under the condition of model uncertainty and disturbance.

[0180] Experiments show that the established hybrid optimization algorithm GPM-AGA-SQP algorithm is effective and feasible, and can avoid multiple forbidden zones. Compared with the improved direct shooting method-AGA algorithm, the performance index value is improved by more than 16%. Compared with GPM-AACO-SQP algorithm, the control variable curve is smoother, the chattering amplitude is smaller, and it is more in line with the actual situation of avoiding forbidden zones. It is closest to the global optimal solution, and the performance index value is improved by more than 2.05%. At the same time, the parameter adjustment amount and difficulty are reduced, and it is easier to get the ideal trajectory that meets the constraint conditions and has the longest range. In addition, considering the model uncertainty and disturbance factors in the actual flight process, the Monte Carlo simulation experiment shows that all kinds of constraint conditions in the skip flight process can still meet the requirements, which shows that the constructed GPM-AGA-SQP algorithm has certain robustness. In terms of application, GPM-AGA-SQP algorithm is essentially an optimization method, which is not only suitable for trajectory optimization of hypersonic vehicles in the skip flight phase, but also suitable for general aircraft trajectory optimization problems, providing certain reference value for the next step of studying robust trajectory optimization problems with uncertainty quantization.

[0181] Simulation

[0182] The application carries out a large number of simulation verification based on the GPM-AGA-SQP hybrid algorithm, and verifies the effectiveness and feasibility of the constructed hybrid algorithm, and compares the simulation results of the GPM-AACO-SQP algorithm and the improved direct shooting method-AGA algorithm, to verify the superiority of the GPM-AGA-SQP algorithm, and finally verifies the robustness of the GPM-AGA-SQP algorithm by using Monte Carlo simulation test.

[0183] (1) Simulation related parameters

[0184] Before performing the trajectory optimization simulation of the high Mach number aircraft jump flight section, the aircraft and earth basic parameters, control constraints, process constraints, boundary constraints, related parameters of the established new AGA algorithm and the terminal time and control quantity initial and final time guess value obtained after the AGA algorithm completes the control quantity optimization need to be given.

[0185] The aircraft and earth basic parameters are shown in Table 1, the control constraints are shown in Table 2, the process constraints are shown in Table 3, and the boundary constraints (initial state and terminal state) are shown in Table 4:

[0186] Table 1 Simulation basic parameters

[0187]

[0188] Table 2 Control constraints

[0189]

[0190] Table 3 Process constraints

[0191]

[0192] Table 4 Boundary constraints (initial state and terminal state)

[0193]

[0194] Note: The initial state height of 50km corresponds to the speed of 10Ma≈3297.99m·s 1 , and the terminal state height of 27km corresponds to the speed of 3Ma≈899.40m·s -1 .

[0195] The main parameters of the adaptive genetic algorithm are shown in Table 5:

[0196] Table 5 Main parameters of adaptive genetic algorithm

[0197]

[0198] The AGA algorithm designed in this patent mainly improves the fitness function of the GA algorithm. Therefore, the other four parameters need to be set in advance. These four parameters are generally set within the following ranges: population size: 20-100; number of iterations: 100-500; crossover probability: 0.4-0.99; mutation probability: 0.0001-0.1. Considering that the AGA algorithm may exhibit premature convergence in actual simulations, missing some better solutions, this can be addressed by increasing the population size and mutation rate. To avoid excessive computation, the population size can be chosen as 500 (a 5-fold increase from the original, requiring trial and error), the number of iterations as 100 (selecting the minimum number of iterations to obtain a near-optimal solution), and the mutation probability as 0.1 (maximum mutation probability). Furthermore, choosing a crossover probability of 0.8±0.5 not only improves the accuracy and speed of the solution but also avoids failing to obtain a good solution due to an improper crossover probability selection. Therefore, the crossover probability can be chosen as 0.8 (median).

[0199] The optimal fitness of the AGA algorithm is as follows: Figure 7 As shown, Figure 7 The graph shows the variation of the optimal fitness value in the population with the number of iterations. As can be seen, the optimal fitness value fluctuates significantly in the early stages, gradually converging with increasing iterations, and stabilizing around generation 10. Slight variations occur between generations 60 and 90, but these are minor. This indicates that the AGA algorithm has good convergence performance for optimizing control variables.

[0200] After the adaptive genetic algorithm (AGA) optimizes the control input, the curves showing the changes in angle of attack and velocity tilt angle over time are shown below. Figures 8-9 As shown, Figures 8-9 It is based on a control quantity dataset constructed within the control constraints, with terminal constraints as the judgment condition, the farthest longitudinal or farthest transverse as the performance index, and considering process constraints (dynamic pressure, overload, heat flow rate, and no-fly zone, etc.). It uses an adaptive genetic algorithm to optimize the function values ​​of the control quantity at piecewise time points, and smooths the curves obtained by cubic spline interpolation of the control quantity-time history.

[0201] from Figures 8-9 As can be seen, using the longest longitudinal range as the performance index, the adaptive genetic algorithm optimizes the control input time history to 495s. The initial and final values ​​of the control input angle of attack are 0.175rad and 0.175rad, respectively, and the initial and final values ​​of the control input velocity tilt angle are -0.2975rad and -0.1575rad, respectively. These are used as the initial values ​​for solving the problem using the GPM-SQP algorithm with the longest longitudinal range as the performance index.

[0202] from Figures 8-9It can be seen that, taking the farthest lateral range as the performance index, the adaptive genetic algorithm optimization control quantity time history is 475s, the control quantity attack angle initial and final time values are 0.14rad and 0.1925rad respectively, and the control quantity speed inclination angle initial and final time values are 0.4025rad and-0.1575rad respectively, which are used as the initial value solved by the GPM-SQP algorithm.(2) GPM-AGA-SQP algorithm effectiveness and feasibility verification

[0203] The present application takes the farthest longitudinal range / farthest lateral range as the performance index, and a large number of simulations are carried out based on the GPM-AGA-SQP hybrid algorithm to verify the effectiveness and feasibility of the hybrid algorithm.

[0204] 1) Farthest longitudinal range trajectory optimization

[0205] In the trajectory optimization problem with the farthest longitudinal range as the target, the longitudinal range corresponds to the change of x, and the performance index is:

[0206] J0=min(-L0)=min(-x(t f )) (35)

[0207] When the GPM-AGA-SQP algorithm is used for farthest longitudinal range trajectory optimization, the GPOPS toolbox is used in combination with SNOPT software for solving, and the simulation result is as shown in Figures 10-24 .

[0208] 2) Farthest lateral range trajectory optimization

[0209] In the trajectory optimization problem with the farthest lateral range as the target, the lateral range corresponds to the change of z, and the performance index is:

[0210] J0=min(-L0)=min(-|z(t f )|) (36)

[0211] When the GPM-AGA-SQP algorithm is used for farthest longitudinal range trajectory optimization, the GPOPS toolbox is used in combination with SNOPT software for solving, and the simulation result is as shown in Figures 25-34 .

[0212] From Figures 16-18 , 31-33, it can be seen that the hybrid algorithm: GPM-AGA-SQP algorithm in the simulation with the farthest longitudinal range / farthest lateral range as the performance index, the control constraint is within the specified range; Figures 20-22 , 35-37, the process constraint (dynamic pressure, overload, heat flow rate) is also within the constraint condition range; Figures 23-24 , 38-39, it can be seen directly that the hybrid algorithm can realize the avoidance of multiple no-fly zones, which meets the actual demand.

[0213] To more intuitively reflect the constructed GPM-AGA-SQP hybrid algorithm meeting the terminal constraints, Table 6 is established for verification.

[0214] Table 6 Terminal state corresponding to terminal constraint condition of GPM-AGA-SQP algorithm

[0215]

[0216] From Figures 10-11 , 14 and 25-26, 29 and Table 6, it can be seen that in the simulation taking the maximum range / maximum range as the performance index, the hybrid algorithm can meet the terminal constraint condition. In summary, the effectiveness and feasibility of the constructed GPM-AGA-SQP hybrid algorithm are illustrated.

[0217] (3) Superiority verification of GPM-AGA-SQP algorithm

[0218] The present application will take the maximum range / maximum range as the performance index, and based on the GPM-AGA-SQP algorithm, the GPM-AACO-SQP algorithm and the improved direct shooting method-AGA algorithm for comparative simulation, in order to verify the superiority of the GPM-AGA-SQP algorithm.

[0219] It should be noted that the adaptive ant colony optimization (AACO) is improved based on the ACO algorithm, and the global search / local search step and the penalty term coefficient in the fitness function are set to be self-adaptive adjustment form with the increase of iteration number. The first stage of the GPM-AACO-SQP algorithm is to obtain a "quasi-optimal solution" close to the global optimal solution by using the AACO algorithm, and the core idea of the second stage is similar to that of the GPM-AGA-SQP algorithm. The improved direct shooting method-AGA algorithm is to optimize the control amount at the discrete point in the improved direct shooting method by AGA algorithm, and to obtain the state variable describing the motion trajectory based on the parameterized control amount. Among them, the improved direct shooting method is proposed for the problem that the direct shooting method has high requirement for initial value and may produce local minimum, which takes the control amount initial and final time value into the optimization design variable, takes the composite constraint term (such as speed and height) in the terminal constraint as the judgment condition, and takes the end time as the terminal time, so as to achieve the purpose of reducing the sensitivity to initial value.

[0220] 1) Performance index comparison

[0221] The simulation comparison results of the performance indexes of the three algorithms are shown in Figures 40-41 , Table 7.

[0222] Table 7 Comparison of optimal performance indexes of three algorithms

[0223]

[0224] From Figures 40-41 , Table 7, it can be seen that in the performance index change trend, GPM-AGA-SQP algorithm and GPM-AACO-SQP algorithm are relatively close in change trend, compared with the improved direct targeting method-AGA algorithm, the change trend in the early stage is basically the same, but the trend in the later stage is quite different, which to some extent shows that the improved direct targeting method-AGA algorithm has the defect of poor local search ability when approaching the global optimal solution; in the performance index value, the GPM-AGA-SQP algorithm improves the farthest longitudinal value by 2.05% and the farthest lateral value by 5.64% compared with the GPM-AACO-SQP algorithm, and the farthest longitudinal value by 16.82% and the farthest lateral value by 18.63% compared with the improved direct targeting method-AGA algorithm.

[0225] 2) Comparison of control variable optimization ability

[0226] The control variable change curves of the three algorithms are shown in Figures 42-45 From Figures 42-45 , it can be seen that the GPM-AGA-SQP algorithm and the GPM-AACO-SQP algorithm are basically consistent in the change trend of the control variable, but compared with the GPM-AACO-SQP algorithm, the control variable change curve of the GPM-AGA-SQP algorithm is smoother and the chattering amplitude is smaller, which is more in line with the actual situation of avoiding the no-fly zone and closest to the global optimal solution; compared with the GPM-AGA-SQP algorithm, the improved direct targeting method-AGA algorithm shows a certain "inertia" in the change trend of the control variable, that is, the small fluctuation of the control variable in the GPM-AGA-SQP algorithm corresponds to the delay fluctuation or large fluctuation of the control variable in the improved direct targeting method-AGA algorithm, which to some extent also reflects the optimization search ability of the control variable of the improved direct targeting method-AGA algorithm, that is, it is strong in large-scale global search ability, but weak in local search ability. It should be noted that the control variable fluctuation is near the "quasi-optimal solution" of the global optimal solution, not random fluctuation, which is also the reason why the terminal time value and the initial and final time value of the control variable obtained by the AGA algorithm can be used as the initial guess value of the GPM-SQP algorithm.

[0227] 3) No-fly zone avoidance ability and test parameters

[0228] The three algorithms avoid the no-fly zone as shown in Figures 46-47 From Figures 46-47From Table 8, it can be seen that the GPM-AGA-SQP algorithm, the GPM-AACO-SQP algorithm, and the improved direct shooting method-AGA algorithm can all achieve evasion of multiple no-fly zones. It should be noted that the GPM-AGA-SQP algorithm and the GPM-AACO-SQP algorithm provide initial and final time control quantities and terminal time guess values through the AGA algorithm and the AACO algorithm, respectively, and the test parameters are mainly concentrated in the AGA and AACO algorithm parts. Among them, the AGA algorithm actually needs 2 test parameters, i.e., the performance index weight coefficient q and the population size Q; the AACO algorithm actually tests 6 parameters, which increases the pheromone evaporation coefficient Rho, the global transfer selection factor P0, and the search step length adjustment factor step (to ensure that the search step length does not exceed the boundary and to make the step length adjustment within a certain accuracy range). The following compares the test parameters of the three algorithms, as shown in Table 8:

[0229] Table 8 Test parameters of three algorithms

[0230]

[0231]

[0232] From Table 8, it can be seen that the constructed hybrid algorithm: GPM-AGA-SQP algorithm has 3 parameters to adjust, among which the terminal time value obtained based on the AGA algorithm can be fine-tuned (simulation verification shows that the fine-tuning range is about ±3s). Compared with the GPM-AACO-SQP algorithm, the GPM-AGA-SQP algorithm has a lower parameter adjustment amount (6 parameters for the GPM-AACO-SQP algorithm) and difficulty (simulation verification shows that the fine-tuning range of the terminal time value based on the AACO algorithm is about ±5s), and it is easier to obtain a trajectory that meets the constraint conditions and has optimal performance (maximum range / maximum lateral range).

[0233] (4) Robustness verification of GPM-AGA-SQP algorithm

[0234] During the jump flight of a high-Mach number aircraft, the influence of model uncertainty and disturbance factors on trajectory optimization performance needs to be considered. A 50-time Monte Carlo simulation test is conducted through the GPM-AGA-SQP hybrid algorithm to verify its robustness. Among them, in terms of model uncertainty, due to factors such as manufacturing shape, measurement error, etc., it is more reasonable to consider the aerodynamic parameters C x and C y as uncertainties. The uncertainty and disturbance are shown in Table 9.

[0235] Table 9 Uncertainty factors and disturbance parameters for Monte Carlo simulation

[0236]

[0237] Under the influence of the above uncertainty parameters and disturbances, the mean value, maximum value and minimum value of the performance index, process constraint term (including heat flow rate, dynamic pressure and overload), terminal constraint term (including terminal velocity, terminal height and terminal ballistic inclination angle) parameters obtained by the GPM-AGA-SQP hybrid algorithm are shown in Tables 10-11, and the Monte Carlo simulation trajectory is shown in Figures 48-49 .

[0238] Table 10 Monte Carlo simulation results (performance index: maximum range)

[0239]

[0240] Table 11 Monte Carlo simulation results (performance index: maximum cross-range)

[0241]

[0242] From Tables 10-11, Figures 48-49 it can be concluded that the hybrid algorithm GPM-AGA-SQP algorithm constructed has certain robustness under the influence of model uncertainty and disturbance factors, and in the simulation taking maximum range / maximum cross-range as the performance index, the control constraints, terminal constraints and process constraints (dynamic pressure, overload, heat flow rate and flight-prohibited area) basically meet the conditions.

[0243] The technical solutions of the present application are not limited to the restrictions of the above specific embodiments, and any technical variations made according to the technical solutions of the present application fall within the protection scope of the present application.

Claims

1. A method for optimizing the trajectory of a high Mach number aircraft during a skip flight segment, characterized in that, include: S1. Ignoring Earth's rotation and curvature, and assuming Earth to be a plane, dynamic and kinematic equations are established in both the ballistic and ground coordinate systems, resulting in a simplified model of the center-of-mass motion equations for a 3-DOF high Mach number aircraft. The simplified model of the center-of-mass motion equations is as follows: （1）， In the formula: V is the velocity of the aircraft, θ and α represents the trajectory inclination angle and trajectory deflection angle, respectively; y represents altitude; x and z represent the longitudinal and transverse strokes, respectively; m represents the aircraft mass; g represents gravitational acceleration; and α represents the angle of attack. α is the velocity tilt angle; B To balance the angle of attack, β B To balance the sideslip angle; Y B and Z B α B β B The corresponding equilibrium lift and equilibrium lateral force, where X is the drag, are calculated using the following expressions: （2）， In the formula: C x C y For drag coefficient and lift coefficient, , where ρ is the air density at the altitude of the high Mach number aircraft, and S is the reference area of ​​the aircraft; S2. Determine the constraints and objective function to be optimized during the hop flight segment, and describe the trajectory optimization problem as a standard dynamic optimal control problem. The constraints include self-control constraints, process constraints during the skip flight segment, and terminal constraints upon reaching the destination. The objective function is selected as the longest flight distance, as detailed below: Control constraints include angle of attack α and velocity tilt angle. The constraint is represented as: （3）， Considering practical engineering applications, the rate of change of angle of attack is typically... Certain restrictions are also necessary, namely: （4）， At this point, the angle of attack change rate can be... It is treated as a virtual control quantity, and the original control quantity angle of attack α is treated as a state variable. Process constraints include heat flow rate constraints, dynamic pressure constraints, overload constraints, and no-fly zone constraints, as detailed below: The heat flow rate constraint is expressed as: （5）， Dynamic pressure constraint is expressed as: （6）， Overload constraints are expressed as: （7）， Assume the no-fly zone that the high Mach number aircraft avoids during its hop flight is cylindrical, with its center point being... The radius is r i The specific expression is as follows: （8）， set up The shortest distance from a high Mach number aircraft to a no-fly zone. If the variable is an infinitesimal normal, then the no-fly zone constraint can be expressed as: In the formula These are constant parameters related to high Mach number aircraft. For the maximum permissible heat flux, For the maximum permissible dynamic pressure, The maximum allowable overload; Terminal constraints include altitude, velocity, and trajectory inclination constraints at the terminal moment, expressed as: （9）， In the formula, , , These are the altitude, velocity, and trajectory inclination at the terminal moment of the jump flight; , , These are the altitude, speed, and trajectory inclination at the predetermined jump-flight terminal moment; , , These are the upper bounds of the errors for the three constraint quantities; The objective function is expressed as: （10）， Where L df Refers to the range of an aircraft; The trajectory optimization problem for the skip flight segment of a high Mach number aircraft can be described as follows: under the conditions of satisfying the equations of motion of the center of mass (1), control constraints (3)~(4), process constraints (5)~(8), and terminal constraints (9), find the control variables angle of attack α and velocity tilt angle. The optimal time series minimizes the objective function (10); a more general optimal control problem is described as follows: finding the control variables This minimizes the Bolza-type performance index, which is a combination of the Mayer and Lagrange types, and is expressed as: （11）， And it satisfies the constraints of the equations of motion of the center of mass. （12）， Boundary condition constraints （13）， and inequality constraints （14）， Where t0 is the starting time, t f For the terminal time; S3. Based on the transformed standard dynamic optimal control problem, we first utilize the strong global search capability of the AGA algorithm to perform global optimization on the control quantity, and obtain a near-optimal solution that is close to the global optimal solution. S4. Use the near-optimal solution as the initial value of the GPM-SQP algorithm, and then use the GPM-SQP algorithm to optimize the control quantity within the conjecture time period to obtain the global optimal solution or the final solution close to the global optimal solution.

2. The trajectory optimization method for the skip flight segment of a high Mach number aircraft according to claim 1, characterized in that, In S3, the Adaptive Genetic Algorithm (AGA) is used to find a near-optimal solution to the trajectory optimization problem, thus addressing the sensitivity of the GPM-SQP algorithm to initial conditions. Specifically, this includes: S31. Population initialization encoding: The feasible solution to the problem to be solved is represented as a chromosome or individual in the genetic algorithm space using real number encoding. S32. Constructing a comprehensive fitness function using the Joints method: This involves weighting the performance metrics at the terminal time step with the terminal constraints to obtain the fitness function, and then calculating the fitness value of each individual in the population using this function. The fitness function is an indicator for evaluating chromosomes and is a key aspect of the "survival of the fittest" principle in the AGA algorithm. The Joints method improves upon the static penalty function by setting the penalty coefficient to be dynamically adaptive, meaning that the penalty coefficient adjusts adaptively as the number of iterations increases. Specifically: （15）， In the formula: f i (x) is the penalty function term for violating the i-th terminal constraint, with the following specific form: （16）， In the formula: p i Let be the weight coefficient of the penalty term for the i-th terminal constraint, and handle the remaining constraints as follows: 1) Control constraints: By constructing a control variable dataset within the constraint conditions, the control constraints can be automatically satisfied; 2) Process constraints: those that do not meet the maximum dynamic pressure requirement will be subject to these constraints. Maximum overload Maximum heat flux And the fitness function value for crossing the no-fly zone is set to infinity, that is, at this time ; In the formula: c, , All are constants. Numerical simulations have verified that c = 0.5 is appropriate. d is the number of iterations, J is the objective function term for trajectory optimization, q is the corresponding weight coefficient, and K is a constant that ensures the fitness is positive. The terminal state of the jump flight segment is strictly constrained and the performance index at the terminal moment is required to meet the optimal requirements. The optimal performance index is under the premise of meeting the terminal constraints. Based on LSP, after logarithmic processing, the differences between terminal constraint terms at the same level are small. The weight coefficients of the penalty terms of the terminal constraints are considered to have the same weight and are the same constant value. On this basis, the weight coefficient q of the next level performance index can be adjusted proportionally. That is, it is necessary to test the weight coefficient q to ensure that the terminal constraint terms are within the trajectory range. S33. Determine if the maximum number of iterations has been reached. If yes, obtain the near-optimal solution and proceed to S4; otherwise, proceed to S34. S34. Genetic operations such as selection, crossover, and mutation are performed on individuals based on fitness values ​​to generate a new generation of population, which then enters S32. Selection genetic operation: Selecting superior individuals from the population with a certain probability to form a new population for reproduction, resulting in the next generation of individuals. The higher the fitness value of an individual, the greater the probability of being selected; selection roulette method, individuals... The probability of being selected is: （17）， Crossover genetic operation: The crossover operation is the most important operation in the AGA algorithm, directly affecting its global search capability. The core idea is to randomly select two individuals from the population and crossover them (or chromosomes), allowing the superior traits of the parent individuals to be inherited by the offspring, producing new superior individuals. Real number crossover is used, i.e., the k-th chromosome a... k and the Chromosomes The crossover operation at position j is as follows: （18）， Mutation genetic operations: To maintain population diversity, an appropriate mutation rate is selected to prevent premature convergence and avoid getting trapped in local optima. A mutation operation involves changing the value of a specific gene in an individual's string within the population. For example, the mutation operation modifies the value of the j-th gene 'a' in the i-th individual. ij The method for performing mutation operations is as follows: （19）， Among them, F i Let be the fitness value of individual i, Q be the number of individuals in the population, b be a random number in the range [0,1], and a be the fitness value of individual i. max a min Gene a ij The upper and lower bounds; r2 is a random number, g0 is the current iteration number, and G max The maximum number of evolutions is r, which is a random number in the range [0,1].

3. The trajectory optimization method for the skip flight segment of a high Mach number aircraft according to claim 2, characterized in that, In S4, the near-optimal solution is used as the initial value for the GPM-SQP algorithm. Then, the GPM-SQP algorithm is used to optimize the control quantity within the conjectured time period to obtain the global optimum or a final solution close to the global optimum. The Gauss-SQP algorithm includes two steps: variable discretization and parameterization, and numerical optimization methods, as detailed below: S41. Variable Discretization and Parameterization – Gaussian Pseudospectral Method, i.e., GPM ① Time Domain Transformation: Let the start time and terminal time of the original dynamic optimal control problem be t0 and t1, respectively. f ,and The discrete points of GPM are chosen in the range [-1, 1]. Therefore, the following variable substitution is required when solving the problem: (expressed as:) （20）， ② Discrete state variables, i.e., state variable approximation: collocation of GPM Given N LG points, with the point distribution interval (-1, 1), find the roots of the Nth-order Legendre polynomial. and As the system's configuration node, N+1 Lagrange polynomial approximation state variables are constructed, represented as follows: （21）， in, The discretized state values ​​at the configuration nodes are exactly equal to the actual state values ​​of the original system, that is... ; ③ Discrete control variables, i.e., control variable approximation: the roots of the Nth-order Legendre polynomial As the system's configuration node, N Lagrange polynomial approximation control variables are constructed, represented as follows: （22）， ④ Convert the differential equation at the collocation point into an algebraic equation, specifically: For approximate state variables, perform an assessment of... Taking the derivative, we get: （23）， The derivative of each Lagrange polynomial at the LG point can be represented by a differential approximation matrix: （24）， In the formula: ,but: （25）， Therefore, the differential equation constraints corresponding to the equations of motion of the aircraft at the collocation point can be converted into algebraic constraints, that is: （26）， ⑤ Discretize the constraints, specifically as follows: Discretization of control constraints: In GPM, the control quantity only needs to satisfy the constraint conditions at the LG point, that is: （27）， Process constraint discretization: For process constraints, it is assumed that as long as the constraints are satisfied at discrete nodes, the process constraints are satisfied throughout the entire process, that is: （28）， Discretization of terminal state constraints: Boundary conditions, i.e., initial state and terminal transition constraints, can be described as follows: （29）， The initial state is fixed and known. Since equation (21) does not consider the terminal time when approximating the state variables, Therefore, the terminal state needs to be supplemented. According to equation (26), we can obtain: （30）， In addition, the system's terminal state variables can be and By Gaussian integration, we obtain: （31）， In the formula: The weights are Gaussian, where ,but: , Point LG; ⑥ Discretize the maximum range performance index function, specifically as follows: Bolza-type performance index function: （32）， The integral term in the formula is approximated by Gaussian integration, yielding the performance index function in GPM: （33）， Therefore, through: （34）， The original dynamic optimal control problem can be transformed into a nonlinear programming problem, and the transformed nonlinear programming problem has a consistent solution with the original control problem. In summary, the trajectory optimization problem can be described as: in [t0,t... f Within a given time period, determine the state variables at discrete points. Control variables Initial time t0, terminal time t f This ensures that, under the above constraints, namely the transformed algebraic equations of motion of the aircraft (26), control constraints (27), process constraints (28), and terminal state constraints (31), the performance index (33) is minimized. S42. Numerical Optimization Method – Sequential Quadratic Programming, i.e., SQP Algorithm For discretized nonlinear programming problems, sequential quadratic programming is a method that directly handles constrained optimization problems. The idea is to determine a descent direction by solving a quadratic programming subproblem in each iteration, thereby reducing the value function to obtain the step size. These steps are repeated until the solution to the original problem is obtained.

Images