A combined power boost reusable spacecraft climbing trajectory optimization method

Abstract

The present application belongs to the field of optimal control technology, and particularly relates to a combined power boost reusable spacecraft climbing segment trajectory optimization method, aiming at the problem that in the existing trajectory optimization of the climbing segment of the reusable spacecraft, the processing of the thrust is idealized, and the speed improvement effect and the height improvement effect are not matched; on the basis of the classical Gauss pseudospectral method, in order to improve the speed of solving the non-smooth problem, in the polynomial degree and collocation point correction, the adaptive decision method is introduced based on the curvature of the control variable on the time period, the traditional Gauss pseudospectral method is combined with the finite element method, and the method of increasing the interpolation polynomial degree or increasing the time interval is used for correction. Under the premise that the turbojet and ram work together, the distribution mode of the action time and the thrust size of the two is proposed, and the distribution scheme of the combined thrust is given.

Description

Technical Field

[0001] This invention belongs to the field of optimal control technology, specifically relating to a method for optimizing the trajectory of a reusable spacecraft during its ascent phase with combined propulsion. Background Technology

[0002] Reusable spacecraft significantly improve the economy and efficiency of space missions. With the development of space engine technology, combined-fuel spacecraft have further enhanced the flight capabilities of traditional spacecraft. Due to their higher maximum altitude and more efficient propulsion systems, combined-fuel spacecraft have become a key focus of spacecraft technology development for various countries. Considering the complex flight environment and wide range of altitudes and speeds, designing and optimizing their flight trajectories is essential for mission success. In 2013, a company announced the SR-72 hypersonic vehicle using a parallel TBCC engine. Similar in appearance to the SR-71 reconnaissance aircraft, it can reach a maximum speed of Mach 6. This project is a continuation of the FALCON program. During its development, it resolved a series of problems, achieved the technology of two engines sharing an air intake and nozzle, and enabled the turbojet engine to operate at speeds above Mach 2.5.

[0003] Existing technologies for integrating ramjet and turbojet engines into trajectory optimization problems have some limitations. Since direct methods offer significant advantages over indirect methods in solving trajectory optimization problems, scholars both domestically and internationally have conducted extensive research on such approaches. Gao studied the periodic cruise trajectory optimization problem, using a pseudospectral method to fit relevant parameters and quickly obtain the optimized trajectory. Hou et al. studied the trajectory optimization problem for the climb-cruise phase of air-breathing vehicles, transforming it into a convex optimization problem with solution efficiency meeting the requirements of online trajectory planning. Elnagar et al. used the Legendre pseudospectral method to solve the ballistic optimization problem and verified that this method can satisfy various constraints and has high accuracy. Rao et al. also used this method to solve the trajectory optimization problem for reentry vehicles. Zhang et al. improved the Gauss pseudospectral method, proposing the whale optimization algorithm, which solves the trajectory optimization problem for hypersonic vehicles under no-no-fly zone conditions and overcomes the sensitivity of the traditional Gauss pseudospectral method to initial conditions.

[0004] Existing trajectory optimization studies have limited research on trajectory optimization for the ascent phase of reusable spacecraft. In studies on trajectory optimization under combined thrust, the handling of thrust is too idealistic, and the effect of velocity improvement does not match the effect of altitude improvement. In addition, it is also necessary to meet the requirement of spacecraft completing the ascent mission quickly under various constraints.

[0005] In summary, existing trajectory optimization methods for the ascent phase of reusable spacecraft powered by a combination of ramjet and turbojet engines suffer from problems such as an overly idealized handling of thrust and a mismatch between the speed and altitude gains. This results in low accuracy of navigation based on the optimized ascent phase trajectory. Summary of the Invention

[0006] The purpose of this invention is to address the problems in existing trajectory optimization methods for the ascent phase of reusable spacecraft, such as idealized thrust handling and a mismatch between velocity and altitude gains, leading to low navigation accuracy based on the optimized ascent phase trajectory. We propose a combined-power-assisted trajectory optimization method for the ascent phase of reusable spacecraft. This method includes:

[0007] Step S1: Establish 6 coordinate systems and 5 coordinate system transformation relationships;

[0008] The six coordinate systems include:

[0009] geocentric inertial coordinate system o e -x I y I z I Navigation coordinate system o n -x n y n z n Geographic coordinate system T -x T y T z T , aircraft body coordinate system o b -x b y b z b Velocity coordinate system o v -x v y v z v Half-velocity coordinate system o h -x h y h z h ;

[0010] The five coordinate system transformation relationships include: the transformation relationship between the navigation coordinate system and the aircraft body coordinate system. Transformation relationship between the aircraft body coordinate system and the velocity coordinate system Transformation relationship between half-velocity coordinate system and velocity coordinate system Transformation between geographic coordinate system and semi-velocity coordinate system Conversion between navigation coordinate system and geographic coordinate system

[0011] Step S2: Determine the spacecraft configuration of the reusable spacecraft, based on the transformation relationships between the geocentric inertial coordinate system, navigation coordinate system, geographic coordinate system, half-velocity coordinate system, and navigation coordinate system and geographic coordinate system. Transformation between geographic coordinate system and semi-velocity coordinate system A three-degree-of-freedom dynamic model was established for the aircraft configuration;

[0012] Step S3: Obtain existing aircraft data, and determine the aircraft's power engine parameters and aerodynamic parameters based on the existing aircraft data;

[0013] Step S4: Based on the three-degree-of-freedom model of the spacecraft, the parameters of the spacecraft's power engine and aerodynamic parameters, the transformation relationship between the spacecraft's body coordinate system, velocity coordinate system, half-velocity coordinate system, the transformation relationship between the spacecraft's body coordinate system and velocity coordinate system, and the transformation relationship between the half-velocity coordinate system and velocity coordinate system, construct the trajectory optimization problem for the climb phase of a reusable spacecraft.

[0014] Step S5: Use the hp adaptive Gauss pseudospectral method to optimize the climb trajectory of the reusable spacecraft and obtain the optimal climb trajectory of the reusable spacecraft.

[0015] The beneficial effects of this invention are as follows:

[0016] Building upon the classical Gaussian pseudospectral method, to improve the solution speed for non-smooth problems, an adaptive decision-making method is introduced based on the curvature of the control variables over a time interval, combining the traditional Gaussian pseudospectral method and the finite element method, when correcting polynomial degree and collocation points. The decision is made by either increasing the interpolation polynomial degree or extending the time interval. Under the premise of turbojet and ramjet operating simultaneously, a method for allocating their interaction time and thrust magnitude is proposed, and a scheme for distributing the combined thrust is presented. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the calculation process of the hp adaptive Gauss pseudospectral method of the present invention;

[0018] Figure 2 This is a flowchart of the spacecraft climb trajectory calculation process of the present invention;

[0019] Figure 3 This is a schematic diagram of the simulated velocity change curve of the present invention;

[0020] Figure 4 This is a schematic diagram of the geocentric distance variation curve of the present invention;

[0021] Figure 5 This is a schematic diagram of the mass change curve of the present invention. Detailed Implementation

[0022] Specific implementation method one: Combining Figures 1 to 5 This invention describes the following:

[0023] Step S1: Establish 6 coordinate systems and 5 coordinate system transformation relationships;

[0024] The six coordinate systems include:

[0025] geocentric inertial coordinate system o e -x I y I z I Navigation coordinate system o n -x n y n z n Geographic coordinate system T -x T y T z T , aircraft body coordinate system o b -x b y b z b Velocity coordinate system o v -x v y v z v Half-velocity coordinate system o h -x h y h z h ;

[0026] The five coordinate system transformation relationships include: the transformation relationship between the navigation coordinate system and the aircraft body coordinate system. Transformation relationship between the aircraft body coordinate system and the velocity coordinate system Transformation relationship between half-velocity coordinate system and velocity coordinate system Transformation between geographic coordinate system and semi-velocity coordinate system Conversion between navigation coordinate system and geographic coordinate system

[0027] Step S2: Determine the spacecraft configuration of the reusable spacecraft, based on the transformation relationships between the geocentric inertial coordinate system, navigation coordinate system, geographic coordinate system, half-velocity coordinate system, and navigation coordinate system and geographic coordinate system. Transformation between geographic coordinate system and semi-velocity coordinate system A three-degree-of-freedom dynamic model was established for the aircraft configuration;

[0028] Step S3: Obtain existing aircraft data, and determine the aircraft's power engine parameters and aerodynamic parameters based on the existing aircraft data;

[0029] Step S4: Based on the three-degree-of-freedom model of the spacecraft, the parameters of the spacecraft's power engine and aerodynamic parameters, the transformation relationship between the spacecraft's body coordinate system, velocity coordinate system, half-velocity coordinate system, the transformation relationship between the spacecraft's body coordinate system and velocity coordinate system, and the transformation relationship between the half-velocity coordinate system and velocity coordinate system, construct the trajectory optimization problem for the climb phase of a reusable spacecraft.

