CN115981372A

CN115981372A - high-Mach-number aircraft jumping flight segment trajectory optimization method

Info

Publication number: CN115981372A
Application number: CN202310127109.0A
Authority: CN
Inventors: 温杰; 周欢; 丛戎飞; 丁智坚; 张旭
Original assignee: General Engineering Research Institute China Academy of Engineering Physics; Institute of Aerospace Technology of China Aerodynamics Research and Development Center
Current assignee: General Engineering Research Institute China Academy of Engineering Physics; Institute of Aerospace Technology of China Aerodynamics Research and Development Center
Priority date: 2023-02-16
Filing date: 2023-02-16
Publication date: 2023-04-18

Abstract

The invention discloses a high-Mach-number aircraft jumping flight segment trajectory optimization method, which relates to the field of aircraft trajectory optimization and comprises the steps of S1, establishing a kinetic equation and a kinematic equation to obtain a centroid motion equation set simplified model; s2, determining constraint conditions and a target function to be optimized, and describing a track optimization problem as a standard dynamic optimal control problem; s3, carrying out global optimization on the control quantity by using an AGA algorithm to obtain a quasi-optimal solution of a global optimal solution, so as to solve the problem that the GPM-SQP algorithm is sensitive to initial value abnormity when solving the nonlinear programming problem; and S4, the quasi-optimal solution is used as an initial value of the GPM-SQP algorithm, then the GPM-SQP algorithm is used for optimizing the control quantity in a guessing time period, and the global optimal solution or the final solution close to the global optimal solution is obtained based on the advantages of high calculation precision and high convergence speed of the GPM-SQP algorithm. The method comprehensively utilizes the advantages of strong global search capability of the AGA algorithm and high convergence speed of the GPM-SQP algorithm, and easily obtains an ideal track which meets constraint conditions and has the farthest flight distance.

Description

high-Mach-number aircraft jumping flight segment trajectory optimization method

Technical Field

The invention relates to the field of optimization of flight vehicle tracks, in particular to a track optimization method for a jumping flight segment of a high-mach-number flight vehicle.

Background

The high-Mach number aircraft has important and gross war values such as high-altitude reconnaissance, high-speed penetration, remote accurate attack and the like, and becomes a hotspot of research in the aerospace field. The 'jumping flight segment' is a section of trajectory of the high Mach number aircraft which decelerates and descends after returning to the atmosphere until reaching a preset target. In order to realize that the high-Mach-number aircraft accurately hits a target to complete a flight task, the optimization of the track of the jumping flight section is the most important one, and the track which meets various constraint conditions and has the optimal performance index is planned by solving the optimal time sequence of control variables and motion states on the basis of an aircraft dynamic model.

The following difficulties exist in solving the optimization problem of the jumping flight segment track: firstly, the motion model of the aircraft is nonlinear and strongly coupled, and certain difficulty exists in model establishment and track optimization; secondly, the flight environment is complex and changeable, the flight track is interfered by model uncertainty, disturbance factors and the like, a certain requirement is put forward on the robustness of a track optimization algorithm; thirdly, not only the constraint limits such as dynamic pressure, overload and heat flow rate need to be considered in the jumping flight process, but also multiple no-fly zones formed by factors such as natural environment, politics and military need to be avoided, and in addition, the strict limit of the terminal state provides certain challenges for the research of the track optimization problem.

The trajectory optimization method can be divided into an indirect method and a direct method according to whether the performance index is directly optimized. The essence of the indirect method is that based on a classical variational method or Pontryagin minimum Value principle, a prime dynamic optimal control Problem is converted into a Hamilton two-point Boundary Value Problem (HBVP) to solve, but for a high-Mach aircraft with complex dynamic characteristics, the Boundary Value Problem is small in solution convergence domain, and initial Value estimation of a covariance variable is difficult. Unlike indirect methods, the direct method is to disperse and parameterize variables in the dynamic optimal control problem, convert the dynamic optimal control problem into a Nonlinear Programming problem (NLP), and further combine a numerical optimization method to solve, without deriving a first-order necessary condition, and the application is relatively wide.

Disclosure of Invention

The invention aims to solve the problems and designs a method for optimizing the track of the jumping flight section of the high-Mach aircraft.

The invention achieves the above purpose through the following technical scheme:

a high-Mach-number aircraft jumping flight segment trajectory optimization method comprises the following steps:

s1, neglecting the earth rotation and curvature influence, assuming the earth as a plane, and establishing a kinetic equation and a kinematic equation in a trajectory coordinate system and a ground coordinate system respectively to obtain a simplified model of a mass center motion equation set of the 3-degree-of-freedom high-Mach aircraft;

s2, determining constraint conditions borne by the jumping flight segment in the flight process and a target function to be optimized, and describing a track optimization problem as a standard dynamic optimal control problem;

s3, based on the converted standard dynamic optimal control problem, firstly, utilizing the advantage of strong global search capability of the AGA algorithm to carry out global optimization on the control quantity to obtain a quasi-optimal solution close to the global optimal solution;

and S4, taking the quasi-optimal solution as an initial value of the GPM-SQP algorithm, and optimizing the control quantity in a guessing time period by using the GPM-SQP algorithm to obtain a global optimal solution or a final solution close to the global optimal solution.

The invention has the beneficial effects that: the method comprehensively utilizes the strong global search capability of an Adaptive Genetic Algorithm (AGA) and the advantage of high convergence speed of the GPM-SQP algorithm, and easily obtains an ideal track which meets constraint conditions and has the farthest voyage. The method is characterized in that a track optimization process is completed off line, although the GPM-AGA-SQP hybrid algorithm provided by the invention is provided by taking a high-Mach aircraft as a research background and aiming at the defect that the GPM-SQP algorithm is sensitive to an initial guess value when solving the problem of track optimization of a jumping flight segment of the high-Mach aircraft under the condition of complex multiple constraints, the algorithm also solves the problems of low convergence speed and weak local optimization capability when the AGA algorithm approaches to a global optimal solution, and simultaneously combines the advantage of strong robustness of the AGA algorithm, so that the constructed GPM-AGA-SQP hybrid algorithm can still complete a post-track optimization task under the conditions of model uncertainty and disturbance interference. Finally, the GPM-AGA-SQP hybrid algorithm is essentially an optimization method, and is not only suitable for optimizing the jumping flight section track of the high-Mach-number aircraft, but also suitable for optimizing the general flight track.

Drawings

FIG. 1 is a drag coefficient;

FIG. 2 is a lift coefficient;

FIG. 3 is a lift-to-drag ratio;

FIG. 4 is a schematic flow chart of the genetic algorithm of the present invention;

FIG. 5 is a flow chart of the "quasi-optimal solution" algorithm for obtaining a globally optimal solution by the AGA algorithm of the present invention;

FIG. 6 is a flow chart of the present invention for rapidly finding a global optimal solution by using the quasi-optimal solution as the initial solution of the GPM-SQP algorithm;

FIG. 7 shows the optimal fitness of the adaptive genetic algorithm when the optimization method for the jumping flight segment trajectory of the high-Mach number aircraft is used for simulation;

FIG. 8 is a time-varying curve of the angle of attack after the adaptive genetic algorithm completes optimizing the control variables during simulation;

FIG. 9 is a time-varying curve of the velocity tilt angle after the adaptive genetic algorithm has completed optimizing the control quantity during simulation;

FIG. 10 is a graph of simulated velocity variations in a trajectory optimization problem targeting the farthest longitudinal voyage;

