A dimension reduction sequence convex optimization trajectory planning method for multi-revolution aerodynamic auxiliary aerospace vehicle skew orbit transfer

Abstract

The application discloses a dimension-reduction sequence convex optimization trajectory planning method for multi-orbit aerodynamic auxiliary skew orbit transfer of aerospace vehicles, and belongs to the technical field of numerical solution of optimal control problems. First, for each atmospheric flight segment, an improved Chebyshev-Picard iteration fixed point equation is established at discrete nodes, and the node state is implicitly determined according to the atmospheric inlet state, the node control quantity and the atmospheric flight time. Second, a convex optimization sub-problem is constructed, and the Jacobian matrix is calculated by using the improved Chebyshev-Picard iteration fixed point equation and the discrete adjoint equation. Third, the trial solution is verified by using the original nonlinear dynamics, the trust region is updated, and the convergence criterion of the iteration of the convex optimization sub-problem is given. Finally, the discrete grid order is adjusted according to the discrete thrust switch function, and the judgment criterion for the convergence of the grid improvement iteration is given. The method can reduce the dimension of the optimization problem, improve the solution efficiency of the multi-orbit fuel optimal trajectory planning, and the trajectory has good feasibility.

Description

Technical Field

[0001] This invention belongs to the field of numerical solution technology for optimal control problems, specifically relating to a dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit changes in aerospace vehicles. The method can be used for overall aerospace vehicle design, orbital maneuver mission planning, aerodynamic-assisted orbital transfer performance evaluation, and nominal trajectory generation for flight control. Background Technology

[0002] Aerodynamic-assisted orbital maneuvering is a type of orbital maneuvering that utilizes the aerodynamic forces generated after an aircraft enters the atmosphere to change the direction of its orbital angular momentum, and combines this with the propulsion system to compensate for the orbital energy loss caused by aerodynamic drag, thereby completing changes in orbital inclination and semi-major axis.

[0003] Limited by minimum flight altitude, maximum heat flux, and maximum thrust, the orbital inclination changes achievable in a single atmospheric pass are typically limited. To accomplish larger non-plane orbital maneuvers, multi-cycle atmospheric-assisted orbital maneuvers are often employed: the spacecraft enters the atmosphere from its initial orbit with a descent velocity pulse; during each atmospheric pass, it is controlled collaboratively by angle of attack, bank angle, and thrust; after exiting the atmosphere, it flies along an exoatmospheric orbit to its apogee; inter-pass adjustment velocity pulses are applied to adjust the atmospheric inlet position for the next pass; and after the final atmospheric pass, terminal injection or circular maneuvers are performed at the apogee to enter the target orbit. The atmospheric flight segments, inter-pass exoatmospheric orbits, initial descent, and terminal injection are all coupled, forming a multi-stage, multi-constraint, and highly nonlinear optimal control problem.

[0004] Existing trajectory planning methods often employ direct collocation methods, including collocation methods based on difference schemes and pseudospectral methods, which transform continuous optimal control problems into nonlinear programming problems for solution. For multi-cycle aerodynamically assisted non-plane trajectory change problems, as the number of atmospheric flight cycles and discrete nodes increases, the scale of nodal state variables and dynamic equality constraints grows rapidly, resulting in large nonlinear programming problem scales and long computation times. After the discrete solution is restored to continuous control, multi-cycle error accumulation occurs.

[0005] Sequential convex optimization methods can approximate nonlinear programming problems as convex optimization problems near a reference trajectory through linearization, and then solve them using mature convex optimization solvers. However, even when retaining all nodal state variables, the convex optimization subproblem still has a large scale. If numerical differencing is used to construct the Jacobian matrix, multiple repeated integrations are required, and the computation time increases with the number of optimization variables. Furthermore, the Jacobian matrix constructed through numerical differencing is affected by the difference step size. For aircraft control in the context of fuel optimization, its thrust typically exhibits a near-pop-pop control structure. If the discretization accuracy of the grid nodes is not high, even after solving the optimization problem on a fixed grid to achieve a small residual, the accuracy may deteriorate after reverting to continuous control.

[0006] The improved Chebyshev-Pika iteration can establish a fixed-point equation based on the Chebyshev polynomial integral matrix, relating the atmospheric flight state to the atmospheric inlet state, nodal control variables, and flight time. This equation can reduce the dimensionality of optimization variables, and further combining it with sequential convex optimization can eliminate dynamic equality constraints. In recent years, some studies have applied the improved Chebyshev-Pika iteration to single-stage optimal control problems. However, to apply it to multi-stage optimal control problems such as multi-cycle aerodynamically assisted non-planar trajectory changes, the following issues need to be addressed: the atmospheric flight segments in each cycle are not simply listed and repeatedly controlled, but rather coupled and influenced by each other. How to quantitatively express the influence of subsequent cycles on previous cycles, i.e., how to calculate the system's Jacobian matrix, is crucial. If the grid error accuracy is low, it will affect the error of the discrete solution, causing the obtained discrete solution to lose physical feasibility under the original nonlinear dynamics due to error accumulation.

[0007] Therefore, in the optimization of multi-cycle aerodynamic-assisted non-plane trajectory changes, there is an urgent need for a trajectory planning method that can reduce the solution scale while taking into account iterative stability, constraint handling capability, and continuous trajectory feasibility, so as to improve the numerical solution efficiency and engineering applicability of such problems. Summary of the Invention

[0008] To address the challenges of numerous nodal state variables, large-scale dynamic equation constraints, strong multi-stage coupling effects, and multi-cycle error accumulation in multi-cycle aerodynamically assisted non-planar trajectory planning, this invention provides a dimensionality-reduced sequential convex optimization trajectory planning method for multi-cycle aerodynamically assisted non-planar trajectory planning of aerospace vehicles. This method is primarily used to calculate offline a physically feasible trajectory that saves fuel compared to pure propulsion-based trajectory changes, given initial orbit, target orbit, spacecraft model, and mission constraints such as altitude, heat flux, and control boundaries. It obtains the angle of attack, bank angle, thrust, flight time, initial descent velocity pulse, inter-cycle adjustment velocity pulse, and final orbit insertion velocity pulse during atmospheric flight. This method reduces the dimensionality of the optimization problem and improves the solution efficiency of multi-cycle fuel-optimal trajectory planning.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] In a first aspect, the present invention provides a dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit changes of aerospace vehicles, comprising:

[0011] Step 1: Obtain relevant parameters for the multi-cycle aerodynamic-assisted orbital maneuvering mission. Based on these parameters, establish a mathematical model for the optimal control problem of the multi-cycle aerodynamic-assisted orbital maneuvering mission, which includes the initial descent, atmospheric flight, extra-atmospheric flight, and final orbit insertion processes. The relevant parameters include initial orbital parameters, target orbital parameters, aerospace vehicle parameters, mission constraint parameters, and the number of atmospheric flight cycles.

[0012] Step 2: Generate the initial reference trajectory for each orbit based on the preset initial descent variables, preset inter-orbit adjustment variables, preset constant angle of attack, and preset constant tilt angle for each orbit; the initial reference trajectory for each orbit includes the initial reference values ​​of atmospheric inlet variables, node status, node control variables, and atmospheric flight time for each orbit.

[0013] Step 3: In the atmospheric flight segment of each orbit, establish an improved Chebyshev-Pika iterative fixed-point equation at discrete nodes, and represent the node state as an implicit function of atmospheric inlet variables, node control variables, and atmospheric flight time. Construct dimension reduction optimization variables based on the node control variables, atmospheric flight time, initial orbit reduction variables, inter-orbit adjustment variables, and atmospheric inlet variables of each orbit.

[0014] Step 4: Calculate the objective function, equality constraint residual, and inequality constraint residual at the current reference trajectory. Make a first-order approximation of the objective function, equality constraint residual, and inequality constraint residual with respect to the dimension reduction optimization variables. Introduce non-negative equality relaxation variables, non-negative inequality relaxation variables, and trust regions to construct a convex optimization subproblem.

[0015] Step 5: Calculate the first-order variation of the improved Chebyshev-Pika iterative fixed-point equation to obtain the relationship between the node state perturbation and the dimension reduction optimization variable perturbation; construct the discrete adjoint equation, and without explicitly calculating the matrix inverse, analytically obtain the Jacobian matrix of the objective function, the equality constraint residual, and the inequality constraint residual with respect to the dimension reduction optimization variable.

[0016] Step 6: Solve the convex optimization subproblem to obtain the trial step, sum the trial step with the current reference trajectory to form the trial solution; interpolate the discrete control quantity of the trial solution, perform numerical integration under the original nonlinear dynamics to obtain the aircraft trajectory corresponding to the trial solution, and then calculate the physical feasibility index, objective function, constraint violation degree and value function corresponding to the trajectory. Then, determine whether to accept the trial step based on whether the trajectory is within the set physical feasibility region, and calculate the ratio between the decrease in the value function predicted by the linearized model and the decrease in the actual value function of the original nonlinear dynamics. Update the size of the trust region based on the ratio.

[0017] Step 7: Based on the physical feasibility of the trajectory corresponding to the trial solution, and the changes in node control variables, atmospheric flight time, initial descent variables, inter-circle adjustment variables, and objective function in the solutions of the two adjacent convex optimization subproblems, determine whether the solution of the convex optimization subproblem has converged.

