Re-entry trajectory optimization method for reusable vehicle based on particle swarm optimization

Abstract

The present invention belongs to the technical field of aircraft trajectory optimization, and relates to a reusable aircraft reentry trajectory optimization method based on a particle swarm algorithm. The present invention establishes a three-degree-of-freedom motion model and an aerodynamic model of an aircraft taking into account the rotation of the earth, and analyzes the aerodynamic characteristics of the aircraft. On this basis, by planning the angle of attack profile, the process constraints are mathematically transformed to form a reentry corridor. The optimization parameters of the reference profile are designed and simplified, the basic principles for optimizing the remaining parameters are proposed, and a multi-objective optimization function for the optimization parameters is designed to convert the trajectory optimization problem into a parameter optimization problem. Finally, the reentry trajectory is optimized by combining the particle swarm algorithm and the trajectory parameter solution algorithm. In view of the shortcomings of the traditional trajectory optimization method, the trajectory optimization method of the present invention randomly generates initial values, solves the problem of initial value sensitivity, avoids the optimization result from falling into the local optimum, and has broad application prospects.

Description

Technical Field

[0001] The invention belongs to the technical field of aircraft trajectory optimization, and relates to a reusable aircraft reentry trajectory optimization method based on a particle swarm algorithm. Background Art

[0002] Reusable vehicle reentry means consuming the huge kinetic energy of the vehicle while satisfying constraints such as heat flow, overload and dynamic pressure, and ensuring its smooth transition to the energy management window in the terminal area by controlling the flight trajectory. This section has the characteristics of high flight altitude, high speed, long range, large flight envelope and strong uncertainty in aerodynamic parameters. In order to successfully complete the reentry flight mission, a reasonable reentry trajectory must be designed.

[0003] According to the difference in design ideas, the design of reentry trajectory can be divided into two categories: trajectory design technology based on flight corridor and trajectory optimization technology. Faced with the huge computational requirements of the trajectory design method based on flight corridor and the limitation that it is difficult to optimize performance through iteration, trajectory optimization technology was born.

[0004] The trajectory optimization problem refers to finding a trajectory that satisfies all flight restrictions and maximizes performance indicators under given flight missions and specific environmental constraints. It is essentially a nonlinear, multi-constrained optimal control problem. Generally, traditional optimization methods are divided into indirect methods based on the maximum principle and direct methods based on nonlinear programming theory, including direct shooting method, differential inclusion method, dynamic inverse method, and point matching method. The indirect method converts the problem into a two-point boundary value problem by deriving the first-order necessary conditions of the optimal control problem. The direct method uses a discretization method to discretize the aircraft equations along the trajectory direction, thereby converting the continuous-time optimal control problem into a nonlinear programming problem with finite discrete points. Since both the direct method and the indirect method use traditional optimization methods, their disadvantages are that they are very sensitive to initial values ​​and have poor real-time performance, and are prone to local convergence. New optimization algorithms are urgently needed to solve the problems of traditional algorithms. Summary of the invention

[0005] Traditional trajectory optimization methods are prone to fall into local optimality when performing trajectory optimization, and it is difficult to obtain the optimal solution. The present invention designs a reusable aircraft reentry trajectory optimization method based on a particle swarm algorithm, which can avoid taking initial values. It mainly uses a particle swarm optimization algorithm to randomly generate a large number of initial values ​​within the activity range of the particle swarm, converts the trajectory optimization problem into an endpoint parameter optimization problem, proposes the basic principles of parameter optimization, simplifies profile parameters, introduces individual optimality and group optimality, and performs iterative solutions. Compared with traditional trajectory optimization methods, this method solves the problem of initial value sensitivity and prevents falling into local optimality during the solution process.

[0006] The technical solution of the present invention:

[0007] A reusable vehicle reentry trajectory optimization method based on particle swarm algorithm, the steps are as follows:

[0008] Step (1) Mathematical modeling and characteristic analysis of reusable aircraft

[0009] (1.1) Establishment of three-degree-of-freedom dynamic model

[0010] The three-degree-of-freedom particle dynamics equation considering the earth's rotation is established in the velocity coordinate system:

[0011]

[0012] Where: is the distance from the center of the earth; is latitude; is the longitude; is the speed of the aircraft relative to the earth; is the mass of the aircraft; is the track inclination; is the track azimuth (north-west is positive); is the gravitational acceleration, where is the standard gravitational acceleration on the surface; is the tilt angle; is the angular velocity of the Earth's rotation; , are the drag and lift of the aircraft respectively.

[0013] To simplify the trajectory calculation, the centrifugal force caused by the rotation of the earth is ignored, and the two-dimensional particle dynamics equation is obtained:

[0014]

[0015] (1.2) Establishing the aerodynamic model of the aircraft