[0030] Step S5: Use the hp adaptive Gauss pseudospectral method to optimize the climb trajectory optimization problem of the reusable spacecraft and obtain the optimal climb trajectory of the reusable spacecraft.

[0031] This invention analyzes flight constraints and selects performance indicators based on the characteristics of the climb phase. The parameters of the turbojet / ramjet combined-engine are modeled to ensure that the power system parameters meet the thrust requirements of the aircraft across a wide range of flight states. In the climb phase trajectory optimization, different engine operating modes are selected, and the hp adaptive Gauss pseudospectral method is used to optimize the climb phase under different engine operating modes.

[0032] Specific Implementation Method Two: The difference between this implementation method and Specific Implementation Method One is that...

[0033] The construction of the geocentric inertial coordinate system o e -x I y I z I The specific process is as follows:

[0034] The origin of the geocentric inertial coordinate system is the center of the Earth, denoted as o. e The coordinate axes o of the geocentric inertial coordinate system e x I， coordinate axis o e y I， coordinate axis o e z I Construct a right-handed rectangular coordinate system;

[0035] The coordinate axes o of the geocentric inertial coordinate system e x I In the equatorial plane, the vernal equinox is pointed to at a fixed time.

[0036] The specified fixed time is usually taken as 12:00 on January 1, 2000.

[0037] The coordinate axes o of the geocentric inertial coordinate system e z IPerpendicular to the equatorial plane; pointing towards the North Pole.

[0038] The construction of the navigation coordinate system o n -x n y n z n The specific process is as follows:

[0039] The navigation coordinate system, also known as the North-Sky-East coordinate system, is a fixed Earth coordinate system.

[0040] Using the launch point of the aircraft as the navigation coordinate system n -x n y n z n The origin of the coordinate system is denoted as o. n ,

[0041] The coordinate axes o of the navigation coordinate system n z n , coordinate axis o n x n , coordinate axis o n y n Construct a right-handed rectangular coordinate system;

[0042] The coordinate axes o of the navigation coordinate system n z n Within the line connecting the Earth's center and the spacecraft's launch point, pointing upwards,

[0043] The coordinate axes o of the navigation coordinate system n x n Within the meridian plane where the aircraft is located, it points north;

[0044] The construction of the geographic coordinate system T -x T y T z T The specific process is as follows:

[0045] A geographic coordinate system is constructed using the intersection of the line connecting the Earth's center and the spacecraft's center of mass with the elliptical surface of the Earth. T -x T y T z T The origin o t ,

[0046] coordinate axes of a geographic coordinate system t z t , coordinate axis o t x t , coordinate axis o t y t Construct a right-handed rectangular coordinate system;

[0047] coordinate axes of a geographic coordinate system ty t The line connecting the Earth's center and the spacecraft's center of mass coincides with the line pointing upwards.

[0048] coordinate axes of a geographic coordinate system t x t Within the meridian plane where the aircraft is located, pointing north,

[0049] The construction of the aircraft body coordinate system o b -x b y b z b The specific process is as follows:

[0050] With the spacecraft's center of mass as the origin of the spacecraft's body coordinate system, denoted as o. b ,

[0051] The coordinate axes of the aircraft body coordinate system b z b , coordinate axis o b x b and coordinate axis o b y b Construct a right-handed rectangular coordinate system;

[0052] The coordinate axes of the aircraft body coordinate system b x b Aligned with the longitudinal axis of the aircraft's fuselage, pointing towards the head,

[0053] The coordinate axes of the aircraft body coordinate system b y b Located within the longitudinal symmetry plane of the aircraft's fuselage, pointing upwards;

[0054] The constructed velocity coordinate system o v -x v y v z v The specific process is as follows:

[0055] With the spacecraft's center of mass as the origin of the velocity coordinate system, denoted as o. v ,

[0056] The coordinate axes o of the velocity coordinate system v z v , coordinate axis o v x v and coordinate axis o v y v Construct a right-handed rectangular coordinate system;

[0057] The coordinate axes o of the velocity coordinate system v x v Coinciding with the direction of the aircraft's velocity vector;

[0058] The coordinate axes o of the velocity coordinate systemv y v Located vertically within the longitudinal symmetry plane of the aircraft, pointing upwards;

[0059] The construction of the half-velocity coordinate system o h -x h y h z h The specific process is as follows:

[0060] With the spacecraft's center of mass as the origin of the half-velocity coordinate system, denoted as o h ,

[0061] The coordinate axis o of the half-velocity coordinate system h z h , coordinate axis o h x h and coordinate axis o h y h Construct a right-handed rectangular coordinate system;

[0062] The coordinate axis o of the half-velocity coordinate system h x h The axis coincides with the direction of the aircraft's velocity vector;

[0063] The coordinate axis o of the half-velocity coordinate system h y h Located in the vertical plane containing the velocity vector, pointing upwards;

[0064] The transformation relationship between the navigation coordinate system and the aircraft body coordinate system Expressed as a formula:

[0065] Define pitch angle Yaw angle ψ, roll angle γ, establish the transformation relationship between the navigation coordinate system and the aircraft body coordinate system:

[0066]

[0067] The transformation relationship between the aircraft body coordinate system and the velocity coordinate system is constructed. The specific process is as follows:

[0068] Define the angle of attack α and the sideslip angle β, and establish the transformation relationship between the aircraft body coordinate system and the velocity coordinate system as follows:

[0069]

[0070] The transformation relationship between the constructed half-velocity coordinate system and the velocity coordinate system is described. The specific process is as follows:

[0071] Define the velocity tilt angle γ ν The transformation relationship between the half-velocity coordinate system and the velocity coordinate system is established as follows:

[0072]

[0073] The transformation relationship between the geographic coordinate system and the semi-velocity coordinate system is constructed. The specific process is as follows:

[0074] Define the velocity deflection angle σ and the velocity tilt angle θ, then the transformation relationship between the geographic coordinate system and the semi-velocity coordinate system is as follows:

[0075]

[0076] The relationship between the navigation coordinate system and the geographic coordinate system is constructed. The specific process is as follows:

[0077] Define the geocentric latitude difference between the navigation coordinate system and the geographic coordinate system as Δφ1, and the longitude difference as Δλ1. Then the transformation relationship between the two is:

[0078]

[0079] The other steps and parameters are the same as in Specific Implementation Method 1.

[0080] Specific Implementation Method Three: The difference between this implementation method and Specific Implementation Method One is that...

[0081] In step S2, the spacecraft configuration of the reusable spacecraft is determined based on the transformation relationships between the geocentric inertial coordinate system, navigation coordinate system, geographic coordinate system, half-velocity coordinate system, and the navigation coordinate system and geographic coordinate system. Transformation between geographic coordinate system and semi-velocity coordinate system A three-degree-of-freedom dynamic model of the aircraft configuration is established; the specific process is as follows:

[0082] S2.1: Determine the aircraft configuration as Lockheed Martin SR-72 and construct a dynamic model of the aircraft in a geocentric inertial coordinate system;

[0083] S2.2: Construct the dynamic model of the aircraft in the navigation coordinate system based on the dynamic model of the aircraft in the geocentric inertial frame:

[0084] S2.3: Based on the transformation relationship between the navigation coordinate system and the geographic coordinate system, and the transformation relationship between the geographic coordinate system and the semi-velocity coordinate system,

[0085] The dynamic model of the aircraft in the navigation coordinate system is decomposed in the half-velocity system to obtain the dynamic equation of the aircraft's center of mass, the kinematic equation of the aircraft's center of mass, and the equation of the change in the aircraft's mass.

[0086] The aircraft's center of mass dynamic equation, the aircraft's center of mass kinematic equation, and the aircraft's mass change equation are used as a three-degree-of-freedom dynamic model. Other steps and parameters are the same as in one of the specific implementation methods one or two.

[0087] Specific Implementation Method Four: This implementation method differs from Specific Implementation Methods One through Four in that...

[0088] The specific process of constructing the dynamic model of the aircraft in the geocentric inertial coordinate system in S2.1 is as follows:

[0089] According to Newton's second law, the dynamic model of the aircraft in the geocentric inertial frame is as follows:

[0090]

[0091] In the formula, m is the mass of the aircraft, t represents time, and r e denoted as the geocentric distance of the spacecraft's location in the geocentric inertial coordinate system, P is the engine thrust of the spacecraft, R is the aerodynamic force, and g is the gravitational acceleration.

[0092] In step S2.2, a dynamic model of the aircraft in the navigation coordinate system is constructed based on the dynamic model of the aircraft in the geocentric inertial frame. The dynamic equations of the aircraft in the navigation coordinate system are derived. According to the definition of the coordinate system, the vector differential rule is used to obtain the following: The specific process is as follows:

[0093]