FIG. 11 is a trajectory inclination change curve for simulation in a trajectory optimization problem targeting the farthest longitudinal voyage;

FIG. 12 is a trajectory deviation angle variation curve for simulation in a trajectory optimization problem targeting the farthest longitudinal voyage;

FIG. 13 is a plot of the course change for the simulation in the trajectory optimization problem targeting the farthest longitudinal course;

FIG. 14 is a graph of simulated height variations in a trajectory optimization problem targeting a furthest longitudinal voyage;

FIG. 15 is a plot of course change for a simulation in a trajectory optimization problem targeting a farthest longitudinal course;

FIG. 16 is a graph of simulated angle of attack variation in a trajectory optimization problem targeting the farthest longitudinal voyage;

FIG. 17 is a plot of simulated rate of change of angle of attack in a trajectory optimization problem targeting a farthest longitudinal voyage;

FIG. 18 is a graph of simulated velocity tilt angle changes in a trajectory optimization problem targeting the farthest longitudinal voyage;

FIG. 19 is a graph of simulated height versus longitudinal variation in a trajectory optimization problem targeting a furthest longitudinal voyage;

FIG. 20 is a graph of simulated heat flow rate variation in a trajectory optimization problem targeting the farthest longitudinal voyage;

FIG. 21 is a dynamic pressure variation curve for simulation in a trajectory optimization problem targeting the farthest longitudinal voyage;

FIG. 22 is a graph of simulated overload variation in a trajectory optimization problem targeting the furthest longitudinal voyage;

fig. 23 is a plan view of a flight avoidance area for simulation in a trajectory optimization problem targeting the farthest longitudinal course;

fig. 24 is a three-dimensional perspective view of a flight-forbidden region evaded by simulation in a trajectory optimization problem with the farthest longitudinal voyage as a target;

FIG. 25 is a graph of simulated velocity variations in a trajectory optimization problem targeting the farthest lateral leg;

FIG. 26 is a graph of ballistic dip variation simulated in a trajectory optimization problem targeting the farthest lateral voyage;

fig. 27 is a trajectory deviation angle change curve for simulation in a trajectory optimization problem targeting the farthest lateral voyage;

FIG. 28 is a graph of the longitudinal variation for simulation in a trajectory optimization problem targeting the farthest lateral voyage;

FIG. 29 is a graph of simulated altitude changes in a trajectory optimization problem targeting the farthest lateral leg;

FIG. 30 is a graph of course change for a simulation in a trajectory optimization problem targeting the farthest lateral course;

FIG. 31 is a graph of simulated angle of attack variation in a trajectory optimization problem targeting the farthest lateral range;

FIG. 32 is a graph of simulated rate of change of angle of attack in a trajectory optimization problem targeting the farthest lateral range;

FIG. 33 is a graph of simulated velocity tilt angle variation in a trajectory optimization problem targeting the farthest lateral voyage;

FIG. 34 is a graph of simulated height-course variation in a trajectory optimization problem targeting the farthest lateral course;

FIG. 35 is a graph of simulated heat flow rate variation in a trajectory optimization problem targeting the farthest lateral voyage;

FIG. 36 is a dynamic pressure variation curve for simulation in a trajectory optimization problem targeting the farthest lateral voyage;

FIG. 37 is a graph of simulated overload variation in a trajectory optimization problem targeting the farthest lateral leg;

fig. 38 is a plan view of an avoidance no-fly zone for simulation in a trajectory optimization problem targeting the farthest lateral voyage;

FIG. 39 is a three-dimensional view of a flight-forbidden region evaded by simulation in a trajectory optimization problem targeting the farthest lateral voyage;

FIG. 40 is a graph showing the comparison of performance indexes of three algorithms (performance index: farthest longitudinal range);

FIG. 41 shows the comparison of the performance indicators of the three algorithms (performance indicator: farthest horizontal run);

FIG. 42 is a graph showing the comparison of the control variables and the attack angles of three algorithms (performance index: farthest longitudinal distance);

FIG. 43 shows a control quantity speed-tilt angle comparison curve (performance index: farthest longitudinal) for three algorithms;

FIG. 44 shows the control quantity-attack angle contrast curves (performance index: farthest traverse);

FIG. 45 is a graph showing the relationship between the control quantity, the speed and the inclination angle (performance index: farthest traverse);

FIG. 46 shows three algorithms for avoiding the top view contrast curve (performance index: farthest longitudinal range) of the no-fly zone;

FIG. 47 shows a comparison curve of top views of a no-fly zone (performance index: farthest horizontal distance) avoided by three algorithms;

FIG. 48 is a graph showing simulated flight trajectories (performance index: farthest longitudinal distance);

FIG. 49 Monte Carlo simulation flight trajectory (performance index: farthest traverse).

Detailed Description

In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention. It is to be understood that the embodiments described are only a few embodiments of the present invention, and not all embodiments. The components of embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations.

Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

The following detailed description of embodiments of the invention refers to the accompanying drawings.

A high Mach number aircraft jumping flight segment trajectory optimization method comprises the following steps:

s1, neglecting the earth rotation and curvature influence, assuming the earth as a plane, and respectively establishing a kinetic equation and a kinematic equation in a trajectory coordinate system and a ground coordinate system to obtain a simplified model of a 3-degree-of-freedom high-Mach aircraft centroid motion equation set:

in the formula: v (t) is the aircraft speed, theta (t) and psi _v (t) ballistic dip and ballistic declination, respectively; y (t) is height, and x (t) and z (t) are respectively longitudinal course and transverse course; m (t) is the aircraft mass, g is the gravitational acceleration, α (t) is the angle of attack, γ _v (t) is the velocity ramp angle. Alpha (alpha) ("alpha") _B To balance angle of attack, beta _B To balance the sideslip angle; y is _B And Z _B Are each alpha _B 、β _B Corresponding balance lift force and balance lateral force, X is resistance, and the calculation expression is as follows:

in the formula: c _x 、C _y Q =0.5 ρ V for drag coefficient and lift coefficient ² Where ρ is the air density at the altitude where the high mach number aircraft is located, and S is the aircraft reference area.

As can be seen from the equation (2), the aerodynamic coefficient (including the coefficient of resistance C) _x And coefficient of lift C _y ) The accuracy degree of the method directly influences whether aerodynamic force (including drag force and lift force) calculation is accurate, the method takes an HTV-like aircraft as a research object, obtains aerodynamic coefficients under different flight states through numerical simulation calculation and evaluation, and is shown in figures 1-3, wherein reference lengths and reference areas in the calculation of the lift coefficient and the drag coefficient are respectively 1.0m and 1.0m ² 。

And S2, determining a constraint condition borne by the jumping flight segment in the flight process and a target function to be optimized, and describing the track optimization problem as a standard dynamic optimal control problem. The constraint conditions comprise self control quantity constraint, process constraint in the process of jumping flight segment flight and terminal constraint of reaching a destination, and the objective function selects the farthest flight path (the farthest longitudinal path/the farthest transverse path) as follows:

the control quantity constraints include an angle of attack alpha and a speed tilt angle gamma _v . To ensure the stable flight of the aircraft, the variable attack angle alpha and the speed inclination angle gamma are controlled _v It needs to be limited within a certain range:

considering the actual application of the project, the angle of attack change rate is usually

Certain limitations also need to be imposed, namely:

at this time, the attack angle change rate can be adjusted

The control quantity is regarded as a virtual control quantity, and the original control quantity attack angle alpha is regarded as a state variable to be processed.