[0016] Aerodynamic force can be decomposed into lift in the velocity coordinate system and resistance , and its calculation formula is:

[0017]

[0018]

[0019] Where: is the atmospheric density, where is the standard atmospheric density, is the current height of the aircraft from the ground, is the atmospheric normalized height; is the wing reference area; and are the lift coefficient and the drag coefficient respectively.

[0020] By fitting the lift-drag characteristics of the aircraft, the and The aerodynamic coefficient fitting model.

[0021]

[0022] Where: , , , , , , , All are aerodynamic fitting coefficients; is the angle of attack of the aircraft.

[0023] (1.3) Analysis of aircraft aerodynamic characteristics

[0024] By observing how lift, drag and lift-to-drag ratio change with angle of attack at the same altitude and different speeds, the aerodynamic characteristics of the aircraft are analyzed and the angle of attack corresponding to the maximum lift-to-drag ratio is determined.

[0025] Step (2) Reentry Corridor

[0026] (2.1) Planning the angle of attack profile

[0027] Referring to the typical angle of attack profile design, by analyzing the impact of the angle of attack on thermal flow constraints, overload constraints, dynamic pressure constraints and other constraints, the size of the angle of attack and the timing of downward pressure are selected, thereby determining the size and shape of the angle of attack profile.

[0028] (2.2) Calculation of reentry corridor

[0029] The reentry process constraints are mainly divided into heat flow constraints, overload constraints, dynamic pressure constraints, and balanced gliding constraints when the roll angle is zero. The expression is:

[0030]

[0031]

[0032]

[0033]

[0034] Where: and are heat flow rate and maximum heat flow rate respectively; is the heat transfer coefficient; and They are pneumatic overload and maximum pneumatic overload respectively; and are the dynamic pressure and the maximum dynamic pressure respectively.

[0035] Through mathematical transformation of process constraints and atmospheric density expressions, multiple process constraints are converted into constraints of a single physical quantity. The three process constraints of heat flow, overload and dynamic pressure constitute the hard constraints of the reentry corridor, and the balanced gliding constraint constitutes the soft constraints of the reentry corridor, which together form the upper and lower boundaries of the reentry corridor, and the altitude-speed corridor ( ) and the drag acceleration-velocity corridor ( ). The expressions for the altitude-speed and drag acceleration-speed reentry corridors are:

[0036]

[0037]

[0038]

[0039]

[0040] Where: Heat flow constraint The corresponding height boundary; Overload constraint The corresponding height boundary; Dynamic pressure constraint The corresponding height boundary; The height boundary corresponding to the balanced gliding constraint; Heat flow constraint The corresponding drag acceleration boundary; Overload constraint The corresponding drag acceleration boundary; Dynamic pressure constraint The corresponding drag acceleration boundary; is the drag acceleration boundary corresponding to the balanced gliding constraint.

[0041] Step (3) Propose an optimization problem

[0042] The reference resistance acceleration-velocity profile is generally composed of four sections: temperature control section, balanced gliding section, constant resistance section and transition section. It is generally designed as a concatenation of four linear functions. The endpoints of each section are considered as parameters to be optimized. The optimization parameters of the reference profile are simplified, and five parameters to be optimized are obtained:

[0043] X = [ a 1 , V 2 , V 3 , Da 3 , V 4 ]

[0044] Where: for location; is the linear function slope of the temperature control section; is the terminal speed of the temperature control section; To balance the terminal speed of the glide segment; To balance the terminal drag acceleration of the gliding segment; is the end resistance acceleration of the constant resistance section.

[0045] From the formula of the reentry corridor in step (2), we can know that the upper and lower boundaries of the reentry corridor are The functions are

[0046]

[0047]

[0048] The reference drag acceleration profile parameters are given by The only certainty can be described as and The function of

[0049]

[0050] According to the principle of continuity of flight mission and reference profile, the basic principle of optimization of remaining parameters is proposed. The final optimization problem is: given the initial conditions, under the conditions of satisfying flight constraints, the trajectory parameters are optimized to reduce the total heating of the process and increase the range, and to establish as large an overload as possible to ensure the maneuverability of the aircraft, so as to plan an optimal trajectory. The expression of the optimization objective function is:

[0051]

[0052] Where: is the target optimization function, i.e., the fitness of the particle swarm algorithm; and The initial and terminal times of the spacecraft reentry; and are the upper and lower boundaries of the drag acceleration reentry corridor; is the total flight distance; , , They are the weight coefficients of total heating, total range and overload respectively; is the absolute value of the slope of the temperature control section; and are the absolute values ​​of the minimum and maximum slopes that the aircraft can achieve in the temperature control section without exceeding the control quantity constraint; and The maximum heat flux decreases and The corresponding speed; is the speed corresponding to the downward pressure of the angle of attack; is the end resistance acceleration of the temperature control section; for The maximum value possible; for The corresponding upper reentry corridor corresponds to a smaller speed; is the threshold value.