[0018] Step 8: When the solution of the convex optimization subproblem under the fixed grid order satisfies the convergence condition, construct the discrete thrust switching function for each discrete node based on the dual variables of the convex optimization subproblem, the terminal augmented gradient of each cycle, and the improved Chebyshev-Pika iterative fixed-point equation; determine the candidate thrust switching position based on the node thrust jump position and the zero of the discrete thrust switching function, and update the discrete grid with the alignment error between the candidate thrust switching position and the candidate new grid node or the midpoint of the adjacent node interval as the index;

[0019] Step 9: Interpolate the converged solution of the convex optimization subproblem corresponding to the previous generation grid to the updated discrete grid as a warm start, and repeat steps 1 to 8; For the converged solution of the convex optimization subproblem obtained in each grid iteration, start from the initial trajectory and use serial integration to obtain the true trajectory. Based on the terminal constraints, path constraints obtained after serial integration, as well as the converged solutions and changes in the objective function of two adjacent grid iterations, determine whether the outer grid iteration has converged; Output the trajectory planning result when the outer grid iteration converges.

[0020] Secondly, the present invention provides an electronic device, comprising: one or more processors; and a memory for storing one or more programs; wherein, when the one or more programs are executed by the one or more processors, the one or more processors implement the above-described method for dimensionality-reduced sequential convex optimization trajectory planning for multi-loop aerodynamic-assisted heterogeneous orbit changes of aerospace vehicles.

[0021] Thirdly, the present invention provides a computer-readable storage medium storing executable instructions thereon, which, when executed by a processor, enable the processor to implement the above-mentioned method for dimensionality-reduced sequential convex optimization trajectory planning for multi-loop aerodynamic-assisted heterogeneous orbit changes of aerospace vehicles.

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

[0023] 1. This invention employs a multi-targeting structure divided according to atmospheric flight cycles, preserving the atmospheric inlet state of each cycle. It utilizes an improved Chebyshev-Pikachu iterative fixed-point equation to implicitly determine the node states of each cycle based on the inlet state, node control variables, and flight time. This reduces independent optimization variables while preserving the state coupling relationships between cycles, decreasing the solution scale of convex optimization subproblems, and improving overall solution efficiency.

[0024] 2. Based on the first-order variational form and discrete adjoint equation of the improved Chebyshev-Pika iterative fixed-point equation, the Jacobian matrix of the objective function, equality constraint function, and inequality constraint function with respect to the dimension-reduced optimization variable can be analytically constructed. This eliminates the dynamic equality constraint and avoids the need to calculate the Jacobian matrix using the finite difference method. It also eliminates the need to explicitly calculate the inverse matrix of the large-scale dense matrix, which helps to reduce the computational amount of the Jacobian matrix and improve the overall solution efficiency.

[0025] 3. Addressing the characteristic of thrust control approaching "bang-bang" control in the optimal fuel trajectory, this invention determines candidate thrust switching regions based on the nodal thrust jump positions and the zero points of the discrete thrust switching function, and adjusts the grid order according to the alignment between grid nodes and candidate switching positions. This method improves the accuracy of discrete grid representation of discontinuous thrust control modes. Attached Figure Description

[0026] Figure 1 This is a flowchart of the dimension-reduced sequence convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to the present invention.

[0027] Figure 2 This is a flowchart of a dimension-reduced sequence convex optimization trajectory planning method for multi-cycle aerodynamic-assisted heterogeneous orbit change of aerospace vehicles, according to an embodiment of the present invention.

[0028] Figure 3 This is a schematic diagram of a multi-turn aerodynamic-assisted non-planar orbit change mission scenario according to an embodiment of the present invention.

[0029] Figure 4 This is a schematic diagram of the output curve and discrete thrust switching function curve obtained by the method of this embodiment of the invention. Detailed Implementation

[0030] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0031] The term “exemplary” as used herein means “serving as an example, embodiment, or illustration.” Any embodiment illustrated herein as “exemplary” is not necessarily to be construed as superior to or better than other embodiments.

[0032] The terms “first”, “second”, etc., are used to distinguish similar objects, not to describe or indicate a specific order or sequence.

[0033] The term "comprising" or any other similar term is intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus / device that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent in such process, method, article, or apparatus / device.

[0034] Furthermore, to better illustrate the present invention, numerous specific details are provided in the following detailed embodiments. Those skilled in the art should understand that the present invention can be practiced without certain specific details. In some instances, methods, means, elements, and circuits well known to those skilled in the art have not been described in detail in order to highlight the spirit of the invention.

[0035] Figure 1 The flowchart of the dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit changes of aerospace vehicles according to the present invention is as follows: Figure 1 As shown, the method includes:

[0036] Step 1: Obtain relevant parameters for the multi-cycle aerodynamic-assisted non-plane orbit change mission, and establish a mathematical model for the optimal control problem of the multi-cycle aerodynamic-assisted non-plane orbit change mission based on the relevant parameters. The model includes the initial orbit descent, atmospheric flight, extra-atmospheric flight, and final orbit insertion process. The relevant parameters for the multi-cycle aerodynamic-assisted non-plane orbit change mission include initial orbit parameters, target orbit parameters, aerospace vehicle parameters, mission constraint parameters, and the number of atmospheric flight cycles.

[0037] Step 2: Generate the initial reference trajectory for each orbit based on the preset initial descent variables, preset inter-orbit adjustment variables, preset constant angle of attack, and preset constant tilt angle. The initial reference trajectory for each orbit includes the initial reference values ​​of atmospheric inlet variables, node status, node control variables, and atmospheric flight time for each orbit.

[0038] For example, the process of generating the initial reference trajectory for each lap includes:

[0039] (1) Apply a preset initial descent variable on the initial orbit and integrate according to the two-body orbit dynamics until the first atmospheric boundary is reached;

[0040] (2) In atmospheric flight, a preset constant angle of attack and a preset constant side angle are used, and the thrust component is set in the form of compensating for the energy consumed by drag. The thrust is integrated according to the dynamics of the atmospheric section until the atmospheric exit is reached.

[0041] (3) Starting from the atmospheric exit, integrate according to the two-body orbital dynamics until the orbital apogee is reached. Apply the preset inter-cycle adjustment variable and integrate according to the two-body orbital dynamics until the next atmospheric entry is reached. Repeat atmospheric flight and extra-atmospheric flight until the last atmospheric flight ends.

[0042] (4) Interpolate the integrated state and control variables to the current discrete node to obtain the initial reference values ​​of atmospheric inlet variables, node states, node control variables and atmospheric flight time for each loop.

[0043] Step 3: In the atmospheric flight segment of each orbit, establish an improved Chebyshev-Pika iterative fixed-point equation at the discrete nodes, and represent the node state as an implicit function of atmospheric inlet variables, node control variables, and atmospheric flight time. Construct dimension reduction optimization variables based on the node control variables, atmospheric flight time, initial orbit reduction variables, inter-orbit adjustment variables, and atmospheric inlet variables of each orbit.

[0044] In the atmospheric flight segment of each orbit, discretization is performed using the Modified Chebyshev-Picard Iteration (MCPI). An implicit state propagation mapping of the nodal states with respect to atmospheric inlet variables, nodal control variables, and atmospheric flight time is established based on the Modified Chebyshev-Picard Iteration fixed-point equations. This represents the nodal states as implicit functions of atmospheric inlet variables, nodal control variables, and atmospheric flight time, no longer as independent optimization variables. Instead, a dimensionality-reduced optimization variable is formed, consisting of nodal control variables, atmospheric flight time, initial descent variables, inter-orbit adjustment variables, and atmospheric inlet variables for each orbit. The dimensionality-reduced optimization variables include the nodal control vectors, atmospheric flight time variables, initial descent velocity pulses, adjustment velocity pulses, and atmospheric inlet variables for each orbit.

[0045] Step 4: Calculate the objective function, equality constraint residual, and inequality constraint residual at the current reference trajectory. Perform a first-order approximation of the objective function, equality constraint residual, and inequality constraint residual with respect to the dimensionality reduction optimization variables. Introduce non-negative equality relaxation variables, non-negative inequality relaxation variables, and a trust region to construct a convex optimization subproblem.

[0046] For example, the objective function could be to maximize the mass of the spacecraft after terminal orbit insertion or circularization, or equivalently to minimize its negative.

[0047] Step 5: Calculate the first-order variation of the improved Chebyshev-Pika iterative fixed-point equation to obtain the relationship between the node state perturbation and the dimension reduction optimization variable perturbation; construct the discrete adjoint equation, and without explicitly calculating the matrix inverse, analytically obtain the Jacobian matrix of the objective function, the equality constraint residual, and the inequality constraint residual with respect to the dimension reduction optimization variable.

[0048] Step 6: Solve the convex optimization subproblem to obtain the trial step, sum the trial step with the current reference trajectory to form the trial solution; interpolate the discrete control quantity of the trial solution, perform numerical integration under the original nonlinear dynamics to obtain the aircraft trajectory corresponding to the trial solution, and then calculate the physical feasibility index, objective function, constraint violation degree and value function corresponding to the trajectory. Then, determine whether to accept the trial step based on whether the trajectory is within the set physical feasibility region, and calculate the ratio between the decrease in the value function predicted by the linearized model and the decrease in the actual value function of the original nonlinear dynamics. Update the size of the trust region based on the ratio.