Process constraints include heat flow rate constraints, dynamic pressure constraints, overload constraints, and no-fly zone constraints. The high mach number aircraft has a significant pneumatic heating effect in the process of jumping flight, and in order to ensure the structural safety of the aircraft, the heat flow rate needs to be constrained:

huge dynamic pressure can be produced to the jump flight in-process, for preventing that the aircraft structure from becoming invalid, need retrain the dynamic pressure:

q＝0.5*ρ*V ² ≤q _dmax (6)

in order to ensure the normal operation of relevant equipment of the aircraft, the overload of the aircraft needs to be limited within a certain range:

the no-fly zone is an area which is formed by natural environment, politics and military factors and is not allowed to fly into, the no-fly zone is a typical regular cylindrical shape, and the central point is (x) _i ,y _i ,z _i ) Radius r _i The specific expression is as follows:

when the high-Mach number aircraft jumps, the track of the high-Mach number aircraft cannot enter the no-fly zone, namely the shortest distance between the aircraft and the no-fly zone must be ensured to be a positive value, L (x (t), t) is the shortest distance between the aircraft and the no-fly zone, and epsilon _n For infinitesimal normal quantities, the no-fly zone constraint may be expressed as:

L(x(t),t)≥ε _n

in the formula K _Q Is a constant parameter associated with high mach number aircraft,

to the maximum allowable heat flow rate, q _dmax For maximum dynamic pressure allowed, n _Lmax Maximum overload allowed;

terminal constraints include altitude, velocity and ballistic dip constraints at the time of the terminal. In order to meet the requirement of middle and terminal guidance shift change, the height, the speed and the trajectory inclination angle of a terminal in a jumping flight section need to be set with constraints:

in the formula, y _f 、V _f 、θ _f Respectively the altitude, the speed and the trajectory inclination angle at the moment of the jumping flight terminal; y is _df 、V _df 、θ _df Respectively the altitude, the speed and the trajectory inclination angle at the moment of the predetermined jump flight terminal; epsilon _h 、ε _v 、ε _θ The upper error bound of the three constraint quantities respectively;

the objective function is a function which enables a certain performance index to be optimal by adjusting the control quantity on the basis that all parameters of the aircraft meet the constraint conditions. According to the specific task selection, the objective function can be generally divided into a maximum flight path (maximum flight path/maximum traverse path), a shortest flight time, a maximum terminal time speed and the like. In the invention, the farthest voyage (farthest voyage/farthest voyage) is taken as a performance index, and then the objective function is as follows:

minJ＝-L _df (10)

wherein L is _df Refers to the aircraft range (longitudinal/lateral).

Based on the above analysis, the high mach number aircraft leap flight segment trajectory optimization problem can be described as: under the conditions of satisfying the mass center motion equation set formula (1), the control constraint formulas (3) to (4), the process constraint formulas (5) to (8) and the terminal constraint formula (9), searching for a control quantity attack angle alpha and a speed inclination angle gamma _v The optimal time series minimizes the objective function (10). A more general optimal control problem is described below: the control variable u (t) is sought so that the Bolza-type performance indicator (i.e., the combination of Mayer-type and Lagrange-type) is minimized

And satisfy the constraint of the mass center motion equation set

Boundary condition constraints

φ(x(t ₀ ),t ₀ ,x(t _f ),t _f )＝0 (13)

And inequality constraints

C(x(t),u(t),t)≤0(14)

Wherein, t ₀ As the starting time, t _f Is the terminal time.

And S3, solving a quasi-optimal solution of a track optimization problem (dynamic optimal control problem) close to a global optimal solution by using an Adaptive Genetic Algorithm (AGA) so as to solve the sensitive defect of the GPM-SQP algorithm to initial values (terminal time and guess values at the initial and final time of a control variable). The Genetic Algorithm (GA) is a random global search optimization method, and the basic principle is to simulate the natural selection mechanism of "superior and inferior" in the biological world, and through population initialization coding, through selection, crossing and variation operations, the genetic algorithm continuously evolves towards the optimal character, and finally generates chromosomes (individuals) meeting the optimization target. The GA algorithm has the advantages of high parallelism, randomness, self-adaptive global optimization probability search and the like, solves a plurality of complex engineering problems (such as multi-objective problem solving), is widely applied to scientific research and engineering practice, and has a specific algorithm flow as shown in FIG. 4. The Adaptive Genetic Algorithm (AGA) can be improved in several aspects such as a fitness function, a cross probability, a mutation probability and the like on the basis of a GA algorithm. Considering that the selection of the fitness function in the GA algorithm is directly related to the convergence rate and whether an optimal solution can be found, the AGA algorithm designed by the present patent is mainly improved as follows for the fitness function:

1) Fitness function adaptive adjustment

In order to increase the population diversity and further quickly and accurately approach the global optimal solution, the fitness function is designed into a form of self-adaptive adjustment along with the increase of the population iteration times. In the initial stage of the algorithm, the individual fitness value difference in the population is large, and in order to prevent the loss of individuals with poor fitness in the initial stage, the fitness function can be kept by changing the fitness function; when the algorithm is close to convergence, the individual fitness value difference in the population is small, and in order to accelerate the convergence speed, the individual fitness value difference can be enlarged by adjusting a fitness function.

2) Logarithmically transforming the fitness function

Considering that the parameter magnitude difference between the performance index and the terminal constraint condition is large, and the parameter magnitude difference cannot be directly used as a fitness function, and transformation such as linear transformation, power law transformation, exponential transformation and the like is required. The patent provides a fitness function based on an inverse function form (namely logarithmic transformation) of exponential transformation, so that the fitness function can better reflect the quality of an individual.

3) Introducing a Multi-object hierarchical Programming method (Lexicgraphing structured Programming, LSP)

What is more difficult to deal with in the fitness function is that the terminal constraint is satisfied and the performance index is optimized at the terminal moment, and a multi-target hierarchical planning method (LSP) can be introduced, namely, a plurality of targets are divided into a plurality of priority levels according to the importance degree. In the optimization of the track of the jump flight segment, the optimal performance index is under the premise of meeting the terminal constraint, namely the importance degree of meeting the terminal constraint is greater than the optimal performance index.

The AGA algorithm optimization process is as follows: firstly, constructing a control quantity data set in a control constraint range, determining a control quantity (an attack angle and a speed inclination angle) segmentation time point, and smoothing the control quantity-time course by using a cubic spline interpolation value; secondly, constructing penalty function processing process constraints (dynamic pressure, overload, heat flow rate, no-fly zone and the like) and terminal constraints, and establishing an objective function by combining a track optimization objective (farthest longitudinal/farthest transverse); and finally, optimizing the control quantity by means of an adaptive genetic algorithm to obtain a track (farthest longitudinal range/farthest transverse range) which meets the constraint condition and has the farthest flight. It should be noted that, in the unpowered jump flight process, the maximum vertical range of the jump flight is closely related to the maximum lift-drag ratio, and on the premise that various constraints and other factors are not considered to be the same, generally speaking, the larger the lift-drag ratio is, the larger the farthest vertical range is. Therefore, when the control quantity segmenting time point is determined, the attack angle corresponding to the maximum lift-drag ratio can be used as a reference quantity to be segmented. S3, firstly, running generations by utilizing an AGA algorithm (the generations are far smaller than generations running when the adaptive genetic algorithm is used independently) to obtain a quasi-optimal solution close to the global optimal solution. The algorithm flow is shown in fig. 5, and specifically includes:

s31, population initialization coding: representing a feasible solution of the problem to be solved into a chromosome or an individual in a genetic algorithm space by adopting a real number coding mode;

s32, calculating the fitness value of each individual in the population through a fitness function, wherein the fitness function is an index for evaluating the chromosome and is also a key point for embodying the 'excellence and disadvantage' of the AGA algorithm. The patent constructs a comprehensive fitness function by using a join method: namely, the performance index of the terminal at the moment and the terminal constraint are weighted to obtain a fitness function. The join method is improved on the basis of a static penalty function, and the penalty coefficient is set to be dynamically and adaptively changed, namely, the penalty coefficient can be adaptively adjusted along with the increase of the iteration number. The concrete form is as follows:

in the formula: r is ₁ ＝(cd) ^η ，f _i (x) The specific form of the penalty function term for violating the ith terminal constraint is as follows:

in the formula: p is a radical of formula _i And constraining the weight coefficient of the penalty term for the ith terminal.

And (4) processing the remaining constraint conditions:

1) And (4) controlling and constraining: constructing a control quantity data set in the constraint condition range to automatically meet the control constraint;

2) Process constraints (dynamic pressure, overload, heat flow rate and no-fly zone): will not satisfy the maximum dynamic pressure q _dmax Maximum overload n _Lmax Maximum heat flow rate

And the fitness function value for crossing the no-fly zone is set to infinity, i.e. when F = inf.

In the formula: c. eta, eta,

All are constants, and are combined with relevant documents and verified by numerical simulation, and c =0.5 is selected and is selected as the value of the blood pressure in the blood vessel>

Is a feasible scheme; d is iteration times, J is a track optimization objective function term, q is a corresponding weight coefficient, and K is a constant (ensuring that the fitness is positive). The terminal state of the jumping flight segment is strictly restricted andthe performance indexes of the terminal at the moment are required to meet the optimal requirement, wherein the optimal performance indexes are that on the premise of meeting terminal constraints, on the basis of LSP and combined with related documents, the difference between terminal constraint items at the same level is small through logarithmic processing, the weight coefficients of penalty items of terminal constraints can be regarded as having the same weight, namely the same constant value, and on the basis, the weight coefficient q of the performance index at the next level is proportionally adjusted (namely the weight coefficient q needs to be tested to ensure that the terminal constraint items are in a track range).

S33, judging whether the maximum iteration times are reached, if so, obtaining a quasi-optimal solution, and entering S4; otherwise, the process goes to S34;

s34, carrying out genetic operations of selection, crossing and mutation on individuals according to the fitness value to generate a new generation of population, and entering S32;

selecting genetic operation: selecting excellent individuals from the population according to a certain probability to form a new population to breed to obtain next generation individuals, wherein the higher the individual fitness value is, the higher the probability of selection is. This patent selects the roulette method, and the probability that individual i was chosen is:

cross genetic manipulation: the cross operation is the most important operation of the AGA algorithm and is directly related to the global searching capability of the algorithm. The core idea is to randomly select two individuals from a group, and cross-combine the two individuals (chromosomes), so that the excellent characteristics of the excellent individuals of the parent are inherited to the offspring, and a new excellent individual is generated. The method adopts a real number intersection method, namely the kth chromosome a _k And the l-th chromosome a _l The interleaving operation at j bits is:

variant genetic manipulation: in order to maintain population diversity, the premature convergence phenomenon is prevented by selecting a proper variation rate, and the local optimum is avoided.The mutation operation is to change the value of one gene of the individual string in the population, and the jth gene a of the ith individual _ij The method for performing mutation operation comprises the following steps:

wherein, F _i Is the fitness value of an individual i, Q is the number of population individuals, b is [0,1]Random number of (a) _max 、a _min Are respectively a gene a _ij Upper and lower bounds of (a); f (g) = r ₂ (1-g ₀ /G _max ) ² ，r ₂ Is a random number, g ₀ Is the current iteration number, G _max For maximum evolutionary number, r is [0,1 ]]The random number of (2).

And S4, taking the quasi-optimal solution as an initial value (including a terminal time and a control quantity initial and final time value) of the GPM-SQP algorithm, and optimizing the control quantity in a guessing time period by using the GPM-SQP algorithm to obtain a global optimal solution or a final solution close to the global optimal solution.

The optimization process of the GPM-SQP algorithm is as follows:

firstly, converting an original optimal control problem into NLP by utilizing a GPM algorithm, and converting a differential equation into an algebraic equation; and then dividing the nonlinear constraint optimization problem into a series of quadratic programming subproblems by utilizing an SQP algorithm to solve. And S4, mainly taking terminal time and control quantity initial and final time values obtained by solving with the AGA algorithm as initial values of the GPM-SQP algorithm, so as to find out a global optimal solution or a final solution close to the global optimal solution, wherein the algorithm flow is shown in figure 6. The Gauss-SQP algorithm comprises two steps of variable dispersion and parameterization and numerical optimization method, and specifically comprises the following steps:

s41, variable discrete and parameterization-Gauss pseudo-spectral method (GPM)

The core idea of Gauss pseudo-spectrum method is as follows: dispersing state variables and control variables on a series of Legendre-Gauss (LG) nodes at the same time, approximating the state variables and the control variables by adopting global polynomial interpolation, converting system differential constraints into algebraic constraints, and then taking the state variables and the control variables at the nodes as optimization design variables so as to convert a track optimization problem (dynamic optimal control problem) into NLP (non-linear programming) for solving, wherein the method comprises the following specific steps of:

(1) time domain conversion: setting the initial time and the terminal time of the dynamic optimal control problem as t ₀ 、t _f And the discrete point of GPM is selected to be [ -1,1]Therefore, the following variables need to be replaced when solving: expressed as:

(2) discrete state variables (state variable approximation): coordination point tau of GPM _k (k =1, 2.. Multidot.n) is N LG points with a distribution interval of (-1, 1) points, the root (τ) of a Legendre polynomial of order N ₁ ,τ ₂ ,...,τ _k ) And τ ₀ = -1 as the configuration node of the system, construct N +1 Lagrange polynomial approximation state variables, expressed as:

wherein the discretized state value at the configuration node is exactly equal to the actual state value of the original system, namely x (tau) _i )＝X(τ _i )。

(3) Discrete control variables (control variable approximation): root of Legendre polynomial of N order (tau) ₁ ,τ ₂ ,...,τ _k ) As configuration nodes of the system, N Lagrange polynomial approximate control variables are constructed and expressed as:

(4) converting a differential equation at the distribution point into an algebraic equation, specifically:

for the approximate state variables, the derivation with respect to τ is made:

and the derivative of each Lagrange polynomial at the LG point can be represented by a differential approximation matrix:

in the formula:

then:

therefore, the differential equation constraint corresponding to the motion equation set of the aircraft at the distribution point can be converted into an algebraic constraint. Namely:

(5) discretizing the constraint condition;

discretizing control constraints: in GPM, the control quantity only needs to satisfy the constraint condition at the LG point, i.e.:

process constraint discretization: for process constraints, it can be considered that as long as the constraints are satisfied at the discrete nodes, the process constraints are satisfied throughout, that is:

C(X(τ _k ),U(τ _k ),τ _k ；t ₀ ,t _f )≤0,k＝1,2,...,N (28)

and (3) discretizing terminal state constraint: the boundary condition (initial state and terminal transition) constraints can be described as follows:

φ(X(t ₀ ),t ₀ ,X(t _f ),t _f )＝0 (29)

the initial state is fixed, knowing x (τ) ₀ )≡X(τ ₀ ) Is equal to X (-1), the terminal time tau is not taken into account due to the approximation of the state variable by equation (21) _f Therefore, the terminal state needs to be supplemented, and the following formula (26) can be obtained:

/>

in addition, the terminal time state variable X (τ) of the system _f ) Can be composed of X _k (k =0,1,2,. Ang., N) and U _k (k =1, 2.., N) is obtained by gaussian integration:

in the formula:

is Gaussian weight, wherein->

Then: />

τ _k Is the LG point.

(6) And discretizing the farthest voyage performance index function.

The Bolza-type performance indicator function:

the integral term in the equation is approximated by gaussian integration, and the performance index function in GPM can be obtained:

thus, by:

the original dynamic optimal control problem can be converted into a nonlinear programming problem (NLP), and the converted nonlinear programming problem and the original control problem have consistent solutions.

In summary, the trajectory optimization problem can be described as: at [ t ] ₀ ,t _f ]Determining the state variable X (tau) at discrete points in time _k ) (k =0,1, 2.. Multidot., N), control variable U (τ) _k ) (k =1, 2.., N), initial time t ₀ Terminal time t _f And minimizing the performance index (33) under the condition that the constraint conditions (the transformed aircraft motion algebraic equation system constraint (26), the control constraint (27), the process constraint (28) and the terminal state constraint (31)) are met.

S42, numerical optimization method-Sequence Quadratic Programming (SQP)

Aiming at the discretized nonlinear programming problem (NLP), the sequential quadratic programming is a method for directly processing the constraint optimization problem, and the basic idea is as follows: in each iteration step, a quadratic programming subproblem is solved to determine a descending direction, so as to reduce the cost function to obtain the step length, and the steps are repeated until the solution of the original problem is obtained.

Conclusion

Aiming at the problem of trajectory optimization of a high-Mach-number aircraft under complex constraint conditions including multiple flight forbidden regions and the like in a jumping flight section, a solution based on a GPM-AGA-SQP hybrid algorithm is provided.

1) In order to further rapidly and accurately approach a quasi-optimal solution of the global optimal solution, three-point improvement is carried out on a fitness function in the genetic algorithm, a logarithmic comprehensive fitness function which can be adaptively adjusted along with the increase of population iteration times is constructed, and a new adaptive genetic algorithm is established;

2) Aiming at the defect that the GPM-SQP algorithm excessively depends on the initial guess value when solving the track optimization problem, the control quantity is globally optimized based on the established adaptive genetic algorithm, and the obtained approximate global optimal solution is used as the initial guess value of the GPM-SQP algorithm, so that the sensitivity of the GPM-SQP algorithm to the initial guess value is reduced;

3) The constructed GPM-AGA-SQP hybrid algorithm solves the problems of low convergence speed and weak local optimization capability of the AGA algorithm when the AGA algorithm approaches to a global optimal solution, and simultaneously combines the advantage of strong robustness of the AGA algorithm, so that the hybrid algorithm can still complete a post-track optimization task under the condition of model uncertainty and disturbance interference.

Experiments show that the established hybrid optimization algorithm: the GPM-AGA-SQP algorithm has effectiveness and feasibility, can realize avoidance of multiple flight-forbidden zones, is improved by more than 16% compared with the performance index value of the improved direct targeting method-AGA algorithm, is smoother in control quantity change curve and smaller in buffeting amplitude compared with the GPM-AACO-SQP algorithm, is more consistent with the actual situation of avoidance of the flight-forbidden zones, is closest to the global optimal solution, is improved by more than 2.05%, and meanwhile, the number and difficulty of parameter adjustment are reduced, so that the ideal track meeting the constraint condition and the farthest voyage is easily obtained; in addition, model uncertainty and disturbance factors in the actual flight process are considered, and Monte Carlo simulation tests show that various constraint conditions in the jumping flight process can still meet requirements, so that the constructed GPM-AGA-SQP algorithm has certain robustness; in the aspect of application, the GPM-AGA-SQP algorithm is essentially an optimization method, is not only suitable for optimizing the track of a jumping flight section of a high-Mach-number aircraft, but also suitable for optimizing the track of a general aircraft, and provides a certain reference value for the next step of researching the uncertainty quantized robust track optimization problem.

Simulation (Emulation)

According to the invention, a large number of simulations are carried out on the basis of the GPM-AGA-SQP hybrid algorithm, so that the effectiveness and the feasibility of the constructed hybrid algorithm are verified, and the comparison is carried out with the GPM-AACO-SQP algorithm and the simulation result of the improved direct targeting method-AGA algorithm so as to verify the superiority of the GPM-AGA-SQP algorithm, and finally, the Monte Carlo simulation test is utilized to verify the robustness of the GPM-AGA-SQP algorithm.

(1) Simulation related parameters

Before the track optimization simulation of the jumping flight segment of the high-Mach number aircraft is carried out, basic parameters of the aircraft and the earth, control constraint, process constraint, boundary constraint, established new AGA algorithm related parameters and guessing values of terminal time and initial and final time of control quantity obtained after the AGA algorithm finishes control quantity optimization are required to be given.

Aircraft and earth basic parameters are shown in table 1, control constraints are shown in table 2, process constraints are shown in table 3, and boundary constraints (initial and terminal states) are shown in table 4:

TABLE 1 simulation basic parameters

TABLE 2 control constraints

TABLE 3 Process constraints

TABLE 4 boundary constraints (initial and terminal states)

Note: speed 10Ma ≈ 3297.99m · s- ¹ The speed 3Ma ≈ 899.40m · s corresponding to the height 27km of the terminal state ^-1 。

The main parameters of the adaptive genetic algorithm are shown in table 5:

TABLE 5 Main parameters of adaptive genetic Algorithm

The AGA algorithm designed by the patent is mainly aimed at the improvement of a GA algorithm fitness function, so that other 4 parameters need to be set in advance. These 4 parameters are generally set within the following ranges: population scale: 20 to 100; iteration times are as follows: 100 to 500; the cross probability: 0.4 to 0.99; mutation probability: 0.0001 to 0.1. Considering that in the actual simulation process, the AGA algorithm may have a premature phenomenon, some more optimal solutions are omitted, the problem can be solved by increasing the population number and increasing the variation rate, and meanwhile, in order to avoid too large calculation amount, the population size can be selected to be 500 (5 times expansion is performed on the original basis, test is needed), the number of iterations is selected to be 100 (the minimum number of iterations is selected to obtain a quasi-optimal solution), and the variation probability is selected to be 0.1 (the maximum variation probability). In addition, the cross probability is selected to be 0.8 +/-0.5, so that the solving accuracy and speed can be improved, and the problem that a better solution cannot be obtained due to improper cross probability selection is avoided, so that the cross probability can be selected to be 0.8 (taking a median).

The optimum fitness of the AGA algorithm is shown in fig. 7, and fig. 7 shows the process that the optimum fitness value in the population changes with the number of iterations, and it can be known from the figure that the optimum fitness value in the population changes greatly in the initial stage, gradually converges with the increase of the number of iterations, and tends to be stable already in about 10 generations, and slightly changes in the later stage in about 60-90 generations, but the change is small. This shows that the AGA algorithm has a good convergence effect on the optimization of the control variable.