[0053] Step (4) Trajectory optimization based on particle swarm algorithm

[0054] (4.1) Planning the activity range of particle swarm

[0055] According to the constraints of the optimization parameters in step (3), the activity range of the particle swarm is planned.

[0056]

[0057] (4.2) Initialize particle swarm

[0058] Initialize the particle swarm parameters, i.e. the number of particles ; Particle Dimension ; Number of iterations ; Inertia weight ; Individual learning factor ; Group learning factor ; Position maximum value and maximum speed ; and randomly generate The position of the particle (i.e. the parameter value to be optimized) and speed , forming a particle swarm , whose expression is:

[0059]

[0060]

[0061] (4.3) Calculate the fitness of the initial particles

[0062] The position values ​​of each particle in the initial particle group will constitute a corresponding independent drag acceleration profile. Through the three-degree-of-freedom dynamic model of step (1), the target optimization function value of each particle in the initialization group can be obtained by numerical integration, that is, the fitness .

[0063] (4.4) Individual optimization and global optimization

[0064] By comparing the fitness of each particle , we can get the individual minimum fitness of the initial particle swarm and optimal location Minimum fitness of the group and optimal location ;

[0065] (4.5) Update particle swarm and cross-bounds check

[0066] Combining individual optimization and group optimization, expressions for updating speed and position are given. In order to prevent the updated particle speed and position from exceeding the constraint range, an out-of-bounds check is required.

[0067] (4.6) Calculate the updated particle fitness and optimize

[0068] The updated position values ​​of the particle swarm will form a corresponding independent drag acceleration profile, which will be passed through the three-degree-of-freedom dynamic model of step (1) again to obtain the first The fitness of each particle after the update And update the individual minimum fitness and optimal location And the minimum fitness of the group and optimal location .

[0069] (4.7) Calculation termination judgment

[0070] If the optimization target value does not meet the requirements and the number of iterations If the maximum value is not reached, return to (4.5); if the requirements are met, end the loop.

[0071] Step (5) Trajectory parameter calculation

[0072] Through step (4), the optimized drag acceleration profile can be obtained. Combined with the three-degree-of-freedom dynamic model in step (1), the drag acceleration-velocity profile is converted into the altitude-velocity profile and the range and its rate of change, sinking rate, trajectory angle and its rate of change, roll angle, speed change rate, etc. at the corresponding speed through the trajectory parameter solution algorithm, and the aircraft state quantity at the corresponding time can be obtained. x = [ H , V , γ ] T and control volume , as shown below:

[0073]

[0074] The beneficial effects of the present invention are:

[0075] The present invention firstly establishes a three-degree-of-freedom motion model and an aerodynamic model of the aircraft considering the rotation of the earth. Since the aerodynamic parameters and lift-drag characteristics determine the trajectory tracking capability and maneuverability of the aircraft, the aerodynamic characteristics of the aircraft are analyzed. On this basis, by planning the angle of attack profile, the process constraints such as heat flow, overload, dynamic pressure and balanced gliding are mathematically transformed into the constraints of a single physical quantity, thereby forming a reentry corridor. According to the continuity principle of the flight mission and the reference profile, the optimization parameters of the reference profile are simplified, and the basic principle of the optimization of the remaining parameters is proposed. The total heat flow, range and overload of the reentry trajectory are taken as the comprehensive optimization targets, and a multi-objective optimization function for the optimization parameters is designed to convert the trajectory optimization problem into a parameter optimization problem. Finally, the reentry trajectory is optimized by combining the particle swarm algorithm and the trajectory parameter solution algorithm. In view of the shortcomings of the traditional trajectory optimization method, the trajectory optimization method of the present invention solves the problem of initial value sensitivity, avoids the optimization result from falling into the local optimum, and has broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 It is a flow chart of a reusable reentry trajectory optimization method based on particle swarm algorithm;

[0077] Figure 2 It is a graph of lift changing with angle of attack;

[0078] Figure 3 It is a graph of drag force versus angle of attack;

[0079] Figure 4 It is the curve of lift-to-drag ratio changing with angle of attack;

[0080] Figure 5 It is the angle of attack profile design diagram;

[0081] Figure 6 It is the altitude-velocity reentry corridor diagram;

[0082] Figure 7 is the drag acceleration-velocity reentry corridor diagram;

[0083] Figure 8 is the reference drag acceleration-velocity profile;

[0084] Fig. 9 This is a flow chart of trajectory optimization based on particle swarm algorithm;

[0085] Fig.10 It is the flow chart of trajectory parameter solution;

