Improved trajectory optimization method combining direct shooting method and adaptive genetic algorithm

Abstract

The application discloses a trajectory optimization method combining improved direct shooting method and adaptive genetic algorithm, relates to the field of aircrafts, and comprises the following steps: S1, a simplified model of a center of mass motion equation set is obtained by establishing a dynamic equation and a kinematic equation; S2, a constraint condition and an objective function are determined, and a trajectory optimization problem is described as a standard dynamic optimal control problem; S3, the optimal control problem is converted into a nonlinear programming problem by using the improved direct shooting method; S4, a control variable data set is constructed within the constraint range of the control variable, an initial state X0 is taken as input, initial and final time values of the control variable are taken into optimization design variables, the control variable at discrete points in the improved direct shooting method is optimized by using the adaptive genetic algorithm, and the control variable-time history is smoothed by using a cubic spline interpolation, numerical integration is carried out based on a fourth-order Runge-Kutta method with the monotonous change amount (such as speed) in the terminal constraint as a judgment condition, and the optimal control u * is calculated through several iterations, so that an ideal trajectory meeting the conditions is obtained.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of aircraft trajectory optimization, in particular to a trajectory optimization method combining improved direct shooting method and adaptive genetic algorithm. BACKGROUND

[0002] The trajectory optimization problem of quasi-equilibrium flight segment of high Mach number aircraft is a hot research topic. During the quasi-equilibrium flight of high Mach number aircraft, not only the constraints of dynamic pressure, overload, heat flux rate, etc. need to be considered, but also the avoidance of multiple no-fly zones (formed by natural environment, political and military factors) is required. The existing research results on no-fly zones are relatively shallow, mainly through trajectory optimization or design of lateral guidance logic to avoid no-fly zones.

[0003] In terms of trajectory optimization, the trajectory optimization problem of quasi-equilibrium flight segment of high Mach number aircraft has high nonlinearity and complexity, and it is difficult to efficiently obtain the optimal solution by using classical variational method and Pontryagin minimum principle, so the direct method is usually used for solving. The direct shooting method is a typical representative method in the direct method, which converts the optimal control problem into a nonlinear programming problem by discretizing the control variable, and has the advantages of good universality and easy implementation. However, the convergence performance of the direct shooting method depends on the terminal time and the initial and final time guess value of the control quantity, and in engineering practical application, the initial value is difficult to estimate, so it needs to be improved. In terms of lateral guidance logic, the methods mainly include the feeler method, the sign inversion method of roll angle, the artificial potential field method, etc., which design the guidance strategy online to avoid no-fly zones, but do not consider the optimality of the terminal time performance index (such as the maximum range, etc.). SUMMARY

[0004] The purpose of the present application is to solve the above problems by designing a trajectory optimization method combining improved direct shooting method and adaptive genetic algorithm.

[0005] The present application achieves the above-mentioned purposes through the following technical solutions:

[0006] The trajectory optimization method combining improved direct shooting method and adaptive genetic algorithm comprises:

[0007] S1, ignoring the effects of earth rotation and curvature, the earth is assumed to be a plane, and the dynamic equation and kinematic equation are established in the trajectory coordinate system and the ground coordinate system respectively, obtaining a 3-DOF high Mach number aircraft center of mass motion equation set simplified model, the center of mass motion equation set simplified model is:

[0008]

[0009] In the formula, V(t) is the speed of the aircraft, θ(t) and ψ v(t) are the ballistic dip and bank angles, respectively; y(t) is the altitude, x(t) and z(t) are the downrange and crossrange, respectively; m(t) is the vehicle mass, g is the gravitational acceleration, a(t) is the angle of attack, g(t) is the velocity tilt angle v (t) are the ballistic dip and bank angles, respectively; y(t) is the altitude, x(t) and z(t) are the downrange and crossrange, respectively; m(t) is the vehicle mass, g is the gravitational acceleration, a(t) is the angle of attack, g(t) is the velocity tilt angle B is the equilibrium angle of attack, b B is the equilibrium sideslip angle; Y B and Z B are the equilibrium lift and side force corresponding to a B and b B , respectively, and X is the drag force, which is calculated as follows:

[0010]

[0011] where C x and C y are the drag and lift coefficients, q = 0.5 * p * V 2 , where p is the air density at the altitude of the high-Mach vehicle, and S is the reference area of the vehicle;

[0012] S2, determine the constraint conditions and objective function during the quasi-equilibrium flight phase, and describe the trajectory optimization problem as finding the optimal time series of the control variables a and g v that minimize the objective function under the constraint conditions including the control variable constraints, process constraints and terminal constraints of reaching the destination; S3, convert the standard dynamic optimal control problem into a nonlinear programming problem by using the improved direct shooting method;

[0013] S3, convert the standard dynamic optimal control problem into a nonlinear programming problem by using the improved direct shooting method;