[0049] Physical feasibility indicators may include perigee altitude, apogee altitude corresponding to atmospheric exit conditions, and maximum heat flux. When the actual trajectory corresponding to the trial solution exceeds the preset physical feasibility region, the trial step is rejected and the trust region is narrowed.

[0050] Step 7: Based on the physical feasibility of the trajectory corresponding to the trial solution, and the changes in node control variables, atmospheric flight time, initial descent variables, inter-circle adjustment variables, and objective function in the solutions of the two adjacent convex optimization subproblems, determine whether the solution to the convex optimization subproblem has converged.

[0051] Step 8: When the solution of the convex optimization subproblem under the fixed grid order satisfies the convergence condition, construct the discrete thrust switching function for each discrete node based on the dual variables of the convex optimization subproblem, the terminal augmented gradient of each cycle, and the improved Chebyshev-Pika iterative fixed-point equation; determine the candidate thrust switching position based on the node thrust jump position and the zero of the discrete thrust switching function, and update the discrete grid with the alignment error between the candidate thrust switching position and the candidate new grid node or the midpoint of the adjacent node interval as the index.

[0052] Step 9: Interpolate the converged solution of the convex optimization subproblem corresponding to the previous generation grid to the updated discrete grid as a warm start, and repeat steps 1 to 8; For the converged solution of the convex optimization subproblem obtained in each grid iteration, start from the initial trajectory and use serial integration to obtain the true trajectory. Based on the terminal constraints, path constraints obtained after serial integration, as well as the converged solutions and changes in the objective function of two adjacent grid iterations, determine whether the outer grid iteration has converged; Output the trajectory planning result when the outer grid iteration converges.

[0053] The serial integration starts from the initial orbit, and after the initial descent process, the atmospheric inlet state of the first orbit is obtained. Subsequently, the atmospheric flight trajectory and the extra-atmospheric flight trajectory of each orbit are obtained by sequential integration. At the inlet of each orbit, the state variables are not reset to the optimized atmospheric inlet variables.

[0054] The method of this invention establishes an improved Chebyshev-Pika iterative fixed-point equation at discrete nodes for each atmospheric flight segment, implicitly determining the node states based on atmospheric inlet conditions, node control variables, and atmospheric flight time, thereby eliminating state optimization variables. Secondly, the objective function and constraints are linearly expanded near the reference trajectory to construct a convex optimization subproblem including constraint relaxation, corresponding penalty terms, and a trust region. The Jacobian matrix in the convex optimization subproblem is analytically calculated using the first-order variational equation and discrete adjoint equation of the improved Chebyshev-Pika iterative fixed-point equation. Then, the trial solution is verified using primal nonlinear dynamics, and the trust region is updated. A convergence criterion for solving the convex optimization subproblem iteratively is given. Finally, the discrete grid order is adjusted according to the discrete thrust switching function, and a criterion for judging the convergence of the improved grid iteration is given. The implementation of this method can reduce the dimensionality of the optimization problem, improve the solution efficiency of multi-cycle fuel optimal trajectory planning, and maintain good trajectory feasibility under primal nonlinear dynamics.

[0055] The following detailed description of the dimension-reduced sequential convex optimization trajectory planning method for multi-cycle aerodynamic-assisted heterogeneous orbit change of aerospace vehicles, with reference to the accompanying drawings and embodiments, is provided below.

[0056] Figure 2 This is a flowchart of a dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted non-plane orbit changes of aerospace vehicles, according to an embodiment of the present invention. Figure 2 As shown, the method includes:

[0057] S101. Mathematical modeling of the optimal control problem for multi-turn aerodynamic-assisted non-plane orbit change.

[0058] First, the multi-turn aerodynamic-assisted non-plane orbit change mission scenario of the present invention will be described. Figure 3 This is a schematic diagram of a multi-turn aerodynamic-assisted non-planar orbit change mission scenario according to an embodiment of the present invention, such as... Figure 3 As shown, the multi-cycle aerodynamic-assisted orbital maneuver of the spacecraft involves four scenarios:

[0059] (1) Initial descent: Starting from the initial orbit, apply an initial descent velocity pulse in the opposite direction of the velocity. This causes the orbital perigee to enter the atmosphere;

[0060] (2) Atmospheric segment flight: in the first During atmospheric flight, the aircraft uses angle of attack Tilt angle and thrust The coordinated control generates out-of-plane lift, changes the orbital inclination, and adjusts the orbital energy, enabling the spacecraft to fly out of the atmosphere again while changing its orbital inclination.

[0061] (3) Extra-atmospheric flight: in front After orbiting the atmosphere, the spacecraft exits the atmosphere and reaches its apogee along a two-body orbit. At the apogee, velocity pulses can be applied in the same or opposite direction as the velocity. It then travels along a two-body orbit to reach the next atmospheric inlet;

[0062] (4) Entering orbit: In the last orbit, i.e., the first orbit... After the spacecraft leaves the atmosphere, the terminal orbit insertion or circularization variables are calculated based on the spacecraft's atmospheric exit state and the target orbit. .

[0063] Let the number of atmospheric revolutions be... , No. The loop state variables and control variables are defined as follows:

[0064] (1)

[0065] (2)

[0066] in, For the first Cyclic state variables, , , , , , , The first The distance from the Earth's center to the Earth, relative atmospheric speed, longitude, latitude, track angle, heading angle, and mass of the spacecraft; For the first Circle control variables, , , The first The tilt angle, angle of attack, and thrust of the ring.

[0067] In this embodiment of the invention, the atmospheric flight segment employs a three-degree-of-freedom particle dynamics model that considers the Earth's rotation:

[0068] (3a)

[0069] (3b)

[0070] (3c)

[0071] (3d)

[0072] (3e)

[0073] (3f)

[0074] (3g)

[0075] in, , , , , , , These are the distance from the Earth's center to the ground, relative atmospheric speed, longitude, latitude, track angle, heading angle, and mass of the aircraft. , , , , , , These are the first derivatives of the distance from the Earth's center to the ground with respect to time, the first derivative of the relative atmospheric velocity with respect to time, the first derivative of longitude with respect to time, the first derivative of latitude with respect to time, the first derivative of the track angle with respect to time, the first derivative of the heading angle with respect to time, and the first derivative of the mass of the aircraft with respect to time, respectively. , , , , These are respectively: tilt angle, angle of attack, thrust, lift, and drag; This is the Earth's rotational angular velocity; It is the acceleration due to gravity; For engine specific impulse; This is standard gravitational acceleration. Lift. and resistance It can be represented as:

[0076] (4)

[0077] (5)

[0078] in, Atmospheric density; For reference area; , These are the lift coefficient and drag coefficient, respectively.

[0079] Stagnant heat flux rate is used as the aerodynamic thermal calculation index, and the calculation is performed in the following form:

[0080] (6)

[0081] in, The stagnation heat flux rate is the current flight state of the aircraft. These are the constant coefficients in the Chapman stagnation heat flow model.

[0082] In other embodiments, equations (3a)-(3g) can be replaced with other dynamic models that can calculate state derivatives and have first-order partial derivatives with respect to state and control, and equation (6) can be replaced with other relevant calculation formulas that can characterize the aerodynamic thermal loads of the aircraft.

[0083] Definition of the first Atmospheric inlet variables for:

[0084] (7)

[0085] Then the first Complete entry status of the circle for:

[0086] (8)

[0087] in, The distance from the Earth's center to the atmospheric boundary; superscript. This indicates transpose.

[0088] Mission constraints include initial descent constraints, control amplitude constraints, minimum flight altitude constraints, upper limit of heat flux constraints, atmospheric exit boundary constraints, inter-orbit connectivity constraints, intermediate apogee altitude constraints, intermediate orbit inclination window constraints, and terminal target orbit constraints.

[0089] The initial orbital elements are denoted as The initial descent mapping is denoted as The initial descent constraint is:

[0090] (9)

[0091] in, This is the residual vector constrained by the initial trajectory reduction equation; For the first atmospheric inlet variable; This is the initial descent velocity pulse.

[0092] No. Atmospheric flight satisfies:

[0093] (10)

[0094] in, For the first The first derivative of the loop state variable with respect to time; This represents the original nonlinear dynamic function; For time variables, from the first Starting from the atmospheric inlet; For the first The end time of atmospheric flight, i.e., the [number]th The duration of flight from atmospheric inlet to atmospheric outlet.

[0095] The control limit constraint can be expressed as:

[0096] (11)

[0097] (12)

[0098] (13)

[0099] in, , These are the preset minimum and maximum tilt angle values, respectively. , These are the preset minimum and maximum angle of attack, respectively; This is the preset maximum thrust.

[0100] The minimum flight altitude constraint is:

[0101] (14)

[0102] The upper limit constraint for heat flux is:

[0103] (15)