[0086] Fig.11 It is the activity range diagram of particle swarm;

[0087] Fig.12 It is the particle swarm iteration convergence diagram;

[0088] Fig.13 It is a graph showing the convergence of the group's optimal fitness with the number of iterations;

[0089] Fig.14 It is the optimization profile result diagram. DETAILED DESCRIPTION

[0090] The specific implementation of the present invention is further described below in conjunction with the accompanying drawings and technical solutions.

[0091] Firstly, the mathematical modeling of the target aircraft was carried out, mainly including the three-degree-of-freedom dynamic model and aerodynamic model considering the rotation of the earth, and the aerodynamic characteristics were analyzed on this basis. The angle of attack profile was planned, and the reentry constraint was converted into a reentry corridor. The basic parameters of the reference profile were designed, the basic principles of parameter optimization were proposed, the profile parameters were simplified, and the optimization problem was proposed. The optimal trajectory was obtained by iterative solution through the particle swarm algorithm. The state quantity of the aircraft during the reentry process was calculated through the trajectory parameter solution algorithm and the dynamic model. A comparative analysis was conducted with the traditional sequential quadratic programming method.

[0092] The process of a reusable reentry trajectory optimization method based on particle swarm algorithm is as follows: Figure 1 As shown, this implementation case will be described using the CAV-H general purpose aircraft, as follows:

[0093] (1) Mathematical modeling and characteristic analysis of reusable aircraft

[0094] (1.1) Establishment of three-degree-of-freedom dynamic model

[0095] Ignoring the earth's oblateness and considering the earth's rotation, the three-degree-of-freedom dynamic equation of the re-entry of a reusable spacecraft is established with time as the independent variable.

[0096] (1)

[0097] To simplify the trajectory algorithm, the centrifugal force caused by the rotation of the earth can be ignored. When , according to formula (1), the two-dimensional particle dynamics equation of the aircraft in the velocity coordinate system is derived.

[0098] (2)

[0099] (3)

[0100] (1.2) Establishing the aerodynamic model of the aircraft

[0101] Aerodynamic force is decomposed into lift in the velocity coordinate system and resistance , and its calculation formula is

[0102] (4)

[0103] (5)

[0104] Sample aircraft and It is only related to the angle of attack. The aerodynamic coefficient fitting model of this example is shown in formula (6).

[0105] (6)

[0106] Where: , , , , , , , .

[0107] (1.3) Analysis of aircraft aerodynamic characteristics

[0108] The aerodynamic characteristics of the aircraft are analyzed to give the changes of lift, drag and lift-to-drag ratio with angle of attack at the same altitude and different speeds, and to determine the angle of attack corresponding to the maximum lift-to-drag ratio. Figure 2 Figure 4 It can be seen from the figure that lift and drag increase with the increase of angle of attack; Figure 4 The angle of attack corresponding to the aircraft's maximum lift-to-drag ratio is approximately 15°, and basically does not change with speed.

[0109] (2) Reentry corridor

[0110] (2.1) Planning angle of attack profile

[0111] Referring to the typical angle of attack profile design, the angle of attack-speed profile ( ), the angle of attack profile does not change during the reentry process, and the real-time angle of attack can be directly obtained through the velocity. Taking into account the constraints of the reentry process, the designed angle of attack-velocity profile is as follows Figure 5 shown.

[0112] Equation (7) gives the expression of the angle of attack-velocity profile.

[0113] (7)

[0114] After the angle of attack-velocity profile is determined, the process constraints can be directly converted into the boundaries of the reentry corridor. As long as the aircraft flies within the range of the reentry corridor, all constraints are met.

[0115] (2.2) Calculation of reentry corridor

[0116] The process constraints are mainly divided into heat flow constraints, overload constraints, dynamic pressure constraints and balanced gliding constraints when the roll angle is set. The expression is

[0117] (8)

[0118] (9)

[0119] (10)

[0120] (11)

[0121] By converting equations (8) to (11) and the atmospheric density model into mathematical form, we obtain the expression of the height-velocity reentry corridor.

[0122] (12)

[0123] (13)

[0124] Figure 6 A schematic diagram of the altitude-velocity reentry corridor is given. , and Together they form the lower boundary of the reentry corridor; It forms the upper boundary of the high velocity corridor.

[0125] If the height-speed reentry corridor is known, the drag acceleration-speed reentry corridor can be calculated, that is,

[0126] (14)

[0127] (15)

[0128] Figure 7 A schematic diagram of the drag acceleration-velocity reentry corridor is given. , and Together they form the lower boundary of the height-velocity corridor; It forms the upper boundary of the height-speed corridor.

[0129] (3) Propose an optimization problem