[0014] S4, construct the control variable a and g v data set within the control variable constraint range, take the initial state X0 as the input, and incorporate the initial and final time values of the control variables into the optimization design variables, optimize the control variables at the discrete points in the improved direct shooting method by using the adaptive genetic algorithm, and smooth the control variable-time history by cubic spline interpolation, based on the judgment condition of the monotonic change quantity (such as velocity) in the terminal constraint, perform numerical integration based on the fourth-order Runge-Kutta method, and calculate the optimal control u * after several iterations, and simultaneously obtain the optimal trajectory of the state variable curve X * .

[0015] The beneficial effect of the present application is that: the control quantity data set is constructed within the control constraint range, the control quantity initial and final time values and terminal time are included in the optimization design variables, the dynamic optimal control problem with the farthest range as the objective function, meeting the control constraint, process constraint and terminal constraint is parameterized as a nonlinear programming problem by using the improved direct shooting method. On this basis, the control quantity parameters are globally optimized by means of the adaptive genetic algorithm, the control quantity-time history is smoothed by cubic spline interpolation, and the numerical integral is carried out by using the fourth-order Runge-Kutta method, so that the ideal trajectory meeting the conditions is obtained. BRIEF DESCRIPTION OF DRAWINGS

[0016] Figure 1 It is the solving idea of the improved direct shooting method in the trajectory optimization method combining the improved direct shooting method and the adaptive genetic algorithm of the present application;

[0017] Figure 2 It is the flow chart of the adaptive genetic algorithm in the trajectory optimization method combining the improved direct shooting method and the adaptive genetic algorithm of the present application;

[0018] Figure 3 In the trajectory optimization problem with the farthest longitudinal range as the target, the best fitness change curve of simulation is shown;

[0019] Figure 4 In the trajectory optimization problem with the farthest longitudinal range as the target, the speed change curve of simulation is shown;

[0020] Figure 5 In the trajectory optimization problem with the farthest longitudinal range as the target, the trajectory inclination angle change curve of simulation is shown;

[0021] Figure 6 In the trajectory optimization problem with the farthest longitudinal range as the target, the trajectory deflection angle change curve of simulation is shown;

[0022] Figure 7 In the trajectory optimization problem with the farthest longitudinal range as the target, the longitudinal range change curve of simulation is shown;

[0023] Figure 8 In the trajectory optimization problem with the farthest longitudinal range as the target, the height change curve of simulation is shown;

[0024] Figure 9 In the trajectory optimization problem with the farthest longitudinal range as the target, the lateral range change curve of simulation is shown;

[0025] Figure 10 In the trajectory optimization problem with the farthest longitudinal range as the target, the attack angle change curve of simulation is shown;

[0026] Figure 11Speed tilt angle variation curve simulated in the trajectory optimization problem with the target of the maximum longitudinal range

[0027] Figure 12 Altitude-range variation curve simulated in the trajectory optimization problem with the target of the maximum longitudinal range

[0028] Figure 13 Heat flux rate variation curve simulated in the trajectory optimization problem with the target of the maximum longitudinal range

[0029] Figure 14 Dynamic pressure variation curve simulated in the trajectory optimization problem with the target of the maximum longitudinal range

[0030] Figure 15 G-Force variation curve simulated in the trajectory optimization problem with the target of the maximum longitudinal range

[0031] Figure 16 Avoiding no-fly zone plan view simulated in the trajectory optimization problem with the target of the maximum longitudinal range

[0032] Figure 17 Avoiding no-fly zone three-dimensional view simulated in the trajectory optimization problem with the target of the maximum longitudinal range

[0033] Figure 18 Best fitness variation curve simulated in the trajectory optimization problem with the target of the maximum lateral range

[0034] Figure 19 Speed variation curve simulated in the trajectory optimization problem with the target of the maximum lateral range

[0035] Figure 20 Trajectory bank angle variation curve simulated in the trajectory optimization problem with the target of the maximum lateral range

[0036] Figure 21 Trajectory angle of attack variation curve simulated in the trajectory optimization problem with the target of the maximum lateral range

[0037] Figure 22 Longitudinal range variation curve simulated in the trajectory optimization problem with the target of the maximum lateral range

[0038] Figure 23 Altitude variation curve simulated in the trajectory optimization problem with the target of the maximum lateral range

[0039] Figure 24 Lateral range variation curve simulated in the trajectory optimization problem with the target of the maximum lateral range

[0040] Figure 25Trajectory optimization problem with the goal of maximum lateral range, the simulated attack angle change curve;

[0041] Figure 26 Trajectory optimization problem with the goal of maximum lateral range, the simulated speed inclination angle change curve;

[0042] Figure 27 Trajectory optimization problem with the goal of maximum lateral range, the simulated height-longitudinal range change curve;

[0043] Figure 28 Trajectory optimization problem with the goal of maximum lateral range, the simulated heat flow rate change curve;