[0094] In the formula, ω e The angular velocity of the navigation coordinate system relative to the inertial coordinate system is represented by r. e =r m +R0, r m R0 represents the position vector of the aircraft in the navigation system, and F represents the vector from the Earth's center to the origin of the navigation coordinate system. c F represents the Coriolis inertial force. e Indicates the inertial force involved;

[0095] In step S2.3, based on the transformation relationships between the navigation coordinate system and the geographic coordinate system, and between the geographic coordinate system and the half-velocity coordinate system, the dynamic model of the aircraft in the navigation coordinate system is decomposed in the half-velocity system to obtain the aircraft's center-of-mass dynamic equation, the aircraft's center-of-mass kinematic equation, and the aircraft's mass change equation. These equations are then used as a three-degree-of-freedom dynamic model. The specific process is as follows:

[0096] S2.3.1: Define the aircraft velocity V, the aircraft velocity tilt angle γv, the aircraft velocity deflection angle б, and the aircraft velocity tilt angle θ. Decompose the dynamic model formula (8) of the aircraft in the navigation coordinate system into the half-velocity system. The decomposition process is based on the coordinate transformation relationship between the navigation coordinate system and the half-velocity system to obtain the dynamic equation of the aircraft's center of mass (force relationship), which is expressed by the formula:

[0097]

[0098] In equation (9), Represents the rate of change of velocity V. This represents the rate of change of the velocity inclination angle γv. This represents the rate of change of the velocity deflection angle б.

[0099] P xh This indicates the thrust along the coordinate axis o of the half-velocity coordinate system. h x h The component, P yh This indicates the thrust along the coordinate axis o of the half-velocity coordinate system. h y h The component, P zh This indicates the thrust along the coordinate axis o of the half-velocity coordinate system. h z h The amount,

[0100] ω θ The rotational angular velocity is represented along the coordinate axis o of the half-velocity coordinate system. h y h The component, ω б The rotational angular velocity is represented along the coordinate axis o of the half-velocity coordinate system. h x h The component, ω V The rotational angular velocity is represented along the coordinate axis o of the half-velocity coordinate system. h x h The amount,

[0101] m represents the mass of the aircraft; X represents the drag force of the aircraft; Y represents the lift force of the aircraft; Z represents the lateral force of the aircraft.

[0102] g θ This represents the Earth's gravitational force along the coordinate axis o of the half-velocity coordinate system. h x h The amount, g V This represents the Earth's gravitational force along the coordinate axis o of the half-velocity coordinate system. h y h The amount, g б This represents the Earth's gravitational force along the coordinate axis o of the half-velocity coordinate system. h z h The components, where g r G represents the component of gravitational acceleration along the Earth's center. ω This represents the component of gravitational acceleration along the direction of Earth's rotation;

[0103] r represents the Earth's radius, α e The flatness of the Earth is represented by φ, latitude by J, and zone harmonic coefficient by μ.E Represents the gravitational constant.

[0104] S2.3.2: Based on the geometric relationships between the variables in the aircraft's center-of-mass dynamic equations, the aircraft's center-of-mass kinematic equations (geometric relationships) are constructed as follows:

[0105]

[0106] In the formula, Represents the rate of change of the geocentric radius vector. Indicates the rate of change of latitude. Indicates the rate of change of longitude;

[0107] S2.3.3: The equation for the change in the mass of the aircraft is constructed as follows:

[0108]

[0109] In the formula, f represents the rate of change of the aircraft's mass. m The function representing the change in mass, where Ma represents the Mach number, h represents altitude, and c represents the engine operating factor, is a time-dependent function. Since fuel mass is constantly decreasing, it is negative. The rate of mass change is related to the Mach number, flight altitude, and engine operating factor; other steps and parameters are the same as in one of the specific implementation methods one to three.

[0110] Specific Implementation Method Five: The difference between this implementation method and Specific Implementation Methods One to Four is that...

[0111] In step S3, existing aircraft data is acquired, and the aircraft's power engine parameters and aerodynamic parameters are determined based on this data. The specific process is as follows:

[0112] S3.1: Select a parallel turbojet / ramjet combined engine as the aircraft's propulsion system.

[0113] The parallel turbojet / ramjet combined-power engine is an aviation propulsion system that combines a turbojet engine and a ramjet engine. It has one turbojet engine and one ramjet engine, which operate in parallel. Its performance is mainly related to flight speed, altitude, and engine duty cycle.

[0114] The fitting formulas for the parameters of the turbojet engine and the ramjet engine are constructed and expressed as follows:

[0115]

[0116] In the formula, P1 represents the thrust of the turbojet engine, P2 represents the thrust of the ramjet engine, dm1 represents the fuel consumption rate per second of the turbojet engine, dm2 represents the fuel consumption rate per second of the ramjet engine, and fp1 f represents the thrust function of a turbojet engine. p2 f represents the thrust function of a ramjet engine. m1 f represents the fuel consumption rate per second function of a turbojet engine. m2 This represents the fuel consumption rate per second function of a ramjet engine.

[0117] S3.2: Obtain existing flight data of turbojet engine and ramjet engine aircraft. Based on this data, obtain the aircraft's propulsion engine parameters, including turbojet engine parameters and ramjet engine parameters, using interpolation fitting methods. The specific process is as follows:

[0118] Based on existing flight data of turbojet engine aircraft and fitting formulas for turbojet engine parameters, the parameters of turbojet-powered engines are determined by interpolation fitting.

[0119] Based on existing data on the flight of ramjet-powered aircraft and the fitting formula for ramjet engine parameters, the parameters of the ramjet-powered engine are determined by interpolation fitting.

[0120] The interpolation fitting method is an important tool in numerical analysis, used to construct a mathematical model or function based on a set of known data points, so that the function can accurately pass through all known points. This invention determines engine parameters using existing aircraft flight data through interpolation fitting.

[0121] The overall parameters of the parallel turbojet / ramjet combined engine are determined based on the parameters of the turbojet engine and the ramjet engine. Then, based on the thrust and fuel consumption rate formulas of the parallel turbojet / ramjet combined engine and the overall parameters of the parallel turbojet / ramjet combined engine, the flight conditions of the parallel turbojet / ramjet combined engine are obtained.

[0122] The flight conditions include: operating altitude h of the combined-propellant engine; range: 5–35 km; Mach number Ma of the combined-propellant engine; range: 0.2–5 Ma.

[0123] The parameters of the turbojet engine and the ramjet engine are determined by interpolation fitting; the specific process is as follows:

[0124] Based on existing flight data, including Mach number, flight altitude, and data on drag coefficient, lift coefficient, and angle of attack, the parameters of the turbojet engine and the ramjet engine are determined using interpolation fitting methods.

[0125] Based on the thrust and fuel consumption rate of the parallel turbojet / ramjet combined engine, as well as the parameters of the turbojet engine and the ramjet engine, the specific process for obtaining the flight conditions of the parallel turbojet / ramjet combined engine is as follows:

[0126] This determines the overall parameters of the combined-power engine. Combining the interpolation results above, the parameters are analyzed to obtain the flight conditions of the parallel turbojet / ramjet combined-power engine, a procedure well-known to those skilled in the art.

[0127] S3.3: Obtain the fixed parameter data of the aircraft; based on the existing fixed parameter data of the aircraft, use an interpolation fitting method to obtain the aerodynamic parameters; the aerodynamic parameters include: aerodynamic force, aerodynamic coefficient, and lift-to-drag ratio coefficient; the specific process is as follows:

[0128] S3.3.1: Construct fitting formulas for aerodynamic forces and aerodynamic coefficients. The specific process is as follows:

[0129]

[0130] In the formula, S represents the area of ​​force application, q represents the dynamic pressure, α represents the angle of attack, β represents the sideslip angle, and f x The drag coefficient function, f y Represents the lift coefficient function, f z Represents the lateral force coefficient function;

[0131] The aerodynamic forces include: drag X, lift Y, and lateral force Z;

[0132] The aerodynamic coefficient includes the drag coefficient C. X Lift coefficient C Y Lateral force coefficient C Z ;

[0133] S3.3.2: Determine the fixed parameters of the aircraft based on the aircraft configuration and obtain the fixed parameter data of the aircraft; the fixed parameters of the aircraft include fixed parameters such as aircraft structural parameters and aircraft environmental parameters;

[0134] Based on the fixed parameter data of the aircraft, the parameters of the aircraft's power engine, and the calculation formula of the aerodynamic coefficient, the aerodynamic coefficient is obtained by fitting through interpolation.

[0135] S3.3.3: Based on the aerodynamic coefficient and its calculation formula, the aerodynamic force is obtained through interpolation fitting, and then the lift-to-drag ratio coefficient is obtained based on the aerodynamic force.

[0136] The lift-to-drag ratio coefficient is known to those skilled in the art based on aerodynamic forces: drag X, lift Y, and lateral force Z. The calculation data of this invention shows that the maximum lift-to-drag ratio of the aircraft can reach more than 7.0, exhibiting good aerodynamic characteristics. Other steps and parameters are the same as in one of the specific embodiments one to four.