After the Adaptive Genetic Algorithm (AGA) finishes optimizing the controlled variable, the time-varying curves of the attack angle and the speed inclination angle are shown in fig. 8-9, and fig. 8-9 are curves obtained by optimizing the function values of the control variable segment time points by using the adaptive genetic algorithm and performing cubic spline interpolation smoothing on the control variable-time history under the conditions that a controlled variable data set constructed in a control constraint range, a terminal constraint is used as a judgment condition, the farthest longitudinal or farthest transverse stroke is used as a performance index and process constraints (dynamic pressure, overload, heat flow rate, no-fly zone and the like) are considered.

It can be seen from fig. 8-9 that the farthest longitudinal is used as the performance index, the time history of the control quantity optimized by the adaptive genetic algorithm is 495s, the initial and final time values of the attack angle of the control quantity are 0.175rad and 0.175rad, respectively, and the initial and final time values of the speed inclination angle of the control quantity are-0.2975 rad and-0.1575 rad, respectively, which are used as the initial values solved by the GPM-SQP algorithm with the farthest longitudinal being used as the performance index.

It can be seen from fig. 8-9 that the farthest transverse distance is used as the performance index, the time history of the control quantity optimized by the adaptive genetic algorithm is 475s, the initial and final time values of the attack angle of the control quantity are 0.14rad and 0.1925rad, respectively, the initial and final time values of the speed inclination angle of the control quantity are 0.4025rad and 0.1575rad, respectively, and the initial values are obtained by using the farthest transverse distance as the performance index and using the GPM-SQP algorithm to solve. (2) Validity and feasibility check of GPM-AGA-SQP algorithm

According to the invention, the farthest longitudinal range/farthest transverse range is taken as a performance index, and a large amount of simulation is carried out based on a GPM-AGA-SQP hybrid algorithm so as to verify the effectiveness and feasibility of the hybrid algorithm.

1) Farthest longitudinal trajectory optimization

In the track optimization problem with the farthest longitudinal voyage as the target, the performance index is as follows if the longitudinal voyage corresponds to the change of x:

J ₀ ＝min(-L ₀ )＝min(-x(t _f )) (35)

when the farthest longitudinal track optimization is carried out by using a GPM-AGA-SQP algorithm, a GPOPS tool box is adopted and SNOPT software is combined to carry out solving, and simulation results are shown in figures 10-24.

2) Farthest traverse trajectory optimization

In the track optimization problem with the maximum lateral voyage as the target, the lateral voyage corresponds to the change of z, and the performance indexes are as follows:

J ₀ ＝min(-L ₀ )＝min(-|z(t _f )|) (36)

and when the maximum longitudinal trajectory optimization is carried out by utilizing a GPM-AGA-SQP algorithm, a GPOPS tool box is adopted and SNOPT software is combined to solve, and simulation results are shown in figures 25-34.

From fig. 16-18, 31-33, it can be seen that the hybrid algorithm is constructed: the GPM-AGA-SQP algorithm controls and restricts in a specified range in the simulation taking the farthest longitudinal travel/the farthest transverse travel as a performance index; in FIGS. 20-22, 35-37, the process constraints (dynamic pressure, overload, heat flow rate) are also within the constraints; fig. 23-24 and 38-39 show that the hybrid algorithm can avoid multiple no-fly zones, and meet actual requirements.

The GPM-AGA-SQP hybrid algorithm constructed for a more intuitive reaction meets the terminal constraint, and the following table 6 is established for verification.

Table 6 terminal states corresponding to terminal constraint conditions of GPM-AGA-SQP algorithm

From fig. 10-11, 14 and 25-26, 29 and table 6, it can be seen that the hybrid algorithm can satisfy the terminal constraints in the simulation with the farthest longitudinal/farthest transverse as the performance index. In conclusion, the effectiveness and feasibility of the constructed GPM-AGA-SQP hybrid algorithm are demonstrated.

(3) GPM-AGA-SQP algorithm superiority checking

The method takes the farthest longitudinal range/farthest transverse range as a performance index, and carries out comparison simulation based on a GPM-AGA-SQP algorithm, a GPM-AACO-SQP algorithm and an improved direct targeting method-AGA algorithm so as to verify the superiority of the GPM-AGA-SQP algorithm.

It should be noted that, an Adaptive Ant Colony Optimization (AACO) is improved based on the ACO algorithm, and the global search/local search step size and the penalty coefficient in the fitness function are set to be in a form that can be adaptively adjusted along with the increase of the number of iterations. The first stage of the GPM-AACO-SQP algorithm is to obtain a quasi-optimal solution close to a global optimal solution by utilizing the AACO algorithm, and the core idea of the second stage is similar to that of the GPM-AGA-SQP algorithm. The improved direct targeting method-AGA algorithm is to optimize and improve the control quantity at discrete points in the direct targeting method through the AGA algorithm, and numerically integrate the motion equation set based on the parameterized control quantity to obtain a state variable describing a motion track. The improved direct target shooting method is provided for solving the problems that the direct target shooting method has high requirement on an initial value and is likely to generate local minimum, the initial and final time values of the controlled variable are brought into an optimization design variable, a composite constraint item (such as speed, height and the like) in a terminal constraint is taken as a judgment condition, and the track ending time is taken as a terminal time, so that the purpose of reducing sensitivity to the initial value is achieved.

1) Comparison of Performance indicators

The results of the performance index simulation comparison of the three algorithms are shown in fig. 40-41 and table 7:

TABLE 7 comparison of the optimal performance indexes of the three algorithms

From the graphs 40-41 and the table 7, it can be seen that, in the variation trend of the performance index, the GPM-AGA-SQP algorithm and the GPM-AACO-SQP algorithm are relatively close to each other in the variation trend, and compared with the improved direct targeting method-AGA algorithm, the variation trend in the early stage is basically consistent, but the difference in the later stage trend is large, which explains to a certain extent the defect of poor local search capability when the improved direct targeting method-AGA algorithm approaches the global optimal solution; on the performance index value, the maximum longitudinal value of the performance index of the GPM-AGA-SQP algorithm is improved by 2.05 percent compared with that of the GPM-AACO-SQP algorithm, the maximum transverse value is improved by 5.64 percent, and the maximum longitudinal value of the performance index of the GPM-AGA-SQP algorithm is improved by 16.82 percent and the maximum transverse value is improved by 18.63 percent compared with that of the performance index of the GPM-AACO-SQP algorithm.

2) Comparison of control quantity optimizing ability

The three algorithm control quantity variation curves are shown in figures 42-45. From the graphs 42 to 45, it can be seen that the GPM-AGA-SQP algorithm and the GPM-AACO-SQP algorithm are basically consistent in the control quantity change trend, but compared with the GPM-AACO-SQP algorithm, the control quantity change curve corresponding to the GPM-AGA-SQP algorithm is smoother, the buffeting amplitude is smaller, the actual situation of avoiding a no-fly zone is better met, and the global optimal solution is closest to; compared with the GPM-AGA-SQP algorithm, the improved direct targeting method-AGA algorithm shows certain inertia on the change trend of the control quantity, namely small fluctuation of the control quantity in the GPM-AGA-SQP algorithm is not as good as the change of the control quantity of the improved direct targeting method-AGA algorithm or the fluctuation range is larger, the optimal searching capability of the improved direct targeting method-AGA algorithm on the control quantity is reflected to a certain extent, namely the optimal searching capability on the large-range global searching capability is stronger, but the local searching capability is weaker. It should be noted that the fluctuation of the control quantity is a fluctuation around the "quasi-optimal solution" of the global optimal solution, rather than a random fluctuation, which is also a reason why the terminal time value obtained by the AGA algorithm and the initial and final time values of the control quantity can be used as the initial guess values of the GPM-SQP algorithm.