[0044] Figure 29 Trajectory optimization problem with the goal of maximum lateral range, the simulated dynamic pressure change curve;

[0045] Figure 30 Trajectory optimization problem with the goal of maximum lateral range, the simulated overload change curve;

[0046] Figure 31 Trajectory optimization problem with the goal of maximum lateral range, the simulated evasion of no-fly zone top view;

[0047] Figure 32 Trajectory optimization problem with the goal of maximum lateral range, the simulated evasion of no-fly zone three-dimensional solid view. DETAILED DESCRIPTION

[0048] To make the objects, technical solutions and advantages of the embodiments of the present application clearer, the following will be a clear and complete description of the technical solutions in the embodiments of the present application in conjunction with the accompanying drawings of the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. The components of the embodiments of the present application described and shown in the accompanying drawings can be arranged and designed in various different configurations.

[0049] Therefore, the following detailed description of the embodiments of the present application provided in the accompanying drawings is not intended to limit the scope of the claimed present application, but only represents selected embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without making creative efforts fall within the scope of protection of the present application.

[0050] The specific embodiments of the present application will be described in detail below in conjunction with the accompanying drawings.

[0051] The improved trajectory optimization method combining the direct shooting method and the adaptive genetic algorithm, comprising:

[0052] S1, ignoring the earth rotation and curvature effects, assuming the earth as a plane, the dynamic equations and kinematic equations are established in the trajectory coordinate system and the ground coordinate system respectively, and the 3-DOF high-Mach number aircraft center of mass motion equation set simplified model is obtained, which is:

[0053]

[0054] In the formula: V(t) is the speed of the aircraft, θ(t) and ψ v (t) are the trajectory inclination angle and the trajectory deflection angle respectively; y(t) is the height, x(t) and z(t) are the longitudinal range and the lateral range respectively; m(t) is the mass of the aircraft, 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 sideslip angle; Y B and Z B are the equilibrium lift and the equilibrium side force corresponding to α B and β B respectively, and X is the drag, the calculation expression is as follows:

[0055]

[0056] In the formula, C x and C y are the drag coefficient and the lift coefficient, q=0.5*ρ*V 2 , wherein ρ is the air density at the height where the high-Mach number aircraft is located, and S is the reference area of the aircraft;

[0057] S2, the constraint conditions and the objective function in the flight process of the quasi-equilibrium flight segment are determined, and the trajectory optimization problem is described as finding the optimal time sequence of the control quantities attack angle α and speed inclination angle γ v to make the objective function minimum under the constraint conditions of the center of mass motion equation set simplified model and the constraint conditions, the constraint conditions include the control quantity constraint, the process constraint and the terminal constraint of reaching the destination;

[0058] In order to ensure the stable flight of the aircraft, the control quantities attack angle α and speed inclination angle γ v are limited in a certain range, the control quantity constraint includes attack angle α and speed inclination angle γ v , and the constraint is expressed as:

[0059]

[0060] The process constraint includes the heat flow rate constraint, the dynamic pressure constraint, the overload constraint and the no-fly zone constraint, which are as follows:

[0061] The heat flow rate constraint, the dynamic pressure constraint and the overload constraint are expressed as:

[0062]

[0063] Suppose the type of the forbidden zone that the high-mach flight vehicle avoids in quasi-equilibrium flight is a cylinder, the center point is (x i ,y i ,z i ), and the radius is r i . The specific expression is as follows:

[0064]

[0065] Suppose L(x(t),t) is the shortest distance from the high-mach flight vehicle to the forbidden zone, and ε n is an infinitesimal greater than 0. Then the forbidden zone constraint can be expressed as L(x(t),t)≥ε n , where K Q is a constant parameter related to the high-mach flight vehicle, is the maximum allowable heat flow rate, q dmax is the maximum allowable dynamic pressure, n Lmax is the maximum allowable overload;

[0066] To meet the requirements of mid-end guidance handover, the height, speed, and trajectory inclination angle at the terminal time of quasi-equilibrium flight need to be set as constraints. Therefore, the terminal constraints include the height, speed, and trajectory inclination angle at the terminal time, which are expressed as:

[0067]

[0068] In the formula, y f , V f , θ f are the height, speed, and trajectory inclination angle at the terminal time of quasi-equilibrium flight, respectively; y df , V df , θ df are the predetermined height, speed, and trajectory inclination angle at the terminal time of quasi-equilibrium flight, respectively; ε h , ε v , ε θ are the upper error bounds of the three constraint quantities, respectively;

[0069] The objective function is selected according to the specific task. Generally, it can be divided into maximum range (maximum longitudinal range / maximum lateral range), shortest flight time, maximum speed at the terminal time, etc. In this paper, the maximum range (maximum longitudinal range / maximum lateral range) is selected as the optimization objective, and the objective function is expressed as:

[0070] min J=-L df (7)

[0071] where L df is the range of the flight vehicle;