[0137] Specific Implementation Method Six: The difference between this implementation method and Specific Implementation Methods One to Five is that...

[0138] In step S4, based on the three-degree-of-freedom model of the spacecraft, the parameter model of the spacecraft's propulsion engine, the aerodynamic parameters of the spacecraft, the transformation relationships between the spacecraft's body coordinate system, velocity coordinate system, half-velocity coordinate system, the transformation relationship between the spacecraft's body coordinate system and velocity coordinate system, and the transformation relationship between the half-velocity coordinate system and velocity coordinate system, a trajectory optimization problem for the climb phase of a reusable spacecraft is constructed; the specific process is as follows:

[0139] S4.1: Select the control variables for the trajectory optimization problem of the climb phase of a reusable spacecraft; based on the transformation relationship between the control variables and the spacecraft's body coordinate system and velocity coordinate system, obtain the transformed three-degree-of-freedom dynamic model of the spacecraft;

[0140] S4.2: Construct a state constraint model for the flight's climb phase based on the transformed three-degree-of-freedom dynamic model of the aircraft, the parameter model of the aircraft's power engine, and the aerodynamic parameters of the aircraft;

[0141] S4.3: Establish performance optimization indices for the spacecraft's climb trajectory; based on the performance optimization indices, control variables, and state constraint model of the spacecraft's climb trajectory, construct the climb trajectory optimization problem for a reusable spacecraft; other steps and parameters are the same as in one of the specific implementation methods one to five.

[0142] Specific Implementation Method Seven: The difference between this implementation method and Specific Implementation Methods One through Six is ​​that...

[0143] The control variables in the climb trajectory optimization problem of the reusable spacecraft in S4.1 include: a first control variable and a second control variable; wherein, the first derivative of the angle of attack is selected as the first control variable, and the first derivative of the velocity tilt angle is selected as the second control variable.

[0144] Based on the transformation relationship between the state variables and the body coordinate system and velocity coordinate system of the aircraft, the transformed three-degree-of-freedom dynamic model of the aircraft is obtained; expressed by the formula:

[0145]

[0146] In the formula, Let α be the first derivative of the angle of attack. For the velocity tilt angle γ v The first derivative, and To control the quantity;

[0147] and For state variables;

[0148] The process of obtaining the transformed three-degree-of-freedom dynamic model of the aircraft based on the transformation relationship between the aircraft body coordinate system and the velocity coordinate system is a well-known process to those skilled in the art.

[0149] In S4.2, a state constraint model for the flight climb phase is constructed based on the transformed three-degree-of-freedom dynamic model of the aircraft, the parameters of the aircraft's power engine, and the aerodynamic parameters of the aircraft.

[0150] The state constraint model for the aircraft's climb phase includes: a turbojet climb phase constraint model, a combined climb phase constraint model, a ramjet climb phase constraint model, and a full-range climb phase constraint model.

[0151] The specific process is as follows:

[0152] The process of constructing the state constraint model of the aircraft's climb phase based on the transformed three-degree-of-freedom dynamic model of the aircraft, the parameters of the aircraft's power engine, and the aerodynamic parameters of the aircraft is well known to those skilled in the art;

[0153] S4.2.1: Construct a constraint model for the turbojet's climb section.

[0154] The constraint model for the turbojet climbing section includes: initial state constraints for the turbojet climbing section and terminal state constraints for the turbojet climbing section.

[0155] The initial state constraints of the turbojet climb section include:

[0156] The initial height of the turbojet's climb phase is hw0 = 10 km, and the initial velocity of the turbojet's climb phase is Vw0 = 408 m / s. 2 The initial velocity inclination angle of the turbojet climb section is θw0 = 0°, the initial velocity deflection angle of the turbojet climb section is σw0 = -21.684°, and the initial mass of the turbojet climb section is mw0 = 12000 kg.

[0157] The terminal state constraints of the turbojet climb section include:

[0158] Terminal velocity Vw of the turbojet climb phase f =544m s, the terminal velocity inclination angle of the turbojet climb section θwf=0°, and the ground center distance increment of the terminal turbojet climb section Δwr1=4500m;

[0159] S4.1.2: Construct a combined climbing section constraint model;

[0160] The combined climbing segment constraint model includes: initial state constraints of the combined power climbing segment and terminal state constraints of the combined power climbing segment;

[0161] The initial state constraints of the combined dynamic climbing section include:

[0162] The initial velocity of the combined power climb segment is Vz0 = 544 m / s, the initial velocity inclination angle of the combined power climb segment is θz0 = 0°, and the initial geocentric distance of the combined power climb segment is r. z =rz0 + Δrz1;

[0163] The combined power climbing section's terminal state constraints include:

[0164] The terminal velocity of the combined power climb section is Vzf = 850 m / s, and the ground center distance increment at the terminal of the combined power climb section is Δzr = 2500 m.

[0165] S4.1.3: Construct a constraint model for the stamping climbing section;

[0166] The constraint model of the stamping climbing section includes: initial state constraints of the stamping climbing section and terminal state constraints of the stamping climbing section;

[0167] The initial state constraints of the stamping climbing section include:

[0168] Initial velocity of the ram-climb section Vc0 = 850 m / s; initial distance from the ground rc of the ram-climb section rc = rc0 + Δrc 1+ Δrc2;

[0169] The terminal state constraints of the stamping climbing section include:

[0170] The terminal speed Vc of the stamping climbing section f =1530m / s, the increase in distance from the center of the earth at the end of the ram-press climbing section Δcr3 = 11500m, and the velocity inclination angle at the end of the ram-press climbing section θc f =0°, angle of attack α at the end of the stamping climb section cf =4.12°;

[0171] S4.1.4: Construct a constraint model for the entire climbing section;

[0172] The constraint model for the entire climbing section includes:

[0173] Maximum overload constraint n throughout the climbing section max =3, maximum dynamic pressure constraint q throughout the climbing section max =100 kPa, maximum heat flux density constraint throughout the climb section The climbing phase has the following constraints: speed tilt angle constraint -0.5° < θ < 10°, angle of attack constraint 1° < α < 10°, speed tilt angle constraint υ = 0°, and control quantity constraints |u1| ≤ 0.1° / s, |u2| ≤ 0° / s; total thrust P and total fuel consumption rate dm during the combined power climbing phase.

[0174] The total thrust P and total fuel consumption rate dm of the combined power climbing phase are expressed by the following formulas:

[0175] P=c1(t)P1(Ma,H)+c2(t)P2(Ma,H)

[0176] dm=c1(t)dm1(Ma,H)+c2(t)dm2(Ma,H) (22)

[0177] Where c1 and c2 are both functions of time t, P1 is the thrust of the ramjet engine, P2 is the thrust of the turbojet engine, dm1 is the fuel consumption rate per second of the ramjet engine, dm2 is the fuel consumption rate per second of the turbojet engine, c1=0.7+0.3(t-t1) / (t2-t1), c2=1-(t-t1) / (t2-t1);

[0178] Simulations were performed based on the established constraint model. The longitude, latitude, longitude, and latitude of the launch point were 0, 30, 15.3, and 53 degrees, respectively, and the initial velocity deflection angle of the launch point was -21.684°.

[0179] The terminal constraints of the climb phase are as follows: initial velocity: 408 m / s; terminal velocity: 1530 m / s; terminal velocity tilt angle: 0 degrees; initial altitude: 10000 m; climb phase ground center increment: 18500 m. Estimate the angle of attack α required for altitude hold flight. f =4.12°;

[0180] The performance optimization index for the aircraft's climb trajectory is established in S4.3; expressed by the formula:

[0181]

[0182] In the formula, R d R is the radius of the aircraft's nose, and C1 is a constant related to the aircraft's characteristics. In this paper, we take R as the radius of the nose. d =0.09m, C1=11000, ρ0 is the standard atmospheric density at sea level, g0 represents the Earth's gravitational pull at sea level, and R0 represents the Earth's average radius. The performance optimization index for this invention is selected as follows: minimizing total heat absorption during the climb phase is used as the optimized performance index. This is combined with the heat flux density calculation formula to obtain the performance optimization index for the aircraft's climb phase trajectory.

[0183] The process of constructing the climb trajectory optimization problem for a reusable spacecraft based on the first and second control variables (performance optimization indices) and the spacecraft climb trajectory constraint model is well known to those skilled in the art.

[0184] In simple terms, the trajectory optimization problem for the climbing section is to find the optimal control variable that maximizes performance and satisfies the constraint model.

[0185] The researchers set their own initial conditions based on the four constraint models, according to the experiments.

[0186] The other steps and parameters are the same as those in any of the specific implementation methods one to six.

[0187] Specific Implementation Method Eight: The difference between this implementation method and Specific Implementation Methods One to Seven is that...