[0104] in, This is the preset minimum distance from the Earth's center for the aircraft; For the first During atmospheric flight The heat flux rate at any given moment; The maximum heat flux rate allowed by the pre-defined aircraft.

[0105] No. The atmospheric outlet is located at the atmospheric boundary, meaning the atmospheric outlet boundary constraint is:

[0106] (16)

[0107] in, for The distance between the spacecraft's center and the Earth at any given moment.

[0108] for , No. The exit state is integrated by atmospheric two-body dynamics and a regulating velocity pulse is applied at the apogee. After that, the first The entry state of the circle and the connection constraints between circles are:

[0109] (17)

[0110] in, For the first Circle and the first The residual vector of the inter-loop connection equation constraint is used to represent the difference between the next loop atmospheric inlet variable in the optimization variables and the next loop atmospheric inlet variable obtained from the atmospheric outlet state of the current loop and the propagation of the inter-loop adjustment velocity pulse. When the residual vector is zero, it indicates that the states of the two loops are continuous. For the first The nonlinear connection mapping between loops, which is based on the first loop... The atmospheric exit state and inter-sphere adjustment variables are taken as inputs. The process includes propagating the exit state to the apogee according to atmospheric exobody orbital dynamics, applying adjustment velocity pulses and updating the corresponding mass at the apogee, and continuing the propagation to the next atmospheric boundary. The output is the [missing information]. Variables at atmospheric inlet; For the first Variables at atmospheric inlet; for The state variable at any given time.

[0111] The intermediate apogee height constraint and the intermediate orbital inclination window constraint are as follows:

[0112] (18)

[0113] (19)

[0114] in, For the first The apogee altitude of the orbit corresponding to the atmospheric exit state; , These are the lower and upper limits of the permissible altitude of the orbital apogee corresponding to the atmospheric exit, respectively. For the reason according to the The orbital inclination angle calculated from the atmospheric exit state. To be assigned to the The orbital inclination of the intermediate target at the atmospheric exit; This represents the allowable absolute deviation between the intermediate actual track inclination angle and the intermediate target track inclination angle; .

[0115] The final target orbit constraints include atmospheric exit boundary constraints, target orbit inclination constraints, and target apogee altitude constraints, which are as follows:

[0116] (20)

[0117] (twenty one)

[0118] (twenty two)

[0119] in, For the first The distance from the Earth's center at the final moment of atmospheric flight; For the first The orbital inclination angle corresponding to the final moment of flight in the atmosphere; For the first The apogee altitude of the orbit corresponding to the final moment of flight in the atmosphere; The inclination angle of the target orbit; The apogee altitude of the target orbit.

[0120] Terminal entry or circularization variable after the last exit lap By mapping Sure:

[0121] (twenty three)

[0122] in, For the first The flight status at the final moment of atmospheric flight.

[0123] According to the Tsiolkovsky equations, the mass of the spacecraft after terminal orbit insertion is:

[0124] (twenty four)

[0125] in, The mass of the spacecraft after it finally completes its orbital insertion (i.e., the remaining mass of the spacecraft after it has applied the orbital insertion velocity pulse and consumed the required fuel, which is the mass of the spacecraft after completing the orbital change mission). For the first The mass of the spacecraft at the final moment of its atmospheric flight.

[0126] In one embodiment, minimization can be adopted. As the objective of fuel optimal trajectory planning The objective function is denoted as .

[0127] S102, Initial reference trajectory generation and discrete grid initialization.

[0128] Step S102 is used to provide an initial reference trajectory for subsequent sequence convex optimization. First, set the... The initial discrete grid order of the loop Preset initial descent velocity pulse (i.e., preset initial descent variables), the first Preset adjustment speed pulse of the circle (i.e., preset inter-circle adjustment variables), preset constant angle of attack and preset constant tilt angle .

[0129] Apply on the initial circular orbit The propagation through atmospheric exodynamics leads to the first arrival at the atmospheric boundary, resulting in the first inlet state. Initial control of atmospheric flight in each cycle can employ a fixed... , The thrust component is set in a form that compensates for the energy consumed by drag:

[0130] (25)

[0131] in, The pre-set circles for the aircraft The magnitude of thrust control at any given moment; for The resistance of time; This is the thrust amplification factor, which is used to compensate for the orbital energy loss caused by the maximum thrust limit. This is a function that takes the minimum value of each component.

[0132] Integrate the atmospheric flight state until crossing the atmospheric boundary again, then integrate the exit state to the apogee and apply... The data is then integrated up to the next atmospheric inlet; the above process is repeated to obtain the inlet status, atmospheric flight time, status history, and control history for each orbit.

[0133] The state and control history obtained from continuous propagation are mapped to the current discrete node using piecewise cubic Hermite interpolation, piecewise linear interpolation, or other interpolation methods that preserve endpoint values, to obtain the _th_ ... Initial reference node state of the loop Initial node control quantity Initial flight time and initial entry interface variables .

[0134] S103. Dimensionality reduction of variables based on MCPI fixed-point equations.

[0135] Step S103 is used to remove the node state from the optimization variables. For the first... Time normalization is performed on the atmospheric flight interval, dividing the physical time interval. Mapping to standard interval Set up Chebyshev-Gauss-Lobata nodes And construct the corresponding improved Chebyshev-Pikachu integral matrix. ;in, For the first The order of the discrete grid of the loop (i.e., the number of discrete subintervals). For the first The first circle The nodes corresponding to each discrete sub-interval .

[0136] No. Circle node state and node control quantity They are represented as follows:

[0137] (26)

[0138] (27)

[0139] in, For the first The first circle The node status at each node; For the first The first circle Node control variables at each node; .

[0140] The improved Chebyshev-Pika iterative fixed-point equation in the 1st... The first circle Node status at each node Represented as:

[0141] (28)

[0142] in, ; For the first Improved Chebyshev-Pikachu integral matrix of the circle The Line 1 Column elements; This is the original nonlinear dynamic function.

[0143] Definition of the first The node state vector of the loop Node control vector and node dynamics vector for:

[0144] (29)

[0145] (30)

[0146] (31)

[0147] in, Vectorization is represented.

[0148] definition and for:

[0149] (32)

[0150] (33)