3) Avoidance capability of no-fly zone and test parameter comparison

Three algorithms avoid the top view of the no-fly zone as shown in fig. 46-47. From FIGS. 46-47, it can be seen that the GPM-AGA-SQP algorithm, the GPM-AACO-SQP algorithm, and the modified direct targeting-AGA algorithm are all able to achieve avoidance of multiple no-fly zones. It should be noted that the GPM-AGA-SQP algorithm and the GPM-AACO-SQP algorithm provide a control quantity at the beginning and end time and a guess value at the end time through the AGA algorithm and the AACO algorithm, respectively, and the test parameters are mainly concentrated in the AGA and AACO algorithm parts. The AGA algorithm actually needs 2 test parameters, namely a performance index weight coefficient Q and a population scale Q; the AACO algorithm has 6 actual test parameters, and compared with the AGA algorithm, the pheromone evaporation coefficient Rho, the global transfer selection factor P0 and the search step length adjusting factor step are added (the search step length is not out of range, and the step length adjustment is enabled to be in a certain precision range). The following compares the three algorithm test parameters, as shown in table 8 below:

TABLE 8 three Algorithm test parameters

From table 8, it can be seen that, the constructed hybrid algorithm: the number of the GPM-AGA-SQP algorithm tuning parameters is 3, wherein the terminal time value obtained based on the AGA algorithm is finely tuned (the fine tuning range is about +/-3 s through simulation verification). Compared with the GPM-AACO-SQP algorithm, the number of the parameters of the GPM-AGA-SQP algorithm (the number of the parameters of the GPM-AACO-SQP algorithm is 6) and the difficulty (the fine tuning range of the terminal time value obtained based on the AACO algorithm is about +/-5 s after simulation verification) are reduced, and the track which meets the constraint condition and has the optimal performance index (the farthest longitudinal range/the farthest transverse range) is easy to obtain.

(4) GPM-AGA-SQP algorithm robustness check

The influence of model uncertainty and disturbance factors on the track optimization performance of the high-Mach aircraft needs to be considered in the jumping flight process, and the robustness of the high-Mach aircraft is verified by carrying out 50 Monte Carlo simulation tests through a GPM-AGA-SQP hybrid algorithm. Wherein, in the aspect of model uncertainty, the pneumatic parameter C is influenced by factors such as manufacturing appearance, measurement error and the like _x And C _y The uncertainty is more reasonable to consider. The uncertainty amount and perturbation are shown in table 9.

TABLE 9 Monte Carlo simulation uncertainty factors and disturbance parameters

Under the influence of the uncertainty parameters and the disturbance, the average value, the maximum value and the minimum value of the performance indexes, the process constraint items (including heat flow rate, dynamic pressure and overload), the terminal constraint items (including terminal speed, terminal height and terminal trajectory inclination angle) parameters obtained by the GPM-AGA-SQP hybrid algorithm are shown in tables 10-11, and Monte Carlo simulation tracks are shown in FIGS. 48-49.

TABLE 10 Monte Carlo simulation results (performance index: farthest longitudinal)

TABLE 11 Monte Carlo simulation results (Performance index: farthest horizontal run)

From tables 10-11, FIGS. 48-49, it can be derived that the hybrid algorithm was constructed: the GPM-AGA-SQP algorithm has certain robustness under the influence of model uncertainty and disturbance factors, and in simulation taking the farthest longitudinal stroke/the farthest transverse stroke as performance indexes, all control constraints, terminal constraints and process constraints (dynamic pressure, overload, heat flow rate and no-fly zone) basically meet conditions.

The technical solution of the present invention is not limited to the limitations of the above specific embodiments, and all technical modifications made according to the technical solution of the present invention fall within the protection scope of the present invention.

Claims

1. A method for optimizing a track of a jumping flight segment of a high-Mach aircraft is characterized by comprising the following steps:

2. The method for optimizing the jumping flight segment trajectory of the high-mach-number aircraft according to claim 1, the method is characterized in that the simplified model of the mass center motion equation set is as follows:

in the formula: v (t) is the aircraft speed, theta (t) and psi _v (t) ballistic dip and ballistic declination, respectively; y (t) is height, x (t) and z (t) are respectivelyA longitudinal course and a transverse course; m (t) is the aircraft mass, g is the gravitational acceleration, α (t) is the angle of attack, γ _v (t) is the velocity ramp angle; alpha is alpha _B To balance angle of attack, beta _B To balance the sideslip angle; y is _B And Z _B Are respectively alpha _B 、β _B Corresponding balance lift force and balance lateral force, X is resistance, and the calculation expression is as follows:

in the formula: c _x 、C _y Q =0.5 ρ V for drag coefficient and lift coefficient ² And rho is the air density of the height of the high Mach number aircraft, and S is the reference area of the aircraft.

3. The method for optimizing the jumping flight segment trajectory of the high mach number aircraft according to claim 2, wherein in S2, the constraint conditions include self control quantity constraint, process constraint applied during the jumping flight segment flight process, and terminal constraint on the destination, and the objective function selects the farthest flight path, which is as follows:

the control quantity constraints include an angle of attack alpha and a speed tilt angle gamma _v The constraint is expressed as:

Certain limitations also need to be imposed, namely:

at this time, the attack angle change rate can be adjusted

The attack angle alpha of the original control quantity is regarded as a state variable for processing;

the process constraints comprise heat flow rate constraints, dynamic pressure constraints, overload constraints and no-fly zone constraints, and are specifically as follows:

the heat flow rate constraint is expressed as:

the dynamic pressure constraint is expressed as:

q＝0.5*ρ*V ² ≤q _dmax (6)

the overload constraint is expressed as:

the type of a no-fly zone avoided by the high-Mach number aircraft in the jumping flight process is cylindrical, and the central point is (x) _i ,y _i ,z _i ) Radius r _i The specific expression is as follows:

l (x (t), t) is the shortest distance from the high-Mach aircraft to the no-fly zone, epsilon _n For infinitesimal normal quantities, the no-fly zone constraint can be expressed as L (x (t), t) ≧ ε _n In the formula, K _Q Is a constant parameter associated with high mach number aircraft,

to the maximum allowable heat flow rate, q _dmax For maximum allowable dynamic pressure, n _Lmax Maximum overload allowed;

the terminal constraints include altitude, velocity and ballistic dip constraints at the terminal time, expressed as:

in the formula, y _f 、V _f 、θ _f Respectively the altitude, the speed and the trajectory inclination angle at the moment of the jumping flight terminal; y is _df 、V _df 、θ _df Respectively the height, the speed and the trajectory inclination angle of a preset jump flight terminal moment; epsilon _h 、ε _v 、ε _θ The upper error bound of the three constraint quantities respectively;

the objective function is expressed as:

minJ＝-L _df (10)

the high mach number aircraft leap flight segment trajectory optimization problem can be described as: under the conditions of satisfying the mass center motion equation set formula (1), the control constraint formulas (3) to (4), the process constraint formulas (5) to (8) and the terminal constraint formula (9), searching for a control quantity attack angle alpha and a speed inclination angle gamma _v The optimal time sequence minimizes the objective function (10); a more general optimal control problem is described below: the control variable u (t) is sought so that the Bolza type performance indicator: combination of Mayer type and Lagrange type, minimum