[0188] In step S5, the Gaussian pseudospectral method and sequential quadratic programming algorithm are used to optimize the ascent trajectory of the reusable spacecraft to obtain the optimal ascent trajectory; the specific process is as follows:

[0189] S5.1: Divide the entire time of the aircraft trajectory optimization problem into S... t A time interval;

[0190] S5.2: Determine the initial number of collocations in each time interval; where the initial number of collocations in the s-th time interval is Ms; s∈S t S t Let S be a positive integer; the spacecraft trajectory optimization problem for each of the St time intervals is discretized to obtain S. t NLP problems within a time interval;

[0191] The NLP problem mentioned is a nonlinear programming problem; the solution to nonlinear programming problems is mature, so it can be solved.

[0192] S5.3: Using NLP algorithms to solve S t Solving NLP problems across time intervals yields S t Calculation results for NLP problems in each time interval;

[0193] The approximation error of each time period in the NLP problem calculation result is calculated. The NLP problem refers to a nonlinear programming problem. The solution method for nonlinear programming problems is very mature and can be solved. Therefore, any NLP solution algorithm can be selected. The calculation of the approximation error of each time period in the NLP problem calculation result is also well known to those in the art. The preferred NLP solution algorithm is the sequential quadratic programming algorithm.

[0194] S5.4: Setting Error Calculate the approximation error of the NLP problem calculation results for S time intervals, and determine whether the approximation error of the NLP problem calculation results for all S time intervals is less than the set error.

[0195] If the approximation error of the NLP problem calculation result in the s-th time interval is ε (s) Greater than the set error Then proceed to S5.5;

[0196] If the approximation error of the NLP problem calculation results in the St time intervals is less than the set error Then stop iterating and combine the NLP problem calculation results of the St time intervals to form the optimal climbing trajectory for the entire time interval;

[0197] S5.5: Set the maximum correction parameter r max Calculate the correction parameters for each time interval based on the NLP problem calculation results for the St time intervals; update each time interval according to the correction parameters for each time interval, and return to S5.2 after updating;

[0198] The other steps and parameters are the same as those in any of the specific implementation methods one to seven.

[0199] Specific Implementation Method Nine: The difference between this implementation method and Specific Implementation Methods One through Eight is that...

[0200] In S5.2, the aircraft trajectory optimization problem for a given time interval is discretized to obtain an NLP problem for that time interval. The specific process is as follows:

[0201] S5.2.1: Transform the time interval of the aircraft trajectory optimization problem to the time interval for applying the Gauss pseudospectral method;

[0202] S5.2.2: Obtain K collocation points in the transformed time interval; use the number of collocation points as the degree of the Largerange interpolation polynomial, where K is a positive integer;

[0203] Where the initial number of points is S t The updated number of points is selected according to the updated quantity; here it simply refers to a general number of K points.

[0204] Based on K collocation points, construct a Larger interpolation polynomial, and use the Larger interpolation polynomial to discretize and represent the state and control variables in the time interval of the aircraft trajectory optimization problem;

[0205] S5.2.3: Convert the discretized representation of the state and control variables of the aircraft trajectory optimization problem in time interval into algebraic form to obtain the algebraic equations satisfied by the state variables on the collocation points; construct an NLP problem in time interval based on the algebraic equations satisfied by the state variables on the collocation points;

[0206] Constructing a time-interval NLP problem based on the algebraic equations satisfied by the state variables at collocation points; This is a well-known concept in the field. By using algebraic equations to constrain the evolution or changes of the state variables, certain constraints are set for the optimization problem, thus transforming it into an NLP problem.

[0207] In S5.2.1, the time interval of the aircraft trajectory optimization problem is transformed to the time interval for applying the Gauss pseudospectral method; the specific process is as follows:

[0208] The time interval for the aircraft trajectory optimization problem is [t0,t]. f The time interval [-11] for applying the Gauss pseudospectral method is converted as follows:

[0209]

[0210] t represents the time interval [t0, t] f In the time variable ], t0 represents the initial time, t f Indicates the terminal time; This represents the time variable t transformed into the time variable of the Gauss pseudospectral method;

[0211] In S5.2.2, the state and control variables of the aircraft trajectory optimization problem are discretized using Larger interpolation polynomials, expressed by the following formula:

[0212]

[0213] in, Represents the state variables of the aircraft. Represents the discretized representation of the aircraft state variables. For matching points The state variables of the aircraft For matching points Discretized representation of aircraft state variables Let K be the i-th collocation point, K be the total number of collocation points, and i be the interpolation sequence number. For aircraft control variables, For discretized representation of aircraft control quantities; Represents the Lagrange basis functions; Indicates matching points Aircraft control quantities; Indicates matching points Discretized representation of aircraft control variables; X represents the aircraft state at the terminal moment; f Let f() represent the discretized terminal state variables of the aircraft, and let f() represent the transformed three-degree-of-freedom dynamic model of the aircraft. Indicates Gaussian weights. Locate the k-th Legendre-Gauss point; t f The terminal time is t0, which is the start time.

[0214] In S5.2.3, the state and control variables of the discretized representation of the aircraft trajectory optimization problem over the time interval are converted into algebraic form to obtain the algebraic equations satisfied by the state variables at the collocation points. The specific process is as follows:

[0215]

[0216] Among them, D ki () represents the weighting coefficient. Let k be the k-th collocation point (k = 1, ..., K). express The derivative, Indicates matching points Discretized representation of aircraft control variables; express The derivative of Indicates matching points Aircraft control quantities; Indicates matching points The state variables of the aircraft;

[0217] Weighting coefficient D ki The differential matrix formed Can be determined offline

[0218] In step S5.5, the correction parameter r for the s-th time interval is calculated based on the calculation results of the NLP problem for the s-th time interval. s The specific process is as follows:

[0219]

[0220] K max (s) For the time interval s, K (s) The maximum value obtained by (τ) For the time interval s, K (k) (τ) is the average value obtained;

[0221] express The derivative of express The derivative, Represents the first time interval in time interval s. Individual aircraft control quantities; Indicates the first control variable. Indicates the second control variable;

[0222] The correction parameter based on the s-th time interval is expressed as r. s The specific process for updating the time interval is as follows:

[0223] When r s ≤r max At that time, a correction strategy that increases the degree of the interpolation polynomial is used to update the degree of the polynomial and the number of collocations in the s-th time interval;

[0224] When r s ≥r max When the s-th time interval is divided into D time intervals, D is a positive integer greater than 2;

[0225] The other steps and parameters are the same as those in one of the specific implementation methods one to eight.

[0226] Specific Implementation Method Ten: The difference between this implementation method and Specific Implementation Methods One through Nine is that...

[0227] The correction strategy of increasing the interpolation polynomial degree is used to update the polynomial degree and collocation number in the s-th time interval; the specific process is as follows:

[0228] S5.5.1: Approximation error ε calculated based on the NLP problem results for the s-th time interval. (s) and setting error The degree of the interpolation polynomial for the s-th time interval is updated, expressed by the formula:

[0229]

[0230] In the formula, the ceil function is the integer function that rounds up the adjacent larger values; p represents the degree of the interpolation polynomial after the update in the s-th time interval. s Let ε represent the degree of the interpolation polynomial before the update of the s-th time interval, where A is a parameter greater than zero. (s) This represents the approximation error of the NLP problem calculation results for the s-th time interval; Indicates the error range;

[0231] The degree update of the interpolation polynomial can be adjusted by changing the value of parameter A;

[0232] S5.5.2: The approximation error ε is calculated based on the results of the NLP problem in the s-th time interval. (s) and setting error The number of coordinate points for the s-th time interval is updated, expressed by the formula:

[0233]

[0234] B is an intermediate parameter; the function density of the intermediate parameter B is expressed by the formula:

[0235] ρ B (τ)=CK(τ) 1 / 3 (17)

[0236] C is an intermediate parameter, and the function density of the intermediate parameter C is expressed by the formula:

[0237]

[0238] The probability distribution function of the intermediate parameter C is expressed by the formula:

[0239]

[0240] ζ represents a probability variable.

[0241] The updated collocation positions are obtained based on the probability distribution function of the intermediate parameter C, and can be expressed by the formula:

[0242]

[0243] v represents the updated allocation point.

[0244] The other steps and parameters are the same as those in any of the specific implementation methods one to nine.

[0245] The above description is merely of preferred embodiments of the present invention. It should be understood that the present invention is not limited to the specific embodiments described above. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent substitutions, and improvements made to the above embodiments without departing from the scope of the present invention, based on the technical essence of the present invention, and within the spirit and principles of the present invention, shall still fall within the protection scope of the present invention.

Claims