[0072] The quasi-equilibrium flight trajectory optimization problem of high Mach number aircraft can be described as: finding the control variables of angle of attack a and velocity tilt angle y, which make the performance index minimum under the conditions of satisfying the equations of mass center motion, control variable constraints, process constraints and terminal constraints v The optimal time sequence makes the objective function minimum; the more general optimal control problem is described as follows: finding the control variable u(t), which makes the Bolza type performance index minimum:

[0073]

[0074] and satisfying the constraints of equations of mass center motion

[0075]

[0076] Boundary condition constraints

[0077] φ(x(t0),t0,x(t f ),t f )=0 (10)

[0078] and inequality constraints

[0079] C(x(t),u(t),t)≤0 (11)

[0080] Where t0 is the initial time, t f is the terminal time.

[0081] S3, the standard dynamic optimal control problem is converted into a nonlinear programming problem by using improved direct shooting method;

[0082] Direct shooting method is a method of only discrete control variables, and the state variables describing the motion trajectory need to be obtained by numerical integration of the motion equations according to the parameterized control variables. For the quasi-equilibrium flight trajectory optimization problem of high Mach number aircraft, fourth-order Runge-Kutta method is generally selected for solving to ensure the accuracy and stability of numerical calculation. It should be noted that direct shooting method has high requirements for initial value, and the accuracy of the algorithm often depends on the accuracy of the initial value guess, which may cause local minimum problem. Therefore, the direct shooting is improved, the control variable initial and final time values are taken as optimization design variables, the monotonic change quantity (such as velocity) in the terminal constraint is taken as the judgment condition, and the end time of trajectory is taken as the terminal time, so as to reduce the sensitivity to initial value;

[0083] S31, in the direct shooting method, time variable is generally discretized to obtain discrete time sequence T:

[0084] t0=t1<t2<···<t k-1 <t k= t f <···< t N-1 < t N (12)

[0085] With the terminal constraint of the monotonic change of the variable (such as speed) as the judgment condition, t N A larger value can be taken, such as t N = 1000s, so that the time when the terminal constraint of the monotonic change of the variable ends is t k = t f ;

[0086] S32, discretization of the control variable: taking the control variable at the discrete node as the design variable, a discrete control sequence U is obtained:

[0087] [u1, u2, ···, u k-1 , u k = u f , ···, u N-1 , u N ] (13)

[0088] The control variable between time nodes is obtained by an interpolation function ψ i (t), so the control variable in the continuous space is approximated as follows:

[0089]

[0090] Where u i is the control variable at the grid node;

[0091] S33, discretization of the state variable: the initial state of the high Mach number aircraft is known, that is, X i , under the premise that the discrete control sequence U is obtained by the optimization algorithm, the state variable sequence X corresponding to the time sequence can be obtained by the method based on the fourth-order Runge-Kutta numerical integration through multiple iterative integrations:

[0092] [X1, X2, ···, X k-1 , X k = X f , ···, X N-1 , X N ] (15)

[0093] Where X k = X f is the terminal state corresponding to the terminal time t f ;

[0094] S34, set the variable set as:

[0095] Y = [X1, u1, t1, X2, u2, t1, ···, X k-1 , u k-1 , t k-1 , X f , u f , t f ] (16)

[0096] The performance index can be expressed as:

[0097]

[0098] The boundary conditions can be expressed as:

[0099] g1(Y) = 0 (18)

[0100] The process constraints can be expressed as:

[0101] g2(Y)≤0 (19)

[0102] The nonlinear programming problem can be expressed as:

[0103]

[0104] S4, constructing control attack angle a and velocity tilt angle g in the control constraint range v data set, taking the initial state X0 as input, taking the control initial and final time values into the optimization design variables, using the adaptive genetic algorithm to optimize and improve the control at the discrete points in the direct shooting method, as shown in Figure 2 , and smoothing the control-time history by cubic spline interpolation, taking the monotonic change amount (such as velocity) in the terminal constraint as the judgment condition, and performing numerical integration based on the fourth-order Runge-Kutta method, and after several iterations, the optimal control u * is calculated, and the optimal trajectory of the state variable curve X * is obtained at the same time.

[0105] S41, population initialization coding: constructing control data set in the control constraint range, selecting several time points according to formula (12) based on the improved direct shooting method, and the two one-dimensional arrays a = (a1, a2, ···, a N-1 , a N ) and g = (g1, g2, ···, g composed of the control attack angle and the velocity tilt angle at each node time are taken as the chromosome of the trajectory individual;

[0106] There are several advantages to choosing this coding:

[0107] 1) The control constraint formula (3) can be automatically satisfied;

[0108] 2) In three degrees of freedom trajectory, when the control variables are the angle of attack and the velocity tilt angle, the complete trajectory can be obtained directly. If other parameters such as the trajectory inclination angle and the height are selected as the control variables, the angle of attack and the velocity tilt angle need to be calculated in the iteration, which increases the time and complexity of the calculation.