[0130] The reference drag acceleration-velocity profile is generally composed of four sections: temperature control section, balanced gliding section, constant drag section and transition section. This design is composed of four sections of linear functions. The physical process and drag acceleration-velocity expression of each stage are given below.

[0131] 1) Temperature control section: ensure that the spacecraft safely passes through the high-temperature heat flow area, starting from the end of the pre-reentry section to the maximum heat flow peak.

[0132] (16)

[0133] 2) Balanced gliding phase: Establish a large overload as quickly as possible.

[0134] (17)

[0135] 3) Constant resistance section: Maintain a large overload to give the aircraft better maneuverability.

[0136] (18)

[0137] 4) Transition section: By adjusting the speed corresponding to the beginning of the transition section, the entire flight can be adjusted more accurately.

[0138] (19)

[0139] Where: , , , , , , These are trajectory design parameters.

[0140] By connecting the trajectories of each stage, the reference resistance acceleration-velocity profile is as follows Figure 8 1 is the temperature control section, 2 is the balanced gliding section, 3 is the constant resistance section, and 4 is the transition section. , , , are the initial resistance acceleration and velocity of each segment respectively; are the terminal resistance acceleration and velocity.

[0141] According to equations (16) to (19), there are 7 design parameters for the reference profile. Direct optimization will be very complicated and not intuitive. If the endpoints of each segment are considered as parameters to be optimized, the optimization difficulty will be greatly simplified.

[0142] Formula (16) - Formula (19) is converted into the following form

[0143] (20)

[0144] A = [ a 1 Da 3 − a 1 V 2 + a 1 V 1 − Da 1 V 3 − V 2 0 Da f − Da 4 V f − V 4 ] (twenty one)

[0145] B = [ − a 1 V 1 + Da 1 − Da 3 − a 1 V 2 + a 1 V 1 − Da 1 V 3 − V 2 V 2 + a 1 V 2 − a 1 V 1 + Da 1 Da 3 − Da f − Da 4 V f − V 4 V f + Da f ] (twenty two)

[0146] Where: , The initial velocity and terminal velocity are known respectively; the resistance acceleration in the constant resistance section remains unchanged, that is .

[0147] Let the optimization parameters for

[0148] X = [ a 1 , V 2 , V 3 , Da 3 , V 4 ] (twenty three)

[0149] From formula (20) and Figure 8 It is found that the 7 parameters can be obtained by The only certainty: given and Through formula (16), we can find and ; Given , , and Through formula (17), we can find , , and ; Given and Through formula (19), we can find and .

[0150] According to the flight mission and characteristics of each segment, the basic principles for optimizing parameter adjustment are:

[0151] 1) The endpoint and the drag acceleration profile determined by the endpoint must be within the reentry corridor;

[0152] 2) Ensure smooth connections at the endpoints where each segment meets another;

[0153] 3) The control variables such as the roll angle do not exceed the control variable constraints;

[0154] 4) Speed ​​at the end of temperature control section It should be after the speed corresponding to the maximum heat flow;

[0155] 5) Considering the uncertainty of flight aerodynamic parameters, the drag acceleration profile should maintain a certain distance from the reentry corridor;

[0156] 6) To ensure the maneuverability of the aircraft, the reference profile should be far away from the lower boundary of the reentry corridor.

[0157] From formula (14), we can see that the upper and lower boundaries of the reentry corridor are The functions are

[0158] (twenty four)

[0159] (25)

[0160] The reference drag acceleration profile parameters are given by The only certainty can be described as and The function of

[0161] (26)

[0162] Taking into account the above principles and the characteristics of the reentry process, the final optimization problem is: given the initial conditions, under the conditions of satisfying the process constraints, control quantity constraints and terminal constraints, the trajectory parameters are optimized with the goal of reducing the total heating of the process and increasing the range, and establishing as large an overload as possible to ensure the maneuverability of the aircraft. , thus planning an optimal trajectory, which is expressed as

[0163] (27)

[0164] (4) Trajectory optimization based on particle swarm algorithm (such as Fig. 9 shown)

[0165] (4.1) Planning the activity range of particle swarm

[0166] The activity range of the particle swarm is planned according to formula (28).

[0167] (28)

[0168] (4.2) Initialize particle swarm

[0169] Initialize the particle swarm parameters and generate The position of the particle and speed , forming a particle swarm , whose expression is:

[0170] (29)

[0171] (30)

[0172] (4.3) Calculate the fitness of the initial particles

[0173] Through the three-degree-of-freedom dynamics model, the numerical integration process obtains the target optimization function value of each particle, that is, the fitness .

[0174] (4.4) Individual optimization and global optimization

[0175] By comparing the fitness of each particle , find the optimal individual fitness of the initial particle swarm and optimal location The optimal fitness of the group and optimal location .