1. A method for optimizing the trajectory of a reusable spacecraft during the ascent phase with combined propulsion boost, characterized in that, Includes the following steps: Step S1: Construct a coordinate system for the reusable spacecraft; Step S2: Determine the spacecraft configuration of the reusable spacecraft, and establish a three-degree-of-freedom dynamic model of the spacecraft based on the coordinate system and spacecraft configuration of the reusable spacecraft; Step S3: Obtain existing spacecraft data, and determine the propulsion engine parameters and aerodynamic parameters of the reusable spacecraft based on the existing spacecraft data; Step S4: Based on the three-degree-of-freedom dynamic model of the aircraft, the parameters of the aircraft's power engine and aerodynamic parameters, and the coordinate system of the reusable spacecraft, construct the trajectory optimization problem for the climb phase of the reusable spacecraft; Step S5: The hp adaptive Gaussian pseudospectral method is used to optimize the climb trajectory of the reusable spacecraft, obtaining the optimal climb trajectory for the reusable spacecraft; the specific process is as follows: S5.1: Divide the entire time of the aircraft trajectory optimization problem into S... t A time interval; S5.2: Determine the initial number of collocations in each time interval; where the initial number of collocations in the s-th time interval is Ms; s∈S t S t Let S be a positive integer; for each S t The spacecraft trajectory optimization problem for each time interval is discretized to obtain S. t NLP problems within a time interval; S5.3: Using NLP algorithms to solve S t Solving NLP problems across time intervals yields S t Calculation results for NLP problems in each time interval; Calculate S t Approximation error of NLP problem calculation results for each time interval; S5.4: Set target error Calculate the approximation error of the NLP problem calculation results for S time intervals, and determine whether the approximation error of the NLP problem calculation results for all S time intervals is less than the target error. If the approximation error of the NLP problem calculation result in the s-th time interval is ε (s) Greater than the target error Then proceed to S5.5; If the approximation error of the NLP problem calculation results in the St time intervals is less than the target error Then stop iterating. Then, the NLP problem calculation results of the St time intervals are used to form the optimal climbing trajectory for the entire time interval; S5.5: Set the maximum correction parameter r max Calculate the correction parameters for each time interval based on the NLP problem calculation results for each time interval; update each time interval according to the correction parameters for each time interval, and return to S5.2 after updating; Specifically, a correction strategy that increases the degree of the interpolating polynomial is used to update the polynomial degree and the number of collocation points in the s-th time interval; the specific process is as follows: S5.5.1: Approximation error ε calculated based on the NLP problem results for the s-th time interval. (s) and setting error The degree of the interpolation polynomial for the s-th time interval is updated, expressed by the formula: In the formula, the ceil function is the integer function that rounds up the adjacent larger values; p represents the degree of the interpolation polynomial after the update in the s-th time interval. s Let ε represent the degree of the interpolation polynomial before the update of the s-th time interval, where A is a parameter greater than zero. (s) This represents the approximation error of the NLP problem calculation results for the s-th time interval; Indicates the error range; S5.5.2: The approximation error ε is calculated based on the results of the NLP problem in the s-th time interval. (s) and target error The number of coordinate points for the s-th time interval is updated, expressed by the formula: B is an intermediate parameter; the function density of the intermediate parameter B is expressed by the formula: r B (τ)=CK(τ) 1 / 3 (17) In the formula, τ represents the time variable, K() represents the curvature function, and C is an intermediate parameter. The function density of the intermediate parameter C is expressed by the formula: The probability distribution function of the intermediate parameter C is expressed by the formula: ζ represents a probability variable. The updated collocation positions are obtained based on the probability distribution function of the intermediate parameter C, and can be expressed by the formula: v represents the updated allocation point.

2. The method for optimizing the ascent trajectory of a reusable spacecraft with combined propulsion boost as described in claim 1, characterized in that, The coordinate system of the reusable spacecraft in S1 includes: 6 coordinate systems and 5 coordinate system transformation relationships; The six coordinate systems include: the geocentric inertial coordinate system o e -x I y I z I Navigation coordinate system o n -x n y n z n Geographic coordinate system T -x T y T z T , aircraft body coordinate system o b -x b y b z b Velocity coordinate system o v -x v y v z v Half-velocity coordinate system o h -x h y h z h ; The five coordinate system transformation relationships include: the transformation relationship between the navigation coordinate system and the aircraft body coordinate system. Transformation relationship between the aircraft body coordinate system and the velocity coordinate system Transformation relationship between half-velocity coordinate system and velocity coordinate system Transformation between geographic coordinate system and semi-velocity coordinate system Conversion between navigation coordinate system and geographic coordinate system The construction of the geocentric inertial coordinate system o e -x I y I z I The specific process is as follows: The origin of the geocentric inertial coordinate system is the center of the Earth, denoted as o. e The coordinate axes o of the geocentric inertial coordinate system e x I， coordinate axis o e y I， coordinate axis o e z I Construct a right-handed rectangular coordinate system; The coordinate axes o of the geocentric inertial coordinate system e x I In the equatorial plane, the vernal equinox is pointed to at a fixed time. The coordinate axes o of the geocentric inertial coordinate system e z I Perpendicular to the equatorial plane; pointing towards the North Pole. The construction of the navigation coordinate system o n -x n y n z n The specific process is as follows: Using the launch point of the aircraft as the navigation coordinate system n -x n y n z n The origin of the coordinate system is denoted as o. n , The coordinate axes o of the navigation coordinate system n z n , coordinate axis o n x n , coordinate axis o n y n Construct a right-handed rectangular coordinate system; The coordinate axes o of the navigation coordinate system n z n Within the line connecting the Earth's center and the spacecraft's launch point, pointing upwards, The coordinate axes o of the navigation coordinate system n x n Within the meridian plane where the aircraft is located, it points north; The construction of the geographic coordinate system T -x T y T z T The specific process is as follows: A geographic coordinate system is constructed using the intersection of the line connecting the Earth's center and the spacecraft's center of mass with the elliptical surface of the Earth. T -x T y T z T The origin o t , coordinate axes of a geographic coordinate system t z t , coordinate axis o t x t , coordinate axis o t y t Construct a right-handed rectangular coordinate system; coordinate axes of a geographic coordinate system t y t The line connecting the Earth's center and the spacecraft's center of mass coincides with the line pointing upwards. coordinate axes of a geographic coordinate system t x t Within the meridian plane where the aircraft is located, pointing north, The construction of the aircraft body coordinate system o b -x b y b z b The specific process is as follows: With the spacecraft's center of mass as the origin of the spacecraft's body coordinate system, denoted as o. b , The coordinate axes of the aircraft body coordinate system b z b , coordinate axis o b x b and coordinate axis o b y b Construct a right-handed rectangular coordinate system; The coordinate axes of the aircraft body coordinate system b x b Aligned with the longitudinal axis of the aircraft's fuselage, pointing towards the head, The coordinate axes of the aircraft body coordinate system b y b Located within the longitudinal symmetry plane of the aircraft's fuselage, pointing upwards; The constructed velocity coordinate system o v -x v y v z v The specific process is as follows: With the spacecraft's center of mass as the origin of the velocity coordinate system, denoted as o. v , The coordinate axes o of the velocity coordinate system v z v , coordinate axis o v x v and coordinate axis o v y v Construct a right-handed rectangular coordinate system; The coordinate axes o of the velocity coordinate system v x v Coinciding with the direction of the aircraft's velocity vector; The coordinate axes o of the velocity coordinate system v y v Located within the longitudinal symmetry plane of the aircraft, pointing upwards; The construction of the half-velocity coordinate system o h -x h y h z h The specific process is as follows: With the spacecraft's center of mass as the origin of the half-velocity coordinate system, denoted as o h , The coordinate axis o of the half-velocity coordinate system h z h , coordinate axis o h x h and coordinate axis o h y h Construct a right-handed rectangular coordinate system; The coordinate axis o of the half-velocity coordinate system h x h The axis coincides with the direction of the aircraft's velocity vector; The coordinate axis o of the half-velocity coordinate system h y h Located in the vertical plane containing the velocity vector, pointing upwards; The transformation relationship between the navigation coordinate system and the aircraft body coordinate system Expressed as a formula: Define pitch angle Yaw angle ψ, roll angle γ, establish the transformation relationship between the navigation coordinate system and the aircraft body coordinate system: The transformation relationship between the aircraft body coordinate system and the velocity coordinate system is constructed. The specific process is as follows: Define the angle of attack α and the sideslip angle β, and establish the transformation relationship between the aircraft body coordinate system and the velocity coordinate system as follows: The transformation relationship between the constructed half-velocity coordinate system and the velocity coordinate system is described. The specific process is as follows: Define the velocity tilt angle γ ν The transformation relationship between the half-velocity coordinate system and the velocity coordinate system is established as follows: The transformation relationship between the geographic coordinate system and the semi-velocity coordinate system is constructed. The specific process is as follows: Define the velocity deflection angle σ and the velocity tilt angle θ, then the transformation relationship between the geographic coordinate system and the semi-velocity coordinate system is as follows: The relationship between the navigation coordinate system and the geographic coordinate system is constructed. The specific process is as follows: Define the geocentric latitude difference between the navigation coordinate system and the geographic coordinate system as Δφ1, and the longitude difference as Δλ1. Then the transformation relationship between the two is:

3. The method for optimizing the ascent trajectory of a reusable spacecraft with combined propulsion boost according to claim 2, characterized in that, In step S2, the spacecraft configuration of the reusable spacecraft is determined, and a three-degree-of-freedom dynamic model of the spacecraft is established based on the coordinate system and configuration of the reusable spacecraft; the specific process is as follows: S2.1: Determine the aircraft configuration as Lockheed Martin SR-72 and construct a dynamic model of the aircraft in a geocentric inertial coordinate system; S2.2: Construct the dynamic model of the aircraft in the navigation coordinate system based on the dynamic model of the aircraft in the geocentric inertial frame: S2.3: Based on the transformation relationship between the navigation coordinate system and the geographic coordinate system, and the transformation relationship between the geographic coordinate system and the semi-velocity coordinate system, The dynamic model of the aircraft in the navigation coordinate system is decomposed in the half-velocity system to obtain the three-degree-of-freedom dynamic model of the aircraft. The three-degree-of-freedom dynamic model of the aircraft includes: the dynamic equation of the aircraft's center of mass, the kinematic equation of the aircraft's center of mass, and the equation of change of the aircraft's mass.

4. The method for optimizing the ascent trajectory of a reusable spacecraft with combined propulsion boost according to claim 3, characterized in that, The specific process of constructing the dynamic model of the aircraft in the geocentric inertial coordinate system in S2.1 is as follows: In the formula, m is the mass of the aircraft, t represents time, and r e denoted as the geocentric distance of the spacecraft's location in the geocentric inertial coordinate system, P is the engine thrust of the spacecraft, R is the aerodynamic force, and g is the gravitational acceleration. In step S2.2, the dynamic model of the aircraft in the navigation coordinate system is constructed based on the dynamic model of the aircraft in the geocentric inertial frame. The specific process is as follows: F e =-mω e ×(ω e ×(R0+r m )) In the formula, ω e The angular velocity of the navigation coordinate system relative to the inertial coordinate system; r e =r m +R0, r m R0 represents the position vector of the aircraft in the navigation system, and F represents the vector from the Earth's center to the origin of the navigation coordinate system. c F represents the Coriolis inertial force. e Indicates the inertial force involved; In step S2.3, based on the transformation relationships between the navigation coordinate system and the geographic coordinate system, and between the geographic coordinate system and the half-velocity coordinate system, the dynamic model of the aircraft in the navigation coordinate system is decomposed in the half-velocity coordinate system to obtain a three-degree-of-freedom dynamic model of the aircraft. The specific process is as follows: S2.3.1: Decompose the dynamic model of the aircraft in the navigation coordinate system into a half-velocity system to obtain the dynamic equation of the aircraft's center of mass, which can be expressed as: In equation (9), V represents the speed of the aircraft, γ v б represents the aircraft's velocity tilt angle, θ represents the aircraft's velocity deflection angle, and θ represents the aircraft's velocity inclination angle. The derivative of the aircraft's velocity V. Indicates the aircraft's velocity tilt angle γ v The derivative, The derivative representing the velocity deflection angle б, P xh This indicates the thrust along the coordinate axis o of the half-velocity coordinate system. h x h The component, P yh This indicates the thrust along the coordinate axis o of the half-velocity coordinate system. h y h The component, P zh This indicates the thrust along the coordinate axis o of the half-velocity coordinate system. h z h The amount, ω θ The rotational angular velocity is represented along the coordinate axis o of the half-velocity coordinate system. h y h The component, ω б The rotational angular velocity is represented along the coordinate axis o of the half-velocity coordinate system. h x h The component, ω V The rotational angular velocity is represented along the coordinate axis o of the half-velocity coordinate system. h x h The amount, m represents the mass of the aircraft; X represents the drag force of the aircraft; Y represents the lift force of the aircraft; Z represents the lateral force of the aircraft. g θ This represents the Earth's gravitational force along the coordinate axis o of the half-velocity coordinate system. h x h The amount, g V This represents the Earth's gravitational force along the coordinate axis o of the half-velocity coordinate system. h y h The amount, g б This represents the Earth's gravitational force along the coordinate axis o of the half-velocity coordinate system. h z h The components, where g r G represents the component of gravitational acceleration along the Earth's center. ω α represents the component of gravitational acceleration along the direction of Earth's rotation; r represents the Earth's radius; α e The flatness of the Earth is represented by φ, latitude by J, and zone harmonic coefficient by μ. E Represents the gravitational constant. S2.3.2: Based on the dynamic equations of the aircraft's center of mass, the kinematic equations of the aircraft's center of mass are constructed as follows: In the formula, Represents the rate of change of the geocentric radius vector. Indicates the rate of change of latitude. Indicates the rate of change of longitude; S2.3.3: The equation for the change in the mass of the aircraft is constructed as follows: In the formula, f represents the rate of change of the aircraft's mass. m The function represents the change in mass, where Ma represents the Mach number, h represents the altitude, and c represents the engine operating coefficient.

5. The method for optimizing the ascent trajectory of a reusable spacecraft with combined propulsion boost according to claim 4, characterized in that, In step S3, existing aircraft data is acquired, and the aircraft's power engine parameters and aerodynamic parameters are determined based on this data. The specific process is as follows: S3.1: Selecting a parallel turbojet / ramjet combined-engine as the aircraft's propulsion system, constructing fitting formulas for the turbojet engine parameters and ramjet engine parameters, expressed as: In the formula, P1 represents the thrust of the turbojet engine, P2 represents the thrust of the ramjet engine, dm1 represents the fuel consumption rate per second of the turbojet engine, dm2 represents the fuel consumption rate per second of the ramjet engine, and f p1 f represents the thrust function of a turbojet engine. p2 f represents the thrust function of a ramjet engine. m1 f represents the fuel consumption rate per second function of a turbojet engine. m2 This represents the fuel consumption rate per second function of a ramjet engine. S3.2: Obtain existing flight data of turbojet engine aircraft and ramjet engine aircraft. Based on the existing flight data of turbojet engine aircraft and ramjet engine aircraft, obtain the engine parameters of the aircraft through interpolation fitting. The engine parameters of the aircraft include: turbojet engine parameters and ramjet engine parameters; the specific process is as follows: Based on existing flight data of turbojet engine aircraft and fitting formulas for turbojet engine parameters, the parameters of turbojet-powered engines are determined by interpolation fitting. Based on existing data on the flight of ramjet-powered aircraft and the fitting formula for ramjet engine parameters, the parameters of the ramjet-powered engine are determined by interpolation fitting. S3.3: Obtain the fixed parameter data of the aircraft; based on the existing fixed parameter data of the aircraft, use an interpolation fitting method to obtain the aerodynamic parameters; the aerodynamic parameters include: aerodynamic force, aerodynamic coefficient, and lift-to-drag ratio coefficient; the specific process is as follows: S3.3.1: Construct fitting formulas for aerodynamic forces and aerodynamic coefficients. The specific process is as follows: In the formula, S represents the area of ​​force application, q represents the dynamic pressure, α represents the angle of attack, β represents the sideslip angle, and f x The drag coefficient function, f y Represents the lift coefficient function, f z Represents the lateral force coefficient function; The aerodynamic forces include: drag X, lift Y, and lateral force Z; The aerodynamic coefficient includes the drag coefficient C. X Lift coefficient C Y Lateral force coefficient C Z ; S3.3.2: Determine the fixed parameters of the aircraft based on the aircraft configuration and obtain the fixed parameter data of the aircraft; based on the fixed parameter data of the aircraft, the parameters of the aircraft's power engine and the calculation formula of the aerodynamic coefficient, fit the aerodynamic coefficient through interpolation fitting method; S3.3.3: Based on the aerodynamic coefficient and its calculation formula, the aerodynamic force is obtained by interpolation fitting, and then the lift-to-drag ratio coefficient is obtained based on the aerodynamic force.

6. The method for optimizing the ascent trajectory of a reusable spacecraft with combined propulsion boost according to claim 5, characterized in that, In step S4, based on the three-degree-of-freedom dynamics model of the spacecraft, the parameters of the spacecraft's propulsion engine and aerodynamic parameters, and the coordinate system of the reusable spacecraft, a trajectory optimization problem for the climb phase of the reusable spacecraft is constructed; the specific process is as follows: S4.1: Select the control variables for the trajectory optimization problem of the climb phase of a reusable spacecraft; based on the transformation relationship between the control variables and the spacecraft's body coordinate system and velocity coordinate system, obtain the transformed three-degree-of-freedom dynamic model of the spacecraft; S4.2: Construct a state constraint model for the flight's climb phase based on the transformed three-degree-of-freedom dynamic model of the aircraft, the parameter model of the aircraft's power engine, and the aerodynamic parameters of the aircraft; S4.3: Establish performance optimization indices for the trajectory of the spacecraft's climb phase; based on the performance optimization indices, control variables, and state constraint model of the spacecraft's climb phase, construct the trajectory optimization problem for the climb phase of a reusable spacecraft.