[0109] S42, the number of angle of attack in the given aerodynamic data is limited, and more data needs to be obtained by interpolation. In addition, the angle of attack value at the node time is not enough, and interpolation calculation is needed. The main interpolation methods include linear interpolation, polynomial interpolation, spline interpolation, etc. This paper uses cubic spline interpolation, which has the following advantages:

[0110] 1) Control constraint (3) can be automatically satisfied;

[0111] 2) The algorithm is easy to implement, and the interpolation error is smaller than linear interpolation. Using spline is easier to evaluate than using high-order polynomial, and it is not affected by the Runge phenomenon;

[0112] 3) The angle of attack obtained by spline interpolation is a smooth curve, not a broken line of linear interpolation, which is beneficial to the design of the aircraft control system.

[0113] S43, the fitness function is an index for evaluating the chromosome, and it is also the key to the "survival of the fittest" of the self-adaptive genetic algorithm. This paper uses Joines method to construct a comprehensive fitness function: based on the multi-objective hierarchical planning method, the performance index at the terminal time and the terminal constraint are weighted to obtain the fitness function. The control constraint has been automatically satisfied, and other process constraints (dynamic pressure, overload, heat flow rate, and no-fly zone) can be set to infinity if they do not satisfy the maximum dynamic pressure q dmax , the maximum overload n Lmax , the maximum heat flow rate and the no-fly zone. Joines method is improved based on the static penalty function, and the penalty coefficient is self-adaptive with the increase of the iteration number. The specific form is:

[0114]

[0115] Where, R1=(cd) λ , f i (x) is the penalty function item of violating the i-th constraint, p i is the weight coefficient of the i-th constraint penalty item, c, λ, are constants, d is the iteration number; J is the trajectory optimization objective function item, q is the corresponding weight coefficient, K is a constant, and numerical simulation verification shows that c=0.5, is a more feasible solution. In addition, the difference between the terminal constraint terms in the comprehensive fitness function is small after logarithmic processing, so the weight coefficient q of the next level performance index can be mainly adjusted;

[0116] S44, if yes, the monotonic change amount (such as speed) in the terminal constraint is taken as the judgment condition, numerical integration is carried out based on the fourth-order Runge-Kutta method, and the optimal control u is calculated through several iterations * , and the optimal trajectory X of the state variable curve can be obtained * , as shown in the figure. Otherwise, go to S45; Figure 1

[0117] S45, the individuals are selected, crossed and mutated by genetic operation according to the fitness value, a new generation population is generated, and S43 is entered.

[0118] Selection genetic operation: the operation of selecting superior individuals from the population and eliminating inferior individuals is called selection. The purpose of selection operation is to directly inherit the optimized individuals to the next generation or to generate new individuals through pairing and crossing and then to inherit them to the next generation. Usually, roulette selection method is used. Let the population size be Q, and the fitness of individual i be F i , then the probability of selecting the ith individual is:

[0119]

[0120] Crossing genetic operation: two chromosomes are crossed and combined according to the crossing probability, so as to expect that the excellent characteristics of the parent individuals are inherited to the offspring and new excellent individuals are generated. In this paper, real number crossing method is used, that is, the crossing operation of the kth chromosome a k and the lth chromosome at j position is:

[0121]

[0122] Mutation genetic operation: the mutation of a gene value of an individual string in the population is carried out. The mutation operation method of the jth gene a ij of the ith individual is:

[0123]

[0124] , 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, and G​max Max is the maximum number of evolution, r is a random number in [0, 1];

[0125] The mutation operation has two purposes, one is to enhance the local optimization ability of adaptive genetic algorithm, and the other is to maintain population diversity to prevent premature convergence.

[0126] Simulation

[0127] The basic parameters of the aircraft and the earth are shown in Table 1:

[0128] Table 1

[0129]

[0130] The control constraints are shown in Table 2:

[0131] Table 2

[0132]

[0133] The process constraints are shown in Table 3:

[0134] Table 3

[0135]

[0136] The boundary constraints (initial state and terminal state) are shown in Table 4:

[0137] Table 4

[0138]

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

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

[0141] Table 5

[0142]

[0143]

[0144] The adaptive genetic algorithm designed in the method is mainly improved for the fitness function in the genetic algorithm, and other parameters can be set in the following ranges: population size: 20-100; iteration number: 100-500; crossover probability: 0.4-0.99; mutation probability: 0.0001-0.1. In order to prevent premature, generally increase the mutation rate or increase the population size to solve, so the maximum value in the range of mutation rate, the population size can be increased appropriately, the minimum value in the range of iteration number is selected to reduce the amount of calculation, the crossover operation is the core of the global search ability of the adaptive genetic algorithm, and a larger value is selected, and the specific values are shown in Table 5.