[0151] in, For the first Improved Chebyshev-Pikachu integral matrix of the circle With the state dimension identity matrix (Right now An identity matrix of order 1. The Kronecker product (where the dimension of the state variables is ); To make the first The entry state of the loop is copied to the initial value mapping matrix of each node; For length is A column vector in which all elements are 1; This is the Kronecker product operator.

[0152] The improved Chebyshev-Pika iterative fixed-point equation is then:

[0153] (34)

[0154] in, For the first Implicit state propagation mapping of the loop.

[0155] Equation (34) is a nonlinear fixed-point equation. Given... , and Subsequently, it can be obtained through pickup truck iterations. .therefore, It is no longer treated as an independent optimization variable, but is mapped by implicit state propagation. Implicit determination.

[0156] For atmospheric flight cycles of The task is to optimize variables by reducing dimensionality. It can be written as:

[0157] (35)

[0158] in, For the first The node control vector of the loop; For the first Atmospheric flight time in the circle; This is the initial descent velocity pulse; For the first The adjustment speed pulse of the loop; For the first Atmospheric inlet variables of the sphere; .

[0159] Compared to the direct discretization method that simultaneously optimizes all node states and node control variables, equation (35) does not include... Each node has its own independent state variables, and there is no need to set dynamic equation constraints for each node in the convex optimization subproblem; among which... In the pneumatically assisted non-planar trajectory changing embodiment given in this invention, the dimension of the state variable is [not specified]. The value is 7.

[0160] S104. Construction of dimensionality reduction convex optimization subproblems.

[0161] Step S104 is used to construct a convex optimization subproblem near the current reference trajectory. Because It is no longer considered an independent optimization variable, so the objective function Equality-constrained residuals and inequality-constrained residuals This can be represented as variables related to dimensionality reduction optimization. The function, that is:

[0162] (36)

[0163] (37)

[0164] (38)

[0165] in, It consists of initial descent constraints, atmospheric exit constraints, inter-orbit connection constraints, and terminal orbit constraints; It consists of control limit constraints, height constraints, heat flow constraints, intermediate far point height constraints, and intermediate tilt angle window constraints.

[0166] In the A first-order expansion is performed at the reference point of the next inner iteration:

[0167] (39)

[0168] (40)

[0169] (41)

[0170] in, To optimize the change in variables for dimensionality reduction; The Jacobian matrix of the objective function with respect to the dimensionality reduction optimization variables; Let be the Jacobian matrix of the equality constraint residuals with respect to the dimension reduction optimization variables; The Jacobian matrix of the inequality-constrained residuals with respect to the dimensionality-reduced optimization variables; superscript Indicates the first The inner-level convex optimization subproblem solution iteration (referred to as inner-level iteration).

[0171] To avoid the linearization constraints becoming unsolvable due to the infeasibility of the initial reference trajectory, nonnegative equality slack variables are introduced. Nonnegative inequality slack variables And set the relaxation penalty weight of the equality constraint in the objective function. Inequality constraint relaxation penalty weight The objective function of the convex optimization subproblem under a fixed mesh can be expressed as:

[0172] (42)

[0173] in, For a first-order term of the objective function; This is a penalty term for slackness of equality constraints; This is a penalty term for inequality constraint relaxation; It is a column vector consisting entirely of 1s.

[0174] The constraints of this convex optimization subproblem include:

[0175] (43)

[0176] (44)

[0177] (45)

[0178] (46)

[0179] (47)

[0180] in, , These are the lower and upper bounds of the variables used for dimensionality reduction optimization, respectively. For the first The trust regions of each physical quantity in the next inner iteration include angle of attack, bank angle, thrust, flight time, initial descent velocity pulse, and inter-orbit adjustment velocity pulse. The penalty term can take the L1 norm form or other non-negative relaxation penalty forms that preserve convexity.

[0181] S105, Analytical Jacobian matrix calculation.

[0182] Step S105 is used for calculation , and In the first In the inner iteration, the first... Circle The state Jacobian matrix of each node and control Jacobian matrix :

[0183] (48)

[0184] (49)

[0185] The state Jacobian matrices of each node are combined into a block diagonal state matrix. The control Jacobian matrices of each node are combined into a block diagonal control matrix. Taking the first variation of equation (34), we get:

[0186] (50)

[0187] in, It is a first-order variational operator.

[0188] Sorted as:

[0189] (51)

[0190] (52)

[0191] (53)

[0192] (54)

[0193] in, This is the mapping matrix from atmospheric inlet variables to the complete inlet state; For the first In the next inner iteration, the improved Chebyshev-Pikka iterative fixed-point equation is used to... The state variation coefficient matrix obtained after linearizing the states of the cyclic nodes is used to represent the state variation coefficient matrix of the first node. The coupling relationship between all node state perturbations in the loop; For the first In the second inner iteration The linear propagation matrix from the control quantity perturbation of a circle node to the state perturbation of all nodes in that circle; For the first In the second inner iteration The column vector of atmospheric flight time perturbation to the state perturbation of all nodes in the loop; It is the identity matrix. Equation (51) shows that the node state disturbance can be represented by disturbances of atmospheric inlet variables, node control, and flight time.

[0194] For the Any output quantity of the loop (i.e. any physical quantity that needs to be calculated). For example, terminal tilt angle, apogee altitude, inter-loop connectivity residuals, altitude, or heat flux, its first variational... for:

[0195] (55)

[0196] in, for Partial derivatives with respect to the node state; for Partial derivative with respect to node control quantity; for Partial derivative with respect to atmospheric flight time.

[0197] Use directly This may require forming a large-scale dense inverse matrix. Therefore, for the output... Solve the discrete adjoint equation:

[0198] (56)

[0199] in, The state variation coefficient matrix The inverse matrix; For the first In the second inner iteration Circle about physical quantities An improved Chebyshev-Pika iterative discrete adjoint matrix is ​​used to represent physical quantities. The first-order derivative information regarding the nodal state is backpropagated to the atmospheric inlet variables, nodal control variables, and atmospheric flight time. When it is a scalar, It is a discrete adjoint vector.

[0200] Therefore, we obtain Jacobian matrices for atmospheric inlet variables, nodal control variables, and atmospheric flight time:

[0201] (57)

[0202] (58)

[0203] (59)

[0204] in, , , They are respectively Jacobian matrix for atmospheric inlet variables, nodal control variables, and atmospheric flight time.

[0205] Will By sequentially taking the objective function, inlet state and inter-loop connection residuals, terminal track residuals, altitude, heat flux, intermediate apogee altitude, and intermediate inclination angle as output quantities, the system can be assembled. , and For those with the same The multiple scalar outputs of the structure can be solved using linear equations with multiple right-hand sides for reuse. Matrix decomposition.

[0206] S106, Trial-step nonlinear verification and trust region update.

[0207] Step S106 is used to avoid excessively poor approximation of the local linear model. According to... Figure 2 The trial step verification and trust region update process shown is obtained by solving equation (47). , by the current reference trajectory (i.e. the first) (Dimensionality reduction optimization variables in the next inner iteration) Trial steps obtained from convex optimization subproblems Adding them together yields a trial solution. :

[0208] (60)

[0209] Define constraint violation degree and value function :

[0210] (61)

[0211] (62)

[0212] in, To constrain the weight of penalties for violations; It is the 1 norm of the vector, which is the sum of the absolute values ​​of the vector's components; This is a function that takes the maximum value of each component.

[0213] right The discrete control is piecewise continuously reconstructed, and the atmospheric flight of each sphere is propagated under the original nonlinear dynamics to calculate the actual trajectory and the actual objective function. Actual equality constraint residuals and actual inequality constraint residuals, actual minimum flight altitude Actual height of intermediate or terminal apogee Actual maximum heat flux and actual value function .

[0214] In one embodiment, the trial step is accepted when the trial solution simultaneously meets the following conditions:

[0215] (63)

[0216] (64)

[0217] (65)

[0218] (66)

[0219] (67)

[0220] in, , and A positive number used to limit the degree to which the trial solution deviates from the physical feasible region; This is the preset minimum flight altitude for the aircraft; This is the preset atmospheric boundary height; This is the preset maximum apogee altitude of the aircraft during its flight segment outside the atmosphere; The preset maximum heat flux rate of the aircraft during atmospheric flight. The value function of the current reference trajectory is obtained by using the first... Dimensionality reduction optimization variables in the next inner iteration Substituting into equation (62) yields the result. If the trial solution fails the physical feasibility check, then let... And shrink the corresponding trust region according to the following formula:

[0221] (68)

[0222] in, For the first The physical quantity in the first... The radius of the trust region in the next inner iteration; For the first The minimum value of the trust region radius corresponding to each physical quantity; Let be the shrinkage rate constant of the trust region radius. ; superscript Indicates the first Second inner iteration.

[0223] If the trial step is accepted, then calculate:

[0224] (69)

[0225] (70)

[0226] (71)

[0227] in, To predict the amount of decline; This represents the actual decrease. Model fit; To test the solution The value function, by The result is obtained by substituting into equation (62).

[0228] when Not less than the first preset threshold When the trust region corresponding to the time expansion is reached, Less than the second preset threshold The corresponding trust region shrinks under certain conditions, while the corresponding trust region remains unchanged under other conditions, i.e.:

[0229] when hour, ;

[0230] when hour, ;

[0231] In other cases, ;

[0232] in, Let be the rate constant for the expansion of the trust region radius. ; For the first The maximum value of the trust region radius corresponding to each physical quantity; Trust domain subscript This indicates that each physical quantity expands or contracts independently according to its own magnitude.

[0233] S107. Iterative convergence criteria for solving convex optimization subproblems with fixed meshes.

[0234] Step S107 is used to determine whether the sequential convex optimization on the current grid has achieved acceptable feasibility and stability. The convergence conditions for solving the convex optimization subproblem include feasibility conditions and stability conditions for adjacent iterations.

[0235] Feasibility conditions may include: terminal inclination error, terminal apogee altitude error, atmospheric exit radius error of each orbit, initial descent residual, inter-orbit connection residual, intermediate inclination window violation, minimum altitude violation, intermediate apogee altitude violation, and heat flow violation, each not exceeding the corresponding feasibility tolerance.

[0236] Stability conditions may include: the changes in nodal control variables, atmospheric flight time, initial descent variables, inter-orbit adjustment variables, and the objective function in the results of two consecutive iterations of solving subproblems in convex optimization do not exceed the corresponding tolerances, i.e.:

[0237] (72)

[0238] (73)

[0239] (74)

[0240] (75)

[0241] (76)

[0242] in, This is the normalized vector for the control quantity; , , , , These represent the nodal control variables, atmospheric flight time, initial descent variables, inter-orbit adjustment variables, and tolerances corresponding to the objective function. Each tolerance is set according to the required accuracy and physical scale. The superscript indicates... The physical quantity represents the first The physical quantities in the innermost iteration (i.e., the first) The physical quantity obtained from solving the subconvex optimization subproblem), with superscript as The physical quantity represents the first The physical quantities in the innermost iteration (i.e., the first) (physical quantities obtained by solving the subconvex optimization subproblem) The vector's infinite norm is defined as the maximum absolute value among its components. Additionally, a maximum number of inner iterations is set as a protective termination condition.

[0243] S108, Construction of Discrete Thrust Switch Function and Mesh Improvement Criteria.

[0244] Step S108 is used to update the mesh order, improving solution accuracy. According to... Figure 2 The grid update process shown defines the first... Circle Normalized thrust of each node for:

[0245] (77)

[0246] in, For the first Circle The thrust of each node.

[0247] After the inner layers of the fixed mesh converge, the dual variables corresponding to the exit radius, intermediate dip angle window, inter-cycle connectivity, intermediate apogee height, terminal trajectory, and objective function are obtained from the convex optimization subproblem. The th... Circle terminal status For discrete augmented Lagrangian functions The gradient is defined as the discrete terminal costate boundary vector. :

[0248] (78)

[0249] For the middle circle, The gradient of the terminal state with respect to the exit radius constraint, intermediate slope window, inter-loop connectivity, and intermediate apogee height constraint is obtained by combining its dual variables; for the last loop, It is obtained by combining the gradient of the objective function and the terminal trajectory constraint with respect to the terminal state with respect to its dual variables.