[0176] (4.5) Update particle swarm and cross-bounds check

[0177] 1) Expressions for particle swarm velocity update and out-of-bounds check

[0178] (31)

[0179] (32)

[0180] (33)

[0181] Where: and for [ 0 , 1 ] A random number between For particles In the The individual's best historical position in the iteration; For the group The historical best position of the group in iterations; For particle swarm The inertia weight of the iteration; is the initial inertia weight; is the inertia weight when the iteration reaches the maximum number; is the maximum number of iterations; For particles No. In the iteration The speed of the dimension; and For particle swarm No. The maximum and minimum speed of the dimension.

[0182] 2) Expressions for particle swarm position update and out-of-bounds check

[0183] (34)

[0184] (35)

[0185] Where: For particles No. In the iteration The location of the dimension; and For particle swarm No. The positional maximum and minimum values ​​of a dimension.

[0186] (4.6) Calculate the updated particle fitness and optimal value

[0187] Calculate the The fitness of each particle after the iteration update And update the individual optimal fitness and optimal location And the optimal fitness of the group and optimal location .

[0188] (4.7) Calculation termination judgment

[0189] If the optimization target value does not meet the requirements and the number of iterations If the maximum value is not reached, return to (4.5); if the requirements are met, end the loop.

[0190] (5) Trajectory parameter calculation

[0191] The trajectory parameter solution algorithm process is as follows Fig.10 As shown in the figure, the drag acceleration-velocity profile can be converted into the altitude-velocity profile, and the range and its rate of change, sinking rate, trajectory angle and its rate of change, roll angle, speed change rate, etc. at the corresponding speed can be obtained, and the aircraft state quantity at the corresponding time can be obtained. x = [ H , V , γ ] T and control volume .

[0192] (5.1) Height corresponding to drag acceleration

[0193] According to equation (15) and the atmospheric density calculation formula, the drag acceleration and altitude can be linked to each other, thus solving for the altitude:

[0194] (36)

[0195] (5.2) Range estimation and sinking rate

[0196] By decomposing the aircraft velocity in the horizontal direction, its horizontal velocity and vertical velocity in the direction of gravity acceleration can be obtained, which correspond to the range change rate and altitude sinking rate, respectively, that is,

[0197] (37)

[0198] (38)

[0199] By integrating the rate of change of the range, we can get the range Estimated formula:

[0200] (39)

[0201] (5.3) Trajectory angle and trajectory angle change rate

[0202] By transforming equation (2) and equation (38), the calculation formula of trajectory angle can be obtained:

[0203] (40)

[0204] When the height-velocity profile is known, the height can be expressed as a function related to the velocity. By taking the velocity difference, the corresponding velocity can be obtained. , then the trajectory angle can be solved, that is,

[0205] (41)

[0206] Where: for The corresponding height; for The corresponding height; is the change in velocity differential.

[0207] Knowing the trajectory angle, the corresponding velocity can be obtained by using the difference method again. ,Right now

[0208] (42)

[0209] Where: for The corresponding height; for The corresponding height.

[0210] According to the chain rule, combined with equation (2), we can get the rate of change of the trajectory angle, which is:

[0211] (43)

[0212] (5.4) Roll angle

[0213] Knowing the rate of change of the trajectory angle, according to formula (3), the size of the roll angle can be expressed as

[0214] (44)

[0215] (6) Comparison between trajectory optimization method based on particle swarm algorithm and traditional sequential quadratic programming

[0216] The basic parameters, initial state, terminal state and constraints of the sample vehicle during reentry are set as shown in Tables 1 and 2.

[0217] Table 1 Simulation parameters

[0218]

[0219] Table 2 Control quantity constraints

[0220]

[0221] The above simulation conditions are brought into the trajectory optimization algorithm and the dynamics model for simulation calculation. The optimization results are shown in Table 3 below. Fig.11 The activity range of the particle swarm is given to ensure that the particle swarm will not exceed the reentry corridor, and the optimized trajectory will be far away from the boundary of the reentry corridor; Fig.12 The particle swarm iteration convergence situation is given. The particle swarm optimization algorithm designed in this paper has good global optimization ability and gradually converges to a fixed value without falling into the local optimum, which further verifies the feasibility of the algorithm. Fig.13 The convergence of the global optimal fitness with the number of iterations is given. When , it basically converges and has the global optimal fitness; Fig.14 The optimized drag-acceleration profile results are given. The temperature control section and the balanced gliding section are smoothly connected, meeting the basic requirements of the optimization parameters. The optimized reference drag-acceleration profile is at a certain distance from the reentry corridor and far away from the lower boundary of the corridor, ensuring that the flight trajectory has a certain stability and good maneuverability.

[0222] Table 3 Optimization results

[0223]