[0145] In order to verify the robustness of the improved direct targeting method-adaptive genetic algorithm constructed, the aerodynamic parameters in the actual flight process are considered as uncertainties, which is more reasonable. Without loss of generality, on the basis of the original aerodynamic parameters, a normal distribution white noise is introduced, and the specific formula is as follows:

[0146]

[0147] Wherein, f(x i ) represents the original aerodynamic coefficient (including the lift coefficient C x and the resistance coefficient C y Lift coefficient), f δ (x i ) is the aerodynamic coefficient added with noise, δ is the relative error level, is a normal distribution random number with mean 0 and variance 1.

[0148] Optimization of the farthest longitudinal trajectory

[0149] 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:

[0150] J0=min(-L0)=min(-x(t f )) (26)

[0151] The original algorithm (i.e. direct targeting method-genetic algorithm, which needs to test the initial value constantly, and the fitness function needs to be logarithmically transformed, otherwise it is difficult to converge), the hybrid algorithm constructed: improved direct targeting method-adaptive genetic algorithm and the white noise introduced on this basis are simulated, and the results are shown in Figures 3-17

[0152] Optimization of the farthest lateral trajectory

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

[0154] J0=min(-L0)=min(-z(t f ​(27)

[0155] Using the original algorithm (i.e., direct shooting method-genetic algorithm, which requires continuous testing of initial values, and the fitness function needs logarithmic transformation, otherwise convergence is difficult), a hybrid algorithm was constructed: an improved direct shooting method-adaptive genetic algorithm, and white noise was introduced into the simulation. The results are as follows. Figures 18-32 As shown.

[0156] Figure 3 and Figure 18 The graph shows the evolution of the optimal fitness value in the population over iterations. As can be seen, the fitness function of the original genetic algorithm converges around generation 40. The optimal fitness value of the constructed adaptive algorithm exhibits significant initial fluctuations but gradually converges with increasing iterations, stabilizing around generation 10. In the later stages, from generation 60 to 90, the optimal fitness function value of the hybrid algorithm for optimizing the control variable shows slight variations, but these are relatively small, indicating a faster convergence speed compared to the original genetic algorithm. This demonstrates that the adaptive genetic algorithm has a better convergence performance for optimizing control variables.

[0157] It should be noted that the original algorithm (direct shooting method - genetic algorithm) obtained its trajectory through extensive trial and error with initial values. In most cases, due to its sensitivity to initial values, it often fails to produce a trajectory that avoids no-fly zones. The constructed hybrid algorithm, however, incorporates the initial and final values ​​of the control variables into the optimization design variables, using monotonically changing quantities in the terminal constraints (such as velocity) as the criteria for determining the trajectory's end time. This eliminates the need for continuous trial and error with initial values, resulting in higher efficiency. Figures 10-11 As can be seen in 25-26, the constructed hybrid algorithm, the improved direct shooting method-adaptive genetic algorithm, controls the constraints within the specified range in simulations where the performance index is the longest longitudinal / longest transverse stroke. Figures 13-15 In sections 28-30, the general constraint terms such as heat flux, overload, and dynamic pressure in the process constraints are also within the constraint range; Figures 16-17 As can be clearly seen in 31-32, the hybrid algorithm can avoid multiple no-fly zones, thus meeting practical needs.

[0158] Finally, considering the model's uncertainties, namely the introduction of normally distributed white noise into the aerodynamic coefficients, the optimal fitness function curve tends to converge around generation 20-30. The adaptive genetic algorithm still has a good convergence effect for optimizing control variables, and the corresponding control constraints and process constraints (heat flow rate, overload, dynamic pressure, and multiple no-fly zones) can still meet the requirements, verifying that the constructed hybrid algorithm has a certain degree of robustness.

[0159] To more intuitively reflect that the constructed hybrid algorithm satisfies the terminal constraints, the following table 6-7 is established:

[0160] Table 6 Terminal constraint corresponding terminal state (performance index: maximum longitudinal range)

[0161]

[0162] Table 7 Terminal constraint corresponding terminal state (performance index: maximum lateral range)

[0163]

[0164] From Figures 4-5 , 7-8 and Figures 19-20 , 23-24 and Tables 6-7, it can be seen that in the simulation with the maximum longitudinal range / maximum lateral range as the performance index, the hybrid algorithm can meet the terminal constraint, and the performance index maximum longitudinal range / maximum lateral range value is larger than that of the original algorithm, closer to the global optimal solution, verifying the effectiveness of introducing the multi-objective hierarchical planning method to ensure that the terminal constraint is met at the terminal time while the performance index is optimal. In addition, in the simulation considering the disturbance of aerodynamic parameters, the terminal constraint can still meet the requirements.

[0165] In summary, the optimal performance index trajectory (maximum longitudinal range / maximum lateral range) obtained based on the improved direct shooting method and adaptive genetic algorithm meets the control constraints, process constraints (dynamic pressure, overload, heat flow rate and flight-prohibited area) and terminal constraint conditions, and can avoid multiple flight-prohibited areas, which illustrates the effectiveness of the constructed hybrid algorithm. Secondly, compared with the original algorithm (direct shooting method-genetic algorithm), the hybrid algorithm does not need to constantly test the initial value, has faster convergence speed, better performance index and can quickly and accurately approach the global optimal solution. Finally, when considering the model uncertainty, i.e. the disturbance of aerodynamic parameters by normal distribution white noise, the requirements can still be met, which illustrates that the constructed hybrid algorithm has certain robustness.

[0166] Conclusion

[0167] The trajectory optimization method for quasi-equilibrium flight of high Mach number vehicles is a key technology in the research of high Mach number vehicles. Considering the complex constraint limitations including multiple flight-prohibited areas during flight, a hybrid optimization solution based on improved direct shooting method and adaptive genetic algorithm is proposed.

[0168] 1) By improving the direct shooting method, the initial and final time values of the control variables are included in the optimization design variables, the monotonic change quantity (such as velocity) in the terminal constraint is used as the judgment condition, and the end time of the trajectory is taken as the terminal time, which reduces the sensitivity to the initial value.

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

[0170] 3) Based on the established adaptive genetic algorithm optimization to improve the control amount at discrete points in the direct targeting method, the control amount-time history is smoothed by cubic spline interpolation, and the fourth order Runge-Kutta method is used for numerical integration, and after several iterations, the optimal control time sequence is calculated, and the ideal trajectory that can avoid multiple forbidden flight zones, meet the constraint conditions and ensure the longest distance is obtained.