[0250] Extracting the matrix using terminal state Discrete terminal costate boundary vector Mapped to the node state space:

[0251] (79)

[0252] in, For the first Discrete augmented Lagrangian function generated by the loop terminal state correlation term Regarding the node status of this circle The gradient vector; express A matrix of order all zeros. express An identity matrix of order 1.

[0253] Solve the first Loop Discrete Adjoint Equation:

[0254] (80)

[0255] in, The improved Chebyshev-Pikka iterative fixed-point equation is used to find the convergent solution of the inner sequence convex optimization under the current fixed grid. The state variation coefficient matrix obtained after linearizing the state of the loop nodes; For the first The discrete adjoint variables of the loop can be obtained by solving equation (80); superscript This indicates transpose.

[0256] definition ,in, For the first Circle The partial derivative of the original nonlinear dynamic function at each node with respect to the normalized thrust is used to characterize the change in the derivative of the state variable when the normalized thrust at that node changes slightly. According to the chain rule of function differentiation, this partial derivative is equal to the partial derivative of the dynamic function with respect to the thrust multiplied by the maximum thrust. For the first Circle The original nonlinear dynamic function value at each node is calculated from the current node state variables and node control variables, i.e., the derivative of the state variable at that node. Circle Discrete thrust switching function of each node for:

[0257] (81)

[0258] in, For the first Circle Discrete adjoint variables of each node; For the first The improved Chebyshev-Pikachu integral matrix of the circle Line 1 The elements of the column.

[0259] Equation (81) reflects the current discrete optimization problem and dual variables, the first... The effect of small changes in normalized thrust at each node on the discrete augmented target can be used as an approximation of the continuous thrust switching function in the current discrete framework.

[0260] The candidate thrust switching position is determined by the position of the nodal thrust transition between shutdown and startup and the zero point position of the discrete thrust switching function. The first thrust can be extracted from the following positions. Candidate thrust switching position in the circle :

[0261] (1) Nodes where normalized thrust changes from near 0 to near 1 or from near 1 to near 0;

[0262] (2) Discrete thrust switching function Nodes close to zero.

[0263] The set of candidate switching positions is denoted as ,in, This represents the number of candidate thrust switching positions. (For candidate grid order) Construct a set consisting of discrete nodes and the midpoints of the intervals between adjacent nodes. :

[0264] (82)

[0265] in, For the first discrete nodes ; For the first The discrete node and the first The midpoint of a discrete node interval .

[0266] Definition of the first Candidate switching positions and sets Normalized alignment error between :

[0267] (83)

[0268] in, This represents the local grid spacing near the candidate switching position.

[0269] Construct the first Circle mesh quality evaluation function :

[0270] (84)

[0271] in, For the first The weights of each candidate switching position; The worst alignment error weight; Penalty weights for grid order growth; It is a function for maximizing the value.

[0272] In the allowed set of candidate grid orders Select to The minimum grid order is used as the first Circle the grid order for the next grid iteration. :

[0273] (85)

[0274] in, The function that takes the minimum value of the function corresponding to the independent variable.

[0275] Interpolate the current inner layer convergence to... The corresponding new node serves as the hot start for the next mesh iteration. The discrete thrust switching function is only used for mesh quality evaluation and mesh updates, and does not replace the inner or outer convergence conditions on its own.

[0276] S109, serial integral verification and grid-improved iterative convergence criteria.

[0277] Step S109 is used to evaluate the feasibility of the complete multi-cycle trajectory under the original continuous dynamics. The serial integration starts from the initial orbit, applies an initial descent velocity pulse, and integrates to the first orbit's atmospheric inlet state. Continuous control obtained from discrete control reconstruction is used to integrate the first orbit's atmospheric flight. The first orbit's exit state is integrated to the apogee, and an inter-orbit adjustment velocity pulse is applied. Integration continues to the second orbit's atmospheric inlet state. Thereafter, each orbit is integrated sequentially, without resetting the inlet state to the discrete optimization state at any atmospheric inlet.

[0278] Therefore, the initial orbit insertion constraints and inter-orbit connection constraints are guaranteed by a numerical integrator in the serial integration. This serial integration can reflect the cumulative impact of atmospheric segment discrete errors and continuous control reconfiguration errors on the multi-orbit maneuver mission. The terminal injection velocity pulse is recalculated based on the exit state of the last orbit, and the post-orbit mass is further calculated.

[0279] The convergence conditions for outer mesh iteration include dynamic feasibility conditions and inter-mesh variation conditions.

[0280] Dynamic feasibility conditions may include: the terminal orbit inclination error obtained by serial integration, the terminal orbit altitude error, the atmospheric exit error for each orbit, the minimum flight altitude violation, the maximum heat flux violation, and the intermediate apogee altitude violation do not exceed the corresponding tolerances, i.e.:

[0281] (86)

[0282] (87)

[0283] (88)

[0284] (89)

[0285] (90)

[0286] (91)

[0287] in, The terminal orbital inclination angle is obtained through serial integration. Let be the orbital inclination angle of the target orbit. The terminal track inclination tolerance parameter is used in the nonlinear verification of the trial step. The terminal orbital height value is obtained through serial integration. The apogee altitude of the target orbit. The terminal track height tolerance parameter is used in the nonlinear verification of the trial step. The terminal geocentric distance value obtained by serial integration. The terminal geocentric distance tolerance parameter is used in the nonlinear verification of the trial step. The minimum flight altitude value obtained by serial integration. The minimum permitted flight altitude for an aircraft. For the flight altitude tolerance parameters in the nonlinear verification of the trial step; To obtain the apogee altitude value of the spacecraft's orbit through serial integration, This represents the maximum permissible apogee altitude for the aircraft. To test the nonlinearity of the orbital apogee height tolerance parameters; To obtain the first step by serial integration The maximum heat flux rate of a spacecraft during atmospheric flight. This is a parameter for the heat flux rate tolerance during the trial nonlinear verification.

[0288] Inter-grid variation conditions are used to describe the requirements between two adjacent grid iterations, including: the changes in the intermediate apogee altitude, inclination angle, atmospheric exit state, atmospheric flight time, initial descent velocity pulse, inter-circle adjustment velocity pulse, and objective function of the converged solutions between two adjacent grid iterations do not exceed the corresponding tolerances.

[0289] The outer mesh is considered to have converged only when both the dynamic feasibility condition and the inter-mesh variation condition are met.

[0290] The final output trajectory planning results include: atmospheric inlet and outlet conditions for each orbit, angle of attack curves, side slope curves and thrust curves for each orbit, discrete thrust function curves, flight time for each orbit, minimum altitude, maximum heat flux, intermediate apogee altitude, initial descent variables, inter-orbit adjustment variables, terminal orbit insertion variables, and post-terminal orbit mass.

[0291] The following set of parameters is given to illustrate the implementation process of the present invention, but the present invention is not limited to this set of parameters.

[0292] This embodiment uses three-circle aerodynamic-assisted non-planar orbit change fuel optimal trajectory planning as an example. The Earth's gravitational constant is taken as... Earth's radius is taken Earth's rotational angular velocity is taken as Atmospheric boundary height is taken The initial mass of the aircraft is taken as... Reference area is taken Engine specific impulse The coefficients of the heat flow model are taken as follows: The initial number of circular orbital elements is set to... In the orbital element vector, the components are, in order: orbital semi-major axis, orbital eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and true anomaly. The target orbital inclination and altitude are respectively taken as... and That is, the corresponding total change in tilt angle is The number of orbital changes is taken .

[0293] The mission constraint parameters are set as follows: the allowable range for the initial descent velocity pulse is... Atmospheric flight time in each circle meets The velocity pulse adjustment between far points of the circle satisfies The range of tilt angles is taken as follows: Angle of attack range Maximum thrust The apogee altitude limit for the middle circle is... The allowable deviation of the intermediate target tilt angle is taken as Minimum flight altitude is taken The maximum heat flux rate is taken .

[0294] In the sequential convex optimization subproblem, the penalty weights for the slack variables of both equality and inequality constraints are taken as... The maximum number of iterations for inner sequence convex optimization is taken as... The maximum number of iterations for the outer mesh is taken as The initial discrete order for each lap is taken as... ,correspond There are Chebyshev-Gauss-Lobato nodes, and the initial number of nodes for all three orbits is [number missing]. The worst alignment error weight in the grid evaluation function is taken as follows: Grid order penalty weights are taken .

[0295] Control variables are normalized vectors The upper and lower bounds of the trust region for the normalized control variables are respectively taken as follows: and The upper and lower limits of the flight time trust region are respectively taken as follows: and The initial descending pulse trust region upper and lower limits are respectively taken as follows: and The upper and lower limits of the inter-loop adjustment pulse trust region are respectively taken as follows: and .

[0296] Trust region expansion and contraction factors are respectively taken as and The threshold for model conformity is taken as follows: and The trial step nonlinear accepts the parameter selection , and .

[0297] The initial trust region radii are respectively set to , , and .