And satisfy the constraint of the mass center motion equation set

Boundary condition constraints

φ(x(t ₀ ),t ₀ ,x(t _f ),t _f )＝0 (13)

And inequality constraints

C(x(t),u(t),t)≤0 (14)

Wherein, t ₀ Is a starting time, t _f Is the terminal time.

4. The method for optimizing the trajectory of the jumping flight segment of the high-mach-number aircraft according to claim 3, wherein in S3, a quasi-optimal solution of the trajectory optimization problem close to a global optimal solution is obtained by using an adaptive genetic algorithm AGA to solve the defect that a GPM-SQP algorithm is sensitive to an initial value, and specifically comprises:

s32, constructing a comprehensive fitness function by using a join method: the method comprises the steps that a fitness function is obtained by weighting a performance index of a terminal at a moment and terminal constraint, and the fitness value of each individual in a population is calculated through the fitness function; the fitness function is an index for evaluating the chromosome and is a key point for reflecting the excellence and the weakness of the AGA algorithm; the join method is improved on the basis of a static penalty function, and the penalty coefficient is set to be dynamically and adaptively changed, namely, the penalty coefficient can be adaptively adjusted along with the increase of iteration times, and the specific form is as follows:

in the formula: r ₁ ＝(cd) ^η ，f _i (x) The specific form of the penalty function term for violating the ith terminal constraint is as follows:

in the formula: p is a radical of _i And (3) constraining the weight coefficient of the penalty term for the ith terminal, and processing the rest constraint conditions:

1) And (3) controlling and constraining: constructing a control quantity data set in the constraint condition range to automatically meet the control constraint;

2) And (3) process constraint: will not satisfy the maximum dynamic pressure q _dmax Maximum overload n _Lmax Maximum heat flow rate

And the fitness function value of the no-fly zone is set to be infinite, namely F = inf;

in the formula: c. eta, b,

Are all constant, and are verified by numerical simulation, and c =0.5 is taken>

Is a feasible scheme; d is iteration times, J is a track optimization objective function term, q is a corresponding weight coefficient, K is a constant and the fitness is guaranteed to be positive; the terminal state of the jumping flight segment is strictly restricted and the performance index of the terminal at the moment is required to meet the optimal performance index, wherein the optimal performance index is based on LSP on the premise of meeting the terminal restriction, the difference between terminal restriction items in the same level is smaller through logarithmic processing, the weight coefficients of penalty items of the terminal restriction can be regarded as having the same weight and being the same fixed value, and on the basis, the weight coefficient q of the next-level performance index is proportionally adjusted, namely the weight coefficient q needs to be tested to ensure that the terminal restriction items are in a track range;

s33, judging whether the maximum iteration times is reached, if so, obtaining a quasi-optimal solution, and entering S4; otherwise, the process goes to S34;

selecting genetic operation: selecting excellent individuals from the population according to a certain probability to form a new population to breed to obtain next generation individuals, wherein the higher the individual fitness value is, the higher the selected probability is; this patent selects the roulette method, and the probability that individual i was chosen is:

cross genetic manipulation: the cross operation is the most important operation of the AGA algorithm and is directly related to the global searching capability of the algorithm, the core idea is that two individuals are randomly selected from a group, and the two individuals or chromosomes are cross-combined, so that the excellent characteristics of parent excellent individuals are inherited to filial generations, and new excellent individuals are generated; the method adopts a real number intersection method, namely the kth chromosome a _k And the l-th chromosome a _l The interleaving operation at j bits is:

variant genetic manipulation: in order to maintain population diversity, the premature convergence phenomenon is prevented by selecting a proper variation rate, and the local optimum is avoided; the mutation operation is to change the value of one gene of the individual string in the population, and the jth gene a of the ith individual _ij The method for performing mutation operation comprises the following steps:

5. The trajectory optimization method for the jumping flight segment of the high-mach-number aircraft according to claim 5, characterized in that in S4, the quasi-optimal solution is used as an initial value of a GPM-SQP algorithm, and then the GPM-SQP algorithm is used to optimize the controlled variable within a guessing time period to obtain a global optimal solution or a final solution close to the global optimal solution; the Gauss-SQP algorithm comprises two steps of variable dispersion and parameterization and numerical optimization method, and specifically comprises the following steps:

s41, variable discretization and parameterization-Gauss pseudo-spectral method, i.e. GPM

(2) discrete state variables, i.e. state variable approximations: collocation point tau of GPM _k (k =1, 2.. Multidot., N) is N LG points, the distribution interval of the collocation points is (-1, 1), and the roots (τ) of the Legendre polynomial of order N are divided into ₁ ,τ ₂ ,...,τ _k ) And τ ₀ = -1 as the configuration node of the system, construct N +1 Lagrange polynomial approximation state variables, expressed as:

wherein the discretized state value at the configuration node is exactly equal to the actual state value of the original system, namely x (tau) _i )＝X(τ _i )；

(3) Discrete controlled variables, i.e. controlled variable approximations: root of Legendre polynomial of N order (tau) ₁ ,τ ₂ ,...,τ _k ) As configuration nodes of the system, N Lagrange polynomial approximate control variables are constructed and expressed as:

the derivation is made with respect to τ for the approximate state variables:

/>

in the formula:

then:

therefore, the differential equation constraint corresponding to the motion equation system of the aircraft at the distribution point can be converted into an algebraic constraint, namely:

(4) discretizing the constraint condition specifically as follows:

C(X(τ _k ),U(τ _k ),τ _k ；t ₀ ,t _f )≤0,k＝1,2,...,N(28)

and (3) discretization of terminal state constraint: the boundary conditions, i.e., initial state and terminal transition constraints, can be described as follows:

φ(X(t ₀ ),t ₀ ,X(t _f ),t _f )＝0 (29)

the initial state is fixed, knowing x (τ) ₀ )≡X(τ ₀ ) Is equal to X (-1), the terminal time tau is not taken into account due to the approximation of the state variable by equation (21) _f Therefore, the terminal status needs to be supplemented, and the following formula (26) can be obtained:

in addition, the terminal time state variable of the system

Can be composed of X _k (k =0,1,2,. Ang., N) and U _k (k =1, 2.., N) is obtained by gaussian integration:

in the formula:

is Gaussian weight, wherein>

Then: />

τ _k Is LG point;

(5) discretizing the farthest voyage performance index function, specifically:

the Bolza-type performance indicator function:

thus, by:

the original dynamic optimal control problem can be converted into a nonlinear programming problem, and the converted nonlinear programming problem has a consistent solution with the original control problem;

in summary, the trajectory optimization problem can be described as: at [ t ] ₀ ,t _f ]Determining the state variable X (tau) at discrete points in time _k ) (k =0,1, 2.. Multidot., N), control variable U (τ) _k ) (k =1, 2.., N), initial time t ₀ Terminal time t _f When the constraint conditions are met, namely the transformed aircraft motion algebraic equation system constraint (26), the control constraint (27), the process constraint (28) and the terminal state constraint (31), the performance index (33) is minimum;

s42, numerical optimization method-sequence quadratic programming, namely SQP algorithm

Aiming at the problem of discretized nonlinear programming, sequential quadratic programming is a method for directly processing a constraint optimization problem, and the basic idea is as follows: in each iteration, a descending direction is determined by solving a quadratic programming subproblem to reduce the cost function to obtain the step size, and the steps are repeated until the solution of the original problem is obtained.