[0171] Experiments show that the established hybrid optimization algorithm does not need to constantly test the initial value, compared with the original algorithm (direct targeting method-genetic algorithm), the convergence speed is faster, the performance index farthest range / longest lateral range value is larger, and the global optimal solution is closest. In addition, considering the uncertainty factors existing in the model in the actual flight process, the normal distribution white noise is introduced into the aerodynamic parameters, and the results show that various constraint conditions can still meet the requirements in the quasi-equilibrium flight process, which shows that the constructed hybrid algorithm has certain robustness. In terms of application, the constructed hybrid optimization algorithm is essentially an optimization method, which is not only suitable for quasi-equilibrium flight trajectory optimization of high Mach number vehicles, but also suitable for general vehicle trajectory optimization problems, and provides certain reference value for the next step of studying the robust trajectory optimization problem of uncertainty quantization.

[0172] The technical scheme of the present application is not limited to the restriction of the above specific embodiments, and any technical deformation made according to the technical scheme of the present application falls within the protection scope of the present application.

Claims

1. An improved ballistic optimization method combining direct target acquisition and adaptive genetic algorithm, 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, respectively, resulting in a simplified model of the center-of-mass motion equations for a 3-DOF high Mach number spacecraft. The simplified model of the center-of-mass motion equations is as follows: ， In the formula: V(t) is the velocity of the aircraft, θ(t) and y(t) represents the trajectory inclination angle and trajectory deflection angle, respectively; y(t) represents altitude; x(t) and z(t) represent the longitudinal and transverse strokes, respectively; m represents the mass of the aircraft; and g represents the acceleration due to gravity. For the angle of attack, The velocity tilt angle; To balance the angle of attack, To balance the sideslip angle; and They are respectively , The corresponding equilibrium lift and equilibrium lateral force, where X is the drag, are calculated using the following expressions: ， In the formula , 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 during the quasi-equilibrium flight phase, and describe the ballistic optimization problem as finding the control angle of attack under the simplified model and constraints of the center of mass motion equations. and velocity tilt angle The standard dynamic optimal control problem is to minimize the objective function by optimizing the time series. The constraints include control quantity constraints, process constraints, and terminal constraints to reach the destination. S3. The standard dynamic optimal control problem is transformed into a nonlinear programming problem by using the improved direct shooting method. S4. Construct the angle of attack of the control quantity within the control quantity constraint range. and velocity tilt angle The dataset will contain the initial state. As input, the initial and final values ​​of the control input are included in the optimization design variables. An adaptive genetic algorithm is used to optimize and improve the control input at discrete points in the direct shooting method. The control input-time history is smoothed by cubic spline interpolation. The monotonic change in the terminal constraint is used as the judgment condition, with the monotonic change being velocity. Numerical integration is performed based on the fourth-order Runge-Kutta method, and the optimal control is calculated after several iterations. At the same time, the optimal trajectory of the state variable curve can be obtained. .