[0298] In the convergence criterion for inner subproblems, the normalized residual tolerances for both the entry constraint and the inter-loop connectivity constraint are taken as... Atmospheric exit radius tolerance is taken The tolerances for terminal tilt angle and terminal apogee height are respectively taken as follows: ° and The tolerances for the inclination angle, height, and heat flux inequality constraint of the middle ring are respectively taken as follows: ° , and The update tolerances for control variables, flight time, initial descent pulse, inter-orbit adjustment pulse, and objective function are respectively set to... , , , and .

[0299] In the trial-step nonlinear verification, the minimum height tolerance is taken as... The tolerances for apogee altitude and atmospheric exit radius are both taken as... Terminal tilt angle tolerance is taken ° Heat flux tolerance is taken In the outer mesh convergence criterion, the tolerance for variations in apogee height and exit radius between adjacent generations is taken as... The tolerance for changes in the initial orbital descent velocity pulse and the inter-cycle adjustment velocity pulse are both taken as... The tolerance for tilt angle variation is taken as follows ° The tolerance for the variation in flight time for each orbit is taken as follows: The tolerance for changes in the objective function is taken as follows: When using trial steps for integration, an adaptive step-size integrator is employed, with the relative and absolute error tolerances of the integrator set as follows: and In serial integration verification, the relative error and absolute error tolerance of the integrator are respectively set as follows: and .

[0300] In the initial reference trajectory construction, the initial bank angle and angle of attack for each orbit are fixed as follows: and The initial thrust value is based on The settings include the thrust amplification factor. Initial descent pulse guessing Initial guess of the apogee adjustment pulse between each circle The orbit is initially induced by a descent pulse, which propagates until the first atmospheric boundary is reached. Subsequently, atmospheric flight, extra-atmospheric flight, and final orbit insertion are completed sequentially. The resulting state and control variables are mapped to the Chebyshev-Gauss-Lobart nodes of each orbit using piecewise cubic Hermite interpolation. The convex optimization subproblem is solved using the CVX framework and the MOSEK solver.

[0301] After calculation using the above parameters, the algorithm reached the set convergence condition after four grid iterations. The obtained control quantities, trajectory quantities, and discrete thrust switching function results can be used as the calculation results diagram of this embodiment. The number of discrete control nodes for the three orbits are 199, 89, and 113, respectively, with a total of 401 nodes. The initial descent velocity pulse is -98.14 m / s, the inter-orbit adjustment velocity pulses are -6.77 m / s and -5.11 m / s, respectively, and the terminal circularization velocity pulse is 91.06 m / s. The mass of the spacecraft after terminal circularization is approximately 52794.18 kg. The terminal tilt angle error is approximately 1.03 × 10⁻⁶ m / s. -4 The terminal apogee altitude error is approximately The atmospheric exit radius error is approximately .

[0302] Figure 4 The diagram shows the output curve and discrete thrust switching function curve obtained using the method of this embodiment of the invention. Figure 4 In the figure, (a), (b), (c), (d), (e), and (f) represent the inclination angle curve, angle of attack curve, thrust normalized variable curve, altitude curve, heat flux rate curve, and right ascension curve of ascending node for orbits 1, 2, and 3, respectively. Figure 4 In this context, (g) represents the thrust normalized variable curve and discrete thrust switching function curve corresponding to the first, second, and third atmospheric flights. For example... Figure 4 As shown, the thrust curves of each revolution switch between 0 and the maximum thrust, exhibiting an approximate on / off control characteristic. The proportion of nodes where the discrete thrust and the discrete thrust switching function are inconsistent decreases from about 21.95% in the initial grid to about 2.49% in the final grid, indicating that the grid improvement method of the present invention improves the grid's ability to discretely represent the thrust switching structure.

[0303] Under the same dynamics, constraints, and initial guesses, compared to a comparative implementation based on Legendre-Gauss-Rado collocations and HP adaptive meshes, the number of discrete control nodes used in this embodiment is reduced from 782 to 401, a reduction of approximately 48.72%; the total computation time is reduced from approximately 848.61 s to approximately 274.78 s, a reduction of approximately 67.62%. These comparative results illustrate that under the stated implementation parameters and initialization conditions, the present invention can reduce the discrete scale and computational overhead while maintaining the feasibility of continuous dynamics, and do not constitute a limitation on performance under all task parameters.

[0304] In other optional embodiments of the present invention, the discrete nodes can be other collocations, as long as a corresponding integral matrix can be constructed and a fixed-point equation can be formed, they can all be used in the present invention.

[0305] In other alternative embodiments of the invention, the convex optimization subproblem can be solved by a solver that supports linear programming, second-order cone programming or other convex programming forms, and the slack variable penalty form and trust region norm can be replaced with equivalents while maintaining convexity.

[0306] In other alternative embodiments of the present invention, different variable-step-size high-precision ordinary differential equation integrators can be used in nonlinear verification and serial integration.

[0307] In other optional embodiments of the present invention, the objective function may be replaced by task time, thermal load or other indicators, and the corresponding output Jacobian matrix of the dimension reduction optimization variables can still be constructed by equations (56) to (59).

[0308] In other optional embodiments of the present invention, different atmospheric density models, aerodynamic models and dynamic models can be used for the atmospheric flight segment. As long as the dynamic function and its first-order partial derivatives with respect to the state and control are available, the variable dimensionality reduction and analytical Jacobian matrix calculation method of the present invention can be applied.

[0309] The present invention also provides an electronic device, comprising: one or more processors; and a memory for storing one or more programs; wherein, when the one or more programs are executed by the one or more processors, the one or more processors implement the above-described method for dimensionality-reduced sequential convex optimization trajectory planning for multi-loop aerodynamic-assisted heterogeneous orbit changes of aerospace vehicles.

[0310] The present invention also provides a computer-readable storage medium storing executable instructions thereon, which, when executed by a processor, enable the processor to implement the above-described method for dimensionality-reduced sequential convex optimization trajectory planning for multi-loop aerodynamic-assisted heterogeneous orbit changes of aerospace vehicles.

[0311] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. The solutions in the embodiments of the present invention can be implemented using various computer languages, such as the object-oriented programming language Java and the interpreted scripting language JavaScript.

[0312] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0313] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0314] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0315] The contents not described in detail in this specification are existing technologies known to those skilled in the art.

[0316] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the present invention. The above descriptions are exemplary and not exhaustive. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted non-planar orbit changes of aerospace vehicles, characterized in that, include: Step 1: Obtain relevant parameters for the multi-cycle aerodynamic-assisted orbital maneuvering mission. Based on these parameters, establish a mathematical model for the optimal control problem of the multi-cycle aerodynamic-assisted orbital maneuvering mission, which includes the initial descent, atmospheric flight, extra-atmospheric flight, and final orbit insertion processes. The relevant parameters include initial orbital parameters, target orbital parameters, aerospace vehicle parameters, mission constraint parameters, and the number of atmospheric flight cycles. Step 2: Generate the initial reference trajectory for each orbit based on the preset initial descent variables, preset inter-orbit adjustment variables, preset constant angle of attack, and preset constant tilt angle for each orbit; the initial reference trajectory for each orbit includes the initial reference values ​​of atmospheric inlet variables, node status, node control variables, and atmospheric flight time for each orbit. Step 3: In the atmospheric flight segment of each orbit, establish an improved Chebyshev-Pika iterative fixed-point equation at discrete nodes, and represent the node state as an implicit function of atmospheric inlet variables, node control variables, and atmospheric flight time. Construct dimension reduction optimization variables based on the node control variables, atmospheric flight time, initial orbit reduction variables, inter-orbit adjustment variables, and atmospheric inlet variables of each orbit. Step 4: Calculate the objective function, equality constraint residual, and inequality constraint residual at the current reference trajectory. Make a first-order approximation of the objective function, equality constraint residual, and inequality constraint residual with respect to the dimension reduction optimization variables. Introduce non-negative equality relaxation variables, non-negative inequality relaxation variables, and trust regions to construct a convex optimization subproblem. Step 5: Obtain the first-order variation of the improved Chebyshev-Pika iterative fixed-point equation to obtain the relationship between the node state perturbation and the dimension reduction optimization variable perturbation. Construct a discrete adjoint equation, and without explicitly calculating the matrix inverse, analytically obtain the Jacobian matrix of the objective function, the equality constraint residual, and the inequality constraint residual with respect to the dimension reduction optimization variable; Step 6: Solve the convex optimization subproblem to obtain the trial step, and sum the trial step with the current reference trajectory to form the trial solution; Interpolating the discrete control quantities of the trial solution and performing numerical integration under the original nonlinear dynamics, the aircraft trajectory corresponding to the trial solution is obtained. Then, the physical feasibility index, objective function, constraint violation degree and value function corresponding to the trajectory are calculated. Then, it is determined whether the trial step is accepted based on whether the trajectory is within the set physical feasibility region. The ratio between the decrease in the value function predicted by the linearized model and the decrease in the actual value function of the original nonlinear dynamics is calculated, and the size of the trust region is updated based on the ratio. Step 7: Based on the physical feasibility of the trajectory corresponding to the trial solution, and the changes in node control variables, atmospheric flight time, initial descent variables, inter-circle adjustment variables, and objective function in the solutions of the two adjacent convex optimization subproblems, determine whether the solution of the convex optimization subproblem has converged. Step 8: When the solution of the convex optimization subproblem under the fixed grid order satisfies the convergence condition, construct the discrete thrust switching function for each discrete node based on the dual variables of the convex optimization subproblem, the terminal augmented gradient of each cycle, and the improved Chebyshev-Pika iterative fixed-point equation; determine the candidate thrust switching position based on the node thrust jump position and the zero of the discrete thrust switching function, and update the discrete grid with the alignment error between the candidate thrust switching position and the candidate new grid node or the midpoint of the adjacent node interval as the index; Step 9: Interpolate the converged solution of the convex optimization subproblem corresponding to the previous generation grid to the updated discrete grid as a warm start, and repeat steps 1 to 8; For the converged solution of the convex optimization subproblem obtained in each grid iteration, start from the initial trajectory and use serial integration to obtain the true trajectory. Based on the terminal constraints, path constraints obtained after serial integration, as well as the converged solutions and changes in the objective function of two adjacent grid iterations, determine whether the outer grid iteration has converged; Output the trajectory planning result when the outer grid iteration converges.

2. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 1, characterized in that, Step 2 includes: The preset initial descent variable is applied on the initial orbit, and integration is performed according to the two-body orbit dynamics until the atmospheric boundary is reached for the first time. In atmospheric flight, the preset constant angle of attack and the preset constant side angle are used, and the thrust component is set in the form of compensating for the energy consumed by drag. The thrust is integrated according to atmospheric segment dynamics until the atmospheric exit is reached. Starting from the atmospheric exit, integration is performed according to the two-body orbital dynamics until the orbital apogee is reached. The preset inter-cycle adjustment variable is applied, and integration is performed according to the two-body orbital dynamics until the next atmospheric inlet is reached. The atmospheric flight and extra-atmospheric flight are repeated until the last atmospheric flight is completed. The integrated state and control variables are interpolated to the current discrete node to obtain the initial reference values ​​of atmospheric inlet variables, node states, node control variables, and atmospheric flight time for each orbit.

3. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 1, characterized in that, The first The atmospheric inlet variable of the circle is denoted as Then the first The complete entry status of the circle for: ； in, The distance from the Earth's center to the atmospheric boundary; superscript. Indicates transpose; No. Circle node state and node control quantity They are represented as follows: ； ； in, For the first The order of the discrete grid of the loop; For the first The first circle The node status at each node; For the first The first circle Node control variables at each node; .

4. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 3, characterized in that, In step 3, the improved Chebyshev-Pika iterative fixed-point equation is in the... The first circle Node status at each node Represented as: ； in, , For the first Atmospheric flight time in the circle; For the first The improved Chebyshev-Pikachu integral matrix of the circle Line number Column elements; This represents the original nonlinear dynamic function; The improved Chebyshev-Pika iterative fixed-point equation is as follows: ； in, For the first The node state vector of the loop; For the first The node control vector of the loop; For the first The dynamic vector of the loop node; For the first The improved Chebyshev-Pika integral matrix of the loop and the Kronecker product of the state dimension identity matrix; To make the first The entry state of the loop is copied to the initial value mapping matrix of each node; For the first Atmospheric inlet variables of the sphere; For the first Implicit state propagation mapping of the loop.

5. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 4, characterized in that, In step 4, the convex optimization subproblem uses the change in the dimensionality reduction optimization variable, the non-negative equality slack variable, and the non-negative inequality slack variable as variables. The constraints of the convex optimization subproblem include: ； ； ； ； ； in, To optimize variables for dimensionality reduction; For equality-constrained residuals; The residuals are constrained by inequality. Let be the Jacobian matrix of the equality constraint residuals with respect to the dimension reduction optimization variables; The Jacobian matrix of the inequality-constrained residuals with respect to the dimensionality-reduced optimization variables; superscript Indicates the first Second inner iteration; To optimize the change in variables for dimensionality reduction; , These are the lower and upper bounds of the variables used for dimensionality reduction optimization, respectively. These are non-negative equality slack variables; These are slack variables for nonnegative inequalities; For the first Trust regions of each physical quantity in the next inner iteration.

6. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 5, characterized in that, Step 5 includes: In the In the inner iteration, the first... Circle the state Jacobian matrix and control Jacobian matrix of each node, and then combine the state Jacobian matrices of each node into a block diagonal state matrix. The control Jacobian matrices of each node are combined into a block diagonal control matrix. ; Taking the first variation of the improved Chebyshev-Pika iterative fixed-point equation and rearranging it, we get: ； in: ； ； ； in, For the first In the next inner iteration, the improved Chebyshev-Pikka iterative fixed-point equation is used to... The state variation coefficient matrix obtained after linearizing the states of the cyclic nodes is used to represent the state variation coefficient matrix of the first cyclic node. The coupling relationship between all node state perturbations in the loop; For the first In the second inner iteration The linear propagation matrix from the control quantity perturbation of a circle node to the state perturbation of all nodes in that circle; For the first In the second inner iteration The column vector of atmospheric flight time perturbation to the state perturbation of all nodes in the loop; It is the identity matrix; , The first In the next inner iteration and ; For any physical quantity that needs to be calculated Calculate separately Partial derivatives of node state, node control variables, and atmospheric flight time , , And solve the discrete adjoint equation: ； in, For the first In the second inner iteration Circle about physical quantities An improved Chebyshev-Pika iterative discrete adjoint matrix is ​​used to represent physical quantities. The first derivative information of the node state is back-propagated to the atmospheric inlet variables, node control variables, and atmospheric flight time; according to , , , respectively obtained Jacobian matrix for atmospheric inlet variables, nodal control variables, and atmospheric flight time.

7. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 6, characterized in that, In step 6, the constraint violation degree is defined. and value function for: ； ； in, To optimize variables for dimensionality reduction; For equality-constrained residuals; The residuals are constrained by inequality. The objective function is... To constrain the weight of penalties for violations; It is the 1-norm of the vector; This is a function that takes the maximum value of each component. When the trial step is accepted, the predicted descent amount is calculated. Actual decrease Model fit : ； ； ； in, The value function of the current reference trajectory; The value function for exploring the solution; It is a function of actual value; when Expand the corresponding trust region when it is not less than the first preset threshold. When the threshold is less than the second preset threshold, the corresponding trust region is shrunk; otherwise, the corresponding trust region remains unchanged. The first preset threshold is greater than the second preset threshold.

8. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 7, characterized in that, In step 7, the convergence conditions for solving the convex optimization subproblem include feasibility conditions and adjacent iteration stability conditions; The feasibility conditions include: terminal inclination error, terminal apogee altitude error, atmospheric exit radius error of each orbit, initial descent residual, inter-orbit connection residual, intermediate inclination window violation, minimum altitude violation, intermediate apogee altitude violation, and heat flow violation do not exceed the corresponding tolerances; The stability conditions include: the changes in nodal control variables, atmospheric flight time, initial descent variables, inter-orbit adjustment variables, and objective function in the solutions to two consecutive convex optimization subproblems do not exceed the corresponding tolerances.

9. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 8, characterized in that, In step 8, No. Circle Discrete thrust switching function of each node for: ； in, For the first Circle Normalized thrust at each node; For the first Circle Discrete adjoint variables at each node; For the first The improved Chebyshev-Pikachu integral matrix of the circle Line number Column elements; , For the first Circle The biased derivative of the original nonlinear dynamic function at each node with respect to the normalized thrust is used to characterize the change in the derivative of the state variables when the normalized thrust at that node changes slightly. For the first Circle The original nonlinear dynamic function values ​​at each node are calculated from the current node state and node control variables; For candidate grid order Construct a set consisting of discrete nodes and the midpoints of the intervals between adjacent nodes. , define the first Candidate thrust switching positions and sets The normalized alignment error between them is and construct the first Circle mesh quality evaluation function : ； in, For the first The weights of each candidate switching position; The worst alignment error weight; Penalty weights for grid order growth; It is a function for maximizing the value; Select the first grid order from the allowed set of candidate grid orders. Circle mesh quality evaluation function The minimum grid order is used as the first Circle the grid order for the next grid iteration.

10. The dimension-reduced sequential convex optimization trajectory planning method for multi-loop aerodynamic-assisted heterogeneous orbit change of aerospace vehicles according to claim 9, characterized in that, In step 9, the serial integration starts from the initial orbit, and after the initial descent process, the first atmospheric inlet state is obtained. Then, the atmospheric flight trajectory and the extra-atmospheric flight trajectory of each orbit are obtained by sequential integration. At the inlet of each orbit, the state variables are not reset to the optimized atmospheric inlet variables. The convergence conditions for the outer mesh iteration include dynamic feasibility conditions and inter-mesh variation conditions; The dynamic feasibility conditions include: the terminal orbit inclination error, terminal orbit altitude error, atmospheric exit error for each orbit, minimum flight altitude violation, maximum heat flux violation, and intermediate apogee altitude violation obtained by serial integration do not exceed the corresponding tolerances; The inter-grid variation conditions include: the changes in the mid-apex altitude, inclination, atmospheric exit state, atmospheric flight time, initial descent velocity pulse, inter-circle adjustment velocity pulse, and objective function of the converged solutions of two adjacent grid iterations do not exceed the corresponding tolerances.

Images