7. The method for optimizing the ascent trajectory of a reusable spacecraft with combined propulsion boost according to claim 6, characterized in that, The control variables in the climb trajectory optimization problem of the reusable spacecraft in S4.1 include: a first control variable and a second control variable; wherein, the first derivative of the angle of attack is selected as the first control variable, and the first derivative of the velocity tilt angle is selected as the second control variable. Based on the transformation relationship between the state variables and the body coordinate system and velocity coordinate system of the aircraft, the transformed three-degree-of-freedom dynamic model of the aircraft is obtained; expressed by the formula: In the formula, Let α be the first derivative of the angle of attack. For the velocity tilt angle γ v The first derivative, and To control the quantity; and For state variables; In S4.2, a state constraint model for the flight climb phase is constructed based on the transformed three-degree-of-freedom dynamic model of the aircraft, the parameters of the aircraft's power engine, and the aerodynamic parameters of the aircraft. The state constraint model for the aircraft's climb phase includes: a turbojet climb phase constraint model, a combined climb phase constraint model, a ramjet climb phase constraint model, and a full-range climb phase constraint model; the specific process is as follows: S4.2.1: Construct a constraint model for the turbojet's climb section. The constraint model for the turbojet climbing section includes: initial state constraints for the turbojet climbing section and terminal state constraints for the turbojet climbing section. The initial state constraints of the turbojet climb section include: The initial height of the turbojet's climb phase is hw0 = 10 km, and the initial velocity of the turbojet's climb phase is Vw0 = 408 m / s. 2 The initial velocity inclination angle of the turbojet climb section is θw0 = 0°, the initial velocity deflection angle of the turbojet climb section is σw0 = -21.684°, and the initial mass of the turbojet climb section is mw0 = 12000 kg. The terminal state constraints of the turbojet climb section include: Terminal velocity Vw of the turbojet climb phase f =544m s, the terminal velocity inclination angle of the turbojet climb section θwf=0°, and the ground center distance increment of the terminal turbojet climb section Δwr1=4500m; S4.1.2: Construct a combined climbing section constraint model; The combined climbing segment constraint model includes: initial state constraints of the combined power climbing segment and terminal state constraints of the combined power climbing segment; The initial state constraints of the combined dynamic climbing section include: Initial velocity V of the combined power climb phase z θ = 544 m / s, initial velocity inclination angle θz0 = 0°, initial geocentric distance r during combined dynamic climb. z =rz0 + Δrz1; The combined power climbing section's terminal state constraints include: Combined power climbing phase terminal speed V z f = 850 m / s, and the ground center distance increment at the end of the combined power climbing section is Δzr = 2500 m; S4.1.3: Construct a constraint model for the stamping climbing section; The constraint model of the stamping climbing section includes: initial state constraints of the stamping climbing section and terminal state constraints of the stamping climbing section; The initial state constraints of the stamping climbing section include: Initial velocity of the ram-climb section Vc0 = 850 m / s; initial distance from the ground rc of the ram-climb section rc = rc0 + Δrc 1+ Δrc2; The terminal state constraints of the stamping climbing section include: The terminal speed Vc of the stamping climbing section f =1530m / s, the increase in distance from the center of the earth at the end of the ram-press climbing section Δcr3 = 11500m, and the velocity inclination angle at the end of the ram-press climbing section θc f =0°, angle of attack α at the end of the stamping climb section cf =4.12°; S4.1.4: Construct a constraint model for the entire climbing section; The constraint model for the entire climbing section includes: Maximum overload constraint n throughout the climbing section max =3 Maximum dynamic pressure constraint q throughout the entire climbing section max =100 kPa, maximum heat flux density constraint throughout the climb section The climbing phase has the following constraints: speed tilt angle constraint -0.5° < θ < 10°, angle of attack constraint 1° < α < 10°, speed tilt angle constraint υ = 0°, and control quantity constraints |u1| ≤ 0.1° / s, |u2| ≤ 0° / s; total thrust P and total fuel consumption rate dm during the combined power climbing phase. The total thrust P and total fuel consumption rate dm of the combined power climbing phase are expressed by the following formulas: P=c1(t)P1(Ma,H)+c2(t)P2(Ma,H) dm=c1(t)dm1(Ma,H)+c2(t)dm2(Ma,H) (22) Where c1 and c2 are both functions of time t, P1 is the thrust of the ramjet engine, P2 is the thrust of the turbojet engine, dm1 is the fuel consumption rate per second of the ramjet engine, dm2 is the fuel consumption rate per second of the turbojet engine, c1=0.7+0.3(t-t1) / (t2-t1), c2=1-(t-t1) / (t2-t1); The performance optimization index for the aircraft's climb trajectory is established in S4.3; expressed by the formula: In the formula, R d Where C is the radius of the aircraft's nose, C1 is a constant, and ρ0 is the standard atmospheric density at sea level. g0 is the Earth's gravitational pull at sea level, and R0 is the Earth's average radius.

8. The method for optimizing the ascent trajectory of a reusable spacecraft with combined propulsion boost according to claim 7, characterized in that, In S5.2, the aircraft trajectory optimization problem for a given time interval is discretized to obtain an NLP problem for that time interval. The specific process is as follows: S5.2.1: Transform the time interval of the aircraft trajectory optimization problem to the time interval for applying the Gauss pseudospectral method; S5.2.2: Obtain K collocation points in the transformed time interval; use the number of collocation points as the degree of the Largerange interpolation polynomial, where K is a positive integer; Based on K collocation points, construct a Larger interpolation polynomial, and use the Larger interpolation polynomial to discretize and represent the state and control variables in the time interval of the aircraft trajectory optimization problem; S5.2.3: Convert the discretized representation of the state and control variables of the aircraft trajectory optimization problem in time interval into algebraic form to obtain the algebraic equations satisfied by the state variables on the collocation points; construct an NLP problem in time interval based on the algebraic equations satisfied by the state variables on the collocation points; In S5.2.1, the time interval of the aircraft trajectory optimization problem is transformed into the time interval for applying the Gauss pseudospectral method; The specific process is as follows: The time interval for the aircraft trajectory optimization problem is [t0,t]. f The time interval [-11] for applying the Gauss pseudospectral method is converted as follows: t represents the time interval [t0, t] f In the time variable ], t0 represents the initial time, t f Indicates the terminal time; This represents the time variable t transformed into the time variable of the Gauss pseudospectral method; In S5.2.2, the state and control variables of the aircraft trajectory optimization problem are discretized using Larger interpolation polynomials, expressed by the following formula: in, Represents the state variables of the aircraft. Represents the discretized representation of the aircraft state variables. For matching points The state variables of the aircraft For matching points Discretized representation of aircraft state variables Let K be the i-th collocation point, K be the total number of collocation points, and i be the interpolation sequence number. For aircraft control variables, For discretized representation of aircraft control quantities; Represents the Lagrange basis functions; Indicates matching points Aircraft control quantities; Indicates matching points Discretized representation of aircraft control variables; X represents the aircraft state at the terminal moment; f Let f() represent the discretized terminal state variables of the aircraft, and let f() represent the transformed three-degree-of-freedom dynamic model of the aircraft. Indicates Gaussian weights. Locate the k-th Legendre-Gauss point; t f The terminal time is t0, which is the start time. In S5.2.3, the state and control variables of the discretized representation of the aircraft trajectory optimization problem over the time interval are converted into algebraic form to obtain the algebraic equations satisfied by the state variables at the collocation points. The specific process is as follows: Among them, D ki () represents the weighting coefficient. Let k represent the k-th collocation point, where k = 1, ..., K. express The derivative, Indicates matching points Discretized representation of aircraft control variables; express The derivative, Indicates matching points Aircraft control quantities; Indicates matching points The state variables of the aircraft; In step S5.5, the correction parameter r for the s-th time interval is calculated based on the calculation results of the NLP problem for the s-th time interval. s The specific process is as follows: K max (s) For the time interval s, K (s) The maximum value obtained by (τ) For the time interval s, K (k) (τ) is the average value obtained; express The derivative, express The derivative, Represents the first time interval in time interval s. Individual aircraft control quantities; Indicates the first control variable. Indicates the second control variable; The correction parameter based on the s-th time interval is expressed as r. s The specific process for updating the time interval is as follows: When r s ≤r max At that time, a correction strategy that increases the degree of the interpolation polynomial is used to update the degree of the polynomial and the number of collocations in the s-th time interval; When r s ≥r max When the s-th time interval is divided into D time intervals, D is a positive integer greater than 2.

Images