2. The ballistic optimization method combining the improved direct firing method and adaptive genetic algorithm according to claim 1, characterized in that, In S2, the control constraints include the angle of attack. and velocity tilt angle The constraint is represented as: ， Process constraints include heat flow rate constraints, dynamic pressure constraints, overload constraints, and no-fly zone constraints, as detailed below: Heat flow rate constraint, dynamic pressure constraint, and overload constraint are expressed as follows: ， Assume that the no-fly zone avoided by the high Mach number aircraft during quasi-equilibrium flight is cylindrical, with its center point being... , radius is The specific expression is as follows: ， set up The shortest distance from a high Mach number aircraft to a no-fly zone. If is an infinitesimal quantity greater than 0, 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: ， In the formula, , , These are the altitude, velocity, and trajectory inclination at the quasi-equilibrium flight terminal moment, respectively. , , These are the altitude, speed, and trajectory inclination at the predetermined quasi-equilibrium flight terminal moment, respectively. , , These are the upper bounds of the errors for the three constraint quantities; The objective function is expressed as: ， in For the range of the aircraft; The trajectory optimization problem for the quasi-equilibrium flight phase of a high Mach number aircraft can be described as: finding the angle of attack of the control variable under the conditions of satisfying the equations of motion of the center of mass, control constraints, process constraints, and terminal constraints. and velocity tilt angle The optimal time series minimizes the objective function: ， And it satisfies the constraints of the equations of motion of the center of mass. ， Boundary condition constraints ， and inequality constraints ， in, At the starting time, This refers to the terminal moment.

3. The ballistic optimization method combining the improved direct firing method and adaptive genetic algorithm according to claim 1, characterized in that, In S3, S31. Discretize the time variable to obtain the discrete time series T: ， Using the monotonic change in the terminal constraint as the judgment condition, and the monotonic change as velocity, we can set the end time of the monotonic change in the terminal constraint to be... ; S32. Discretize the control variables: Using the control variables at discrete nodes as design variables, obtain the discrete control sequence. : ， The control quantities between time points are obtained through interpolation functions. If obtained, the control variables in the continuous space can be approximated as follows: ， in These are the control variables at the grid nodes; S33. Discretize the state variables: The initial state of the high Mach number aircraft is known, i.e. Discrete control sequences are obtained through optimization algorithms. Under the premise that this is the case, the state variable sequence corresponding to the time series can be obtained through multiple iterations of numerical integration based on the fourth-order Runge-Kutta method. : ， In equation (15) For terminal time Corresponding terminal status; S34. Let the set of variables be: ， Performance metrics can be expressed as: ， The boundary conditions can be expressed as: ， Process constraints can be expressed as: ， The nonlinear programming problem is represented as: 。 4. The ballistic optimization method combining the improved direct firing method and adaptive genetic algorithm according to claim 1, characterized in that, Specifically, S4 includes: S41. Population Initialization Encoding: Construct a control quantity dataset within the control quantity constraints to ensure that the control constraints are met. Then, based on the discrete time series method in the improved direct shooting method, select several time points. For each node, the control quantity angle of attack and velocity tilt angle form two one-dimensional arrays. , Chromosomes as a ballistic individual; S42. Given aerodynamic data, the number of angles of attack is finite, and the angles of attack values ​​at noon are insufficient. Cubic spline interpolation is required to meet the requirements. S43. Calculate the fitness value of each individual in the population: Based on the Joints method and multi-objective hierarchical programming, the performance index at the terminal time and the terminal constraints are weighted to obtain the comprehensive fitness function. The control constraints have been satisfied, and the process constraints can be set to allow for the maximum dynamic pressure not to be satisfied. Maximum overload Maximum heat flux The fitness function value for crossing no-fly zones is set to infinity, and the fitness function is expressed as: ， in, , For the penalty function term that violates the i-th terminal constraint, Let c, λ, and c be the weight coefficients of the i-th terminal constraint penalty term. All are constants, d is the number of iterations; J is the objective function term for ballistic optimization, q is the corresponding weight coefficient, and K is a constant; S44. Determine if the maximum number of iterations has been reached. If so, use the monotonic change in the terminal constraint as the criterion, perform numerical integration based on the fourth-order Runge-Kutta method, and calculate the optimal control after several iterations. At the same time, the optimal trajectory of the state variable curve can be obtained. Conversely, it enters S45; S45. 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 S43.

5. The ballistic optimization method combining the improved direct firing method and adaptive genetic algorithm according to claim 4, characterized in that, In S45: Selection genetic operation: Superior individuals are selected from the population with a certain probability to form a new population for reproduction, resulting in the next generation of individuals. Let the population size be Q, and the fitness of individual i be . Then the probability of the i-th individual being selected is: ， Crossover inheritance: This operation combines two chromosomes based on crossover probability, aiming to pass on superior traits from the parent to the offspring, thus creating new superior individuals. In the real number crossover method, the k-th chromosome... and the The crossover operation at position j of each chromosome is as follows: ， Genetic manipulation of variation: the j-th gene of the i-th individual The method for performing mutation operations is as follows: ， in, Let be the fitness value of individual i, Q be the number of individuals in the population, and b be a random number in the range [0,1]. , Genes The upper and lower bounds; , It is a random number. It is the current iteration number. The maximum number of evolutions is r, which is a random number in the range [0,1].

Images