[0224] To enhance the contrast, the traditional sequential quadratic programming method is used to optimize the parameters By taking different initial values ​​multiple times, different optimization results can be obtained, and compared with the optimization results based on the particle swarm algorithm, as shown in Table 4. It can be seen from the table that the optimization parameters of sequential quadratic programming are sensitive to the initial values. Different initial values ​​will result in large differences in convergence speed and optimization results, and even fall into local convergence, making it difficult to obtain the optimal solution in the end. The particle swarm algorithm can avoid taking initial values ​​and introduce individual and group cognitive terms, which makes the algorithm have good global optimization ability and fast convergence speed.

[0225] Table 4 Comparison of particle swarm optimization and sequential quadratic programming results

[0226]

Claims

1. A reusable aircraft reentry trajectory optimization method based on particle swarm algorithm, characterized in that: Here are the steps: Step (1) Mathematical modeling and characteristic analysis of reusable aircraft (1.1) Establishment of three-degree-of-freedom dynamic model The three-degree-of-freedom particle dynamics equation considering the earth's rotation is established in the velocity coordinate system: Where: r is the distance from the center of the earth; λ is the latitude; l is the longitude; V is the speed of the aircraft relative to the earth; m is the mass of the aircraft; γ is the track inclination; φ is the track azimuth; is the gravitational acceleration, where g0 is the standard gravitational acceleration on the surface; σ is the tilt angle; ω e is the angular velocity of the earth's rotation; D and L are the drag and lift of the aircraft respectively; To simplify the trajectory calculation, the centrifugal force caused by the rotation of the earth is ignored, and the two-dimensional particle dynamics equation is obtained: (1.2) Establishing the aerodynamic model of the aircraft Aerodynamics can be decomposed in the velocity coordinate system to obtain lift L and drag D, and the calculation formula is: L=ρV 2 C L S / 2 D=ρV 2 C D S / 2 Where: ρ = ρ0e -H / hs is the atmospheric density, where ρ0 is the standard atmospheric density, H is the current height of the aircraft from the ground, hs is the atmospheric normalized height; S is the wing reference area; C L and C D are lift coefficient and drag coefficient respectively; By fitting the lift-drag characteristics of the aircraft, the C L and C D Aerodynamic coefficient fitting model; Where: C L0 , C L1 , C L2 , C L3 , C D0 , C D1 , C D2 , C D3 are aerodynamic fitting coefficients; α is the angle of attack of the aircraft; (1.3) Analysis of aircraft aerodynamic characteristics By observing the changes of lift, drag and lift-to-drag ratio with angle of attack at the same altitude and different speeds, the aerodynamic characteristics of the aircraft are analyzed to determine the angle of attack corresponding to the maximum lift-to-drag ratio; Step (2) Re-enter the corridor (2.1) Planning the angle of attack profile Refer to the typical angle of attack profile design, analyze the impact of the angle of attack on thermal flow constraints, overload constraints, and dynamic pressure constraints, select the angle of attack and the timing of downward pressure, and thus determine the size and shape of the angle of attack profile; (2.2) Calculation of reentry corridor The reentry process constraints are mainly divided into heat flow constraints, overload constraints, dynamic pressure constraints, and balanced gliding constraints when the roll angle is zero. The expression is: n z =(Lcosα+Dsinα) / mg≤n z·max qbar=ρV 2 / 2≤qbar max gV 2 / rL≤0 Where: and are heat flow rate and maximum heat flow rate respectively; K Q is the heat transfer coefficient; n z and n z·max They are pneumatic overload and maximum pneumatic overload; qbar and qbar max are dynamic pressure and maximum dynamic pressure respectively; Through mathematical transformation of process constraints and atmospheric density expressions, multiple process constraints are converted into constraints of a single physical quantity. The three process constraints of heat flow, overload and dynamic pressure constitute the hard constraints of the reentry corridor, and the balanced gliding constraint constitutes the soft constraints of the reentry corridor, which together form the upper and lower boundaries of the reentry corridor. The height-velocity corridor HV and the drag acceleration-velocity corridor Da-V are calculated; the expressions of the height-velocity and drag acceleration-velocity reentry corridors are: Where: Heat flow constraint The corresponding height boundary; H min (n z ) is the overload constraint n zmax The corresponding height boundary; H min (qbar) is the dynamic pressure constraint qbar max The corresponding height boundary; H max (eg) is the height boundary corresponding to the equilibrium gliding constraint; Heat flow constraint The corresponding drag acceleration boundary; Da max (n z ) is the overload constraint n zmax The corresponding drag acceleration boundary; Da max (qbar) is the dynamic pressure constraint qbar max The corresponding drag acceleration boundary; Da min (eg) is the drag acceleration boundary corresponding to the equilibrium gliding constraint; Step (3) Propose an optimization problem The reference resistance acceleration-velocity profile consists of four sections: temperature control section, balanced gliding section, constant resistance section, and transition section. It is composed of four linear functions. The endpoints of each section are considered as parameters to be optimized. The optimization parameters of the reference profile are simplified, and the five parameters to be optimized are obtained as follows: X=[a1,V2,V3,Da3,V4] Where: X is the position; a1 is the linear function slope of the temperature control section; V2 is the terminal velocity of the temperature control section; V3 is the terminal velocity of the balanced gliding section; Da3 is the terminal resistance acceleration of the balanced gliding section; V4 is the terminal resistance acceleration of the constant resistance section; From the formula of the reentry corridor in step (2), we can know that the Da-V functions of the upper and lower boundaries of the reentry corridor are: Yes up =f up (V) Yes down =f down (V) The reference drag acceleration profile parameter is uniquely determined by X and described as a function of X and V, that is, Da = f(X, V) According to the principle of continuity of flight mission and reference profile, the basic principle of optimization of remaining parameters is proposed. The final optimization problem is: given the initial conditions, under the conditions of satisfying flight constraints, the trajectory parameters are optimized to reduce the total heating of the process and increase the range, and to establish as large an overload as possible to ensure the maneuverability of the aircraft, so as to plan an optimal trajectory. The expression of the optimization objective function is: Where: J(X) is the target optimization function, that is, the fitness of the particle swarm algorithm; t0 and t f Da is the initial and final time of the spacecraft reentry; up and Da down are the upper and lower boundaries of the drag acceleration reentry corridor; is the total flight distance; c1, c2, c3 are the weight coefficients of total heating, total flight distance and overload respectively; |a1| is the absolute value of the slope of the temperature control section; |a1| min and |a1| max are the absolute values ​​of the minimum and maximum slopes that the aircraft can achieve in the temperature control section without exceeding the control quantity constraint; V Qm and V Qn are the speeds corresponding to the maximum heat flux decreasing by m% and n% respectively; V alpha↓ is the speed corresponding to the downward pressure of the attack angle; Da2 is the end resistance acceleration of the temperature control section; is the maximum value that Da3 can take; Da3 corresponds to the smaller speed in the upper corridor of reentry; ε is the threshold; Step (4) Trajectory optimization based on particle swarm algorithm (4.1) Planning the activity range of particle swarm According to the constraints of the optimization parameters in step (3), plan the activity range of the particle swarm; (4.2) Initialize particle swarm Initialize the particle swarm parameters, i.e., the number of particles n; particle dimension d; number of iterations k; inertia weight ω; individual learning factor m1; group learning factor m2; maximum position X max and maximum speed V max ; and randomly generate the position X and velocity V of n particles within the activity range of the particle swarm to form a particle swarm N, whose expression is: (4.3) Calculate the fitness of the initial particles The position values ​​of each particle in the initial particle group will constitute a corresponding independent drag acceleration profile. Through the three-degree-of-freedom dynamic model of step (1), the target optimization function value of each particle in the initialization group is obtained by numerical integration, that is, the fitness J(X i ); (4.4) Individual optimization and global optimization By comparing the fitness of each particle J(X i ), and obtain the individual minimum fitness f of the initial particle swarm Pbest,i and the optimal position X Pbest,i With the minimum fitness of the group f Gbest and the optimal position X Gbest ; (4.5) Update particle swarm and cross-boundary check Combining individual optimality and group optimality, the expressions for updating speed and position are given. In order to prevent the updated particle speed and position from exceeding the constraint range, an out-of-bounds check is required. (4.6) Calculate the updated particle fitness and optimize The updated position values ​​of the particle group will form a corresponding independent resistance acceleration profile, and pass through the three-degree-of-freedom dynamic model of step (1) again to obtain the fitness J (X) of each particle after the k+1th update i )(k+1) and update the individual minimum fitness f Pbest,i (k+1) and the optimal position X Pbest,i (k+1) and the minimum fitness of the group f Gbest (k+1) and the optimal position X Gbest (k+1); (4.7) Calculation termination judgment If the optimization target value does not meet the requirements and the number of iterations k does not reach the maximum, return to (4.5); if it meets the requirements, end the loop; Step (5) Trajectory parameter calculation The optimized drag acceleration profile is obtained through step (4). Combined with the three-degree-of-freedom dynamic model in step (1), the drag acceleration-velocity profile is converted into the altitude-velocity profile and the range and its rate of change, sink rate, trajectory angle and its rate of change, roll angle, and velocity change rate at the corresponding speed through the trajectory parameter solution algorithm to obtain the aircraft state x=[H, V, γ] at the corresponding time. T and the control quantity σ, as shown below:

Images