Robust trajectory optimization method for rarefied atmosphere pneumatic auxiliary orbit reduction
By converting random ordinary differential equation systems into deterministic ordinary differential equation systems and using orthogonal point selection strategies and pseudo-spectral methods, the problem of high computational cost in robust trajectory optimization is solved, and efficient robust trajectory generation is achieved.
Patent Information
- Application Number
- CN202510349138.0
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-03-24
- Publication Date
- 2025-07-04
AI Technical Summary
The prior art has high computational cost in robust trajectory optimization, resulting in low optimization efficiency, which can lead to significant errors when using linear covariance methods in high nonlinear scenarios.
Through Riemann–Stieltjes integral, a random system of differential equations is converted into deterministic system of ordinary differential equations, an extended process constraint and cost function is established, and an orthogonal point selection strategy is proposed, the integral function is discrete, and a standard trajectory optimization model suitable for pseudo-spectral method is constructed, and this model is solved using GPOPS tool.
The calculation cost of optimization problems under high-dimensional uncertain parameters is reduced, the efficiency of robust trajectory optimization is improved, and a highly robust flight trajectory is generated.
Smart Images
Figure CN120258270A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the field of space exploration, and particularly to a robust trajectory optimization method for aerodynamic-assisted orbit descent in rarefied atmosphere. Background Art
[0002] The basic idea of solving the robust trajectory optimization problem is to transform the stochastic optimal control problem dominated by stochastic ordinary differential equations (ODEs) into a deterministic problem dominated by deterministic ODEs, which is called dynamic uncertainty propagation. Linear covariance analysis (LCA) is an effective method for solving the uncertainty propagation problem and has the advantage of high computational efficiency. In previous studies, LCA and Legendre pseudospectral method were used in the robust trajectory design of the reentry atmosphere process, and the different robust performances of the aircraft with feedback when using open-loop covariance and closed-loop covariance for shaping were analyzed. LCA and non-dominated sorting genetic algorithm were used to solve the time-fixed linearized impulsive rendezvous optimization problem with control uncertainty. The idea of the LCA method is to linearize the dynamic system, so it may lead to significant errors in highly nonlinear scenarios. To solve this problem, a framework based on the generalized polynomial chaos (PC) theory was developed to solve the optimal trajectory generation problem with uncertain parameters in dynamic systems. Although the PC-based extended method has been widely used in the robust trajectory optimization problem, its computational cost is always high, resulting in low optimization efficiency. Summary of the Invention
[0003] To solve the above problems, the present invention provides a robust trajectory optimization method for aerodynamic-assisted orbit descent in rarefied atmosphere. By using Riemann–Stieltjes integral, the stochastic ordinary differential equation system is transformed into a deterministic ordinary differential equation system, an extended process constraint and cost function are established, and a point selection strategy is proposed to further discretize the integral function in the model, and a standard trajectory optimization model suitable for the pseudospectral method is constructed. Finally, the GPOPS tool is used to solve this standard optimization model. This method is applicable to the robust trajectory generation of dynamic systems with multi-dimensional parameter uncertainties in the scenario of aerodynamic-assisted orbit transfer.
[0004] To achieve the above object, the present invention provides a robust trajectory optimization method for aerodynamic-assisted orbit descent in rarefied atmosphere, including the following steps:
[0005] Step S1: Establish a deterministic dynamics model;
[0006] Step S2: Determine the optimization objective, which is to maximize the energy change passing through the sphere of influence of the planet;
[0007] Step S3: Establish deterministic constraint functions, including altitude constraint and heating rate constraint:
[0008] Altitude constraint: The altitude of the spacecraft is not lower than the minimum allowable altitude:
[0009] H = r - R0 ≥ H min
[0010] where r is the Cartesian position of the spacecraft in the inertial frame of the planet's center, and R0 is the radius of the planet;
[0011] Heating rate constraint: The total heating rate of the spacecraft does not exceed the allowed maximum value:
[0012]
[0013] where is the total heating rate, and are the convective and radiative heating rates respectively, and are calculated as follows:
[0014]
[0015] where C c represents a constant dependent on the atmospheric composition, r n represents the stagnation region radius, C r is a constant dependent on the planet's atmosphere, F(v) is a look-up function dependent on the magnitude of the flight speed and the atmospheric composition, and the exponents b and c are determined by ρ and v;
[0016] Step S4: Establish a deterministic standard optimization model;
[0017] Step S5: Establish an optimization model including uncertainties based on the RS integral theory;
[0018] Step S6: Establish a standard robust trajectory optimization model based on the orthogonal sampling strategy;
[0019] Step S7: Solve the optimization model in Step S6 using the pseudospectral method.
[0020] Preferably, in Step S1, the dynamic equation for planar atmospheric flight without considering thrust and planetary rotation is:
[0021]
[0022] where the orbital state x = [r, v], u is the control variable, and the lift coefficient C L is used as the control variable. The specific calculation formula is as follows:
[0023]
[0024] where r = (x, y) and v = (v x , v yare the Cartesian position and velocity states of the spacecraft in the planet-centered inertial frame, and μ represents the planet's mass multiplied by the gravitational constant. L and D represent the accelerations generated by lift and drag, respectively, expressed as:
[0025]
[0026] where S is the aerodynamic reference area of the spacecraft, m is the mass of the spacecraft, v n = (-v y , v x ) is the normal vector perpendicular to the velocity vector v, v is the magnitude of the velocity ||v||, and σ is the motion direction index of the entry orbit relative to the planet: prograde and retrograde correspond to 1 and -1, respectively;
[0027] The drag coefficient C D and the lift coefficient C L have the following relationship:
[0028]
[0029] where:
[0030]
[0031] (L / D) max represents the maximum lift-to-drag ratio, is the lift coefficient at (L / D) max ;
[0032] Regarding the atmospheric modeling, the density ρ is approximately an exponential function within the atmosphere, defined by the equation ρ = ρ0exp(-βH), where β represents the reciprocal of the standard altitude of the atmosphere, ρ0 is the reference density at sea level, and H is the altitude of the spacecraft.
[0033] Preferably, in step S2, the energy of the spacecraft in the sun-planet system is defined as:
[0034]
[0035] where V is the magnitude of the heliocentric velocity of the spacecraft, μ s is the gravitational constant multiplied by the mass of the sun, r s represents the distance between the spacecraft and the sun, and the energy change through the sphere of influence of the planet is expressed as:
[0036] ΔE = E + - E -
[0037] where E + and E - represent the energy values at the exit and entrance of the sphere of influence of the planet, respectively.
[0038] Preferably, in step S4, a deterministic AGA trajectory optimization model is constructed:
[0039] Min.J[x,u]=-|ΔE|
[0040]
[0041] H min -H(x(t),C L (t))≤0。
[0042] Preferably, in step S5, the uncertain dynamic system is described as:
[0043]
[0044] where p represents the M-dimensional uncertain parameter of the dynamic model, and the objective function is defined as the integral of the original cost function with respect to the uncertain parameter p:
[0045] ∫…∫-|ΔE(p)|dα(p)
[0046] where ΔE(p) represents the cost function of energy change, which is a function of the value of the uncertain parameter, is the known joint probability distribution function CDF of p, the multiplicity of the integral corresponds to the dimension of the uncertainty parameter, and the integration region is on the probability distribution domain of p;
[0047] Inspired by the construction of the objective function, the Riemann–Stieltjes optimization RSO model of the dynamic system containing uncertain parameters is constructed as follows:
[0048] Min.J[x(·,·),C L (·)]=∫ dom(p) -|ΔE(p)|dα(p)
[0049]
[0050] ∫ dom(p) (H min -H(x(t),C L (t);p))dα(p)≤C2
[0051] where C1, C2 represent constants, dom(p) represents the probability distribution domain of p, and ∫ dom(p) is a shorthand notation for multiple integrals.
[0052] Preferably, in step S6, the RSO model obtained in step S5 is discretized. First, the integral of the original cost function of the objective function with respect to the uncertain parameter p is approximated as a finite sum form, and the optimized objective function is obtained as follows:
[0053]
[0054] Among them n represents the set of values of the uncertain parameters with a definite probability density function PDF, and W i represents the weight coefficient after product discretization, is the joint CDF with respect to the uncertain parameter p, and its differential dα(p) is the joint PDF of p. W i is the value of the joint PDF of the uncertain parameters corresponding to The constraint function in fractional form is discretized into the following form:
[0055]
[0056] In the AGA optimization problem, two independent random variables ρ0 and φ = S / m are considered. The point selection strategy is introduced with a two-dimensional Gaussian distribution. For the two-dimensional uncertain parameter p = (p1, p2) that satisfies the Gaussian distribution, three points are uniformly selected on the axis of the uncertain parameter p k centered on the mean value μ k with the standard deviation σ k as the boundary, where k = 1, 2,..., M, M = 2, to obtain 2M + 1 = 5 sampling points. The weight coefficients of the sampling points are obtained through the joint PDF of the uncertain parameters, as follows:
[0057]
[0058] where μ represents the mean value of the uncertain parameters, and Σ represents its covariance matrix. The AGA optimization problem model of dynamic uncertainty is rewritten in the standard form:
[0059]
[0060] Preferably, in step S7, the continuous optimal control problem solver GPOPS-II is used to solve the AGA optimization problem model of dynamic uncertainty obtained in step S6.
[0061] Therefore, the present invention adopts the above-mentioned method for robust trajectory optimization of aeroassisted orbit descent in a rarefied atmosphere. By introducing the uncertainty of parameters into the deterministic optimization model through the properties of Riemann-Stieltjes, a highly robust flight trajectory can be obtained starting from the original optimization problem. And it is discretized into the standard form through the orthogonal uniform point selection strategy, reducing the computational cost of the optimization problem under high-dimensional uncertain parameters and improving the optimization efficiency.
[0062] Next, through the accompanying drawings and embodiments, the technical solutions of the present invention will be further described in detail. Description of the Drawings
[0063] Figure 1 It is a flowchart of a robust trajectory optimization method for aeroassisted orbit lowering in a tenuous atmosphere according to the present invention;
[0064] Figure 2 It is a schematic diagram of the trajectory of AGA in the sphere of influence of a planet in an embodiment of the present invention;
[0065] Figure 3 It is a schematic diagram of the two-dimensional uncertain parameter selection point strategy in an embodiment of the present invention;
[0066] Figure 4 It is the control law C obtained from the RSO result in an embodiment of the present invention L ;
[0067] Figure 5 It is a diagram showing the dispersion of the trajectory altitude and speed under the RSO result in an embodiment of the present invention. Detailed implementation manners
[0068] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.
[0069] Unless otherwise defined, the technical terms or scientific terms used in the present invention should have the ordinary meanings understood by those of ordinary skill in the art to which the present invention belongs.
[0070] The terms "including" or "comprising" and the like used in the present invention are intended to mean that the elements before this word cover the elements listed after this word, and do not exclude the possibility of also covering other elements. The orientation or positional relationship indicated by the terms "inside", "outside", "above", "below", etc. is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation of the present invention. When the absolute position of the object being described changes, the relative positional relationship may also change accordingly. In the present invention, unless otherwise clearly defined and limited, terms such as "attached" should be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or integrated; it can be directly connected, or indirectly connected through an intermediate medium, and can be the internal communication of two elements or the interaction relationship between two elements. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.
[0071] Embodiment
[0072] As Figure 1 shown, a robust trajectory optimization method for aeroassisted orbit lowering in a tenuous atmosphere takes the aero-gravitational assist orbit in the Mars scenario as an example. Table 1 lists the data of the constants, spacecraft characteristics and boundary values used to optimize the Mars AGA orbit.
[0073] Table 1
[0074]
[0075] Specifically, it includes the following steps:
[0076] Step S1: Establish a deterministic dynamics model;
[0077] As Figure 2 shown, in Step S1, the dashed trajectory represents entering a hyperbolic orbit, which is an orbit without aerodynamic or propulsive forces. The entry orbit when the spacecraft approaches the planet is a hyperbolic orbit, where ψ p represents the phase angle between the hyperbolic axis of the entry orbit and the line connecting the sun and the planet, R0 is the radius of the planet, and the outermost circular dashed line represents the radius of the sphere of influence (SOI) of the planet. The dynamic equation for planar atmospheric flight without considering thrust and planetary rotation is:
[0078]
[0079] where the orbital state x = [r, v], u is the control variable, and the lift coefficient C L is used as the control variable. The specific calculation formula is as follows:
[0080]
[0081] where r = (x, y) and v = (v x , v y ) are the Cartesian position and velocity states of the spacecraft in the inertial frame of the planet's center respectively, μ represents the planet's mass multiplied by the gravitational constant. L and D represent the accelerations generated by lift and drag respectively, and are expressed as:
[0082]
[0083] where S is the aerodynamic reference area of the spacecraft, m is the mass of the spacecraft, v n = (-v y , v x ) is the normal vector perpendicular to the velocity vector v, v is the magnitude of the velocity ||v||, and σ is the motion direction index of the entry orbit relative to the planet: prograde and retrograde correspond to 1 and -1 respectively;
[0084] The drag coefficient C D and the lift coefficient C L have the following relationship:
[0085]
[0086] where:
[0087]
[0088] (L / D) max represents the maximum lift-to-drag ratio, and is (L / D) max when the lift coefficient is;
[0089] Regarding the atmospheric modeling, the density ρ is approximately an exponential function within the atmosphere and is defined by the equation ρ = ρ0exp(-βH), where β represents the reciprocal of the standard height of the atmosphere, ρ0 is the reference density at sea level, and H is the altitude of the spacecraft.
[0090] Step S2: Determine the optimization objective, which is to maximize the energy change through the planetary sphere of influence; in Step S2, the energy of the spacecraft in the sun - planet system is defined as:
[0091]
[0092] where V is the magnitude of the heliocentric velocity of the spacecraft, μ s is the gravitational constant multiplied by the mass of the sun, r s represents the distance between the spacecraft and the sun, Figure 2 and the energy change through the planetary sphere of influence is expressed as:
[0093] ΔE = E + -E -
[0094] where E + and E - represent the energy values at the exit P3 and the entrance P0, respectively.
[0095] Step S3: Establish deterministic constraint functions, including altitude constraint and heating rate constraint:
[0096] Altitude constraint: The altitude of the spacecraft is not lower than the allowed minimum altitude:
[0097] H = r - R0 ≥ 50 km
[0098] where r is the Cartesian position of the spacecraft in the inertial frame of the planet's center, and R0 is the radius of the planet;
[0099] Heating rate constraint: The total heating rate of the spacecraft does not exceed the allowed maximum value:
[0100]
[0101] where, is the total heating rate, and are the convective and radiative heating rates, respectively, and are calculated as follows:
[0102]
[0103] where C c represents a constant dependent on the atmospheric composition, and r n represents the radius of the stagnation zone, and C r is a constant dependent on the planetary atmosphere, F(v) is a look-up table function dependent on the flight speed magnitude and atmospheric composition, and the exponents b and c are determined by ρ and v;
[0104] Step S4: Establish a deterministic standard optimization model; in Step S4, construct a deterministic AGA trajectory optimization model:
[0105] Min.J[x,u] = -|ΔE|
[0106]
[0107] H min -H(x(t), C L (t)) ≤ 0.
[0108] Step S5: Establish an optimization model containing uncertainties based on the RS integral theory; in Step S5, the dynamic system with uncertainties is described as:
[0109]
[0110] where p represents the M-dimensional uncertain parameters of the dynamic model, and the objective function is defined as the integral of the original cost function with respect to the uncertain parameter p:
[0111] ∫…∫ -|ΔE(p)| dα(p)
[0112] where ΔE(p) represents the cost function of the energy change, which is a function of the uncertain parameter values, is the known joint probability distribution function CDF of p, the multiplicity of the integral corresponds to the dimension of the uncertainty parameter, and the integration region is over the probability distribution domain of p;
[0113] Inspired by the construction of the objective function, construct the Riemann–Stieltjes optimization RSO model of the dynamic system containing uncertain parameters as follows:
[0114] Min.J[x(·,·), C L (·)] = ∫ dom(p) -|ΔE(p)| dα(p)
[0115]
[0116] ∫ dom(p) (H min -H(x(t), C L (t); p)) dα(p) ≤ 0
[0117] Among them, dom(p) represents the probability distribution domain of p, and ∫ dom(p) is a shorthand symbol for multiple integral.
[0118] Step S6: Establish a standard robust trajectory optimization model based on the orthogonal point selection strategy;
[0119] In step S6, the RSO model obtained in step S5 is discretized. First, the integral of the original cost function of the objective function with respect to the uncertain parameter p is approximated as a finite sum form, and the optimization objective function is obtained as follows:
[0120]
[0121] Among them n represents the value set of the uncertain parameter with a definite probability density function PDF, and W i represents the weight coefficient after integral discretization, is the joint CDF with respect to the uncertain parameter p, and its differential dα(p) is the joint PDF of p. W i is the value of the joint PDF of the corresponding uncertain parameter. The constraint function in the fractional form is discretized into the following form:
[0122]
[0123] In the AGA optimization problem, two independent random variables ρ0 and φ = S / m are considered. The point selection strategy is introduced with a two-dimensional Gaussian distribution. As Figure 3 shown, for the two-dimensional uncertainty parameter p = (ρ0, φ) that satisfies the Gaussian distribution, three points are uniformly selected on each uncertainty parameter axis. With its mean μ = (ρ0, φ) as the center and the standard deviation σ = 3×10 -5 ·(ρ0, φ) as the boundary, the following 5 sampling points are obtained:
[0124]
[0125] The weight coefficient W of the sampling point i is obtained through the joint PDF of the uncertain parameter, as follows:
[0126]
[0127] Among them, Σ represents its covariance matrix, and the specific formula is as follows:
[0128]
[0129] Rewrite the AGA optimization problem model with dynamic uncertainty into the standard form:
[0130]
[0131] Step S7: Solve the optimization model in Step S6 using the pseudospectral method. In Step S7, use the continuous optimal control problem solver GPOPS-II to solve the AGA optimization problem model of dynamic uncertainty obtained in Step S6. The nominal orbit control law of the optimization result is as Figure 4 shown, and the scatter of the orbit altitude and speed under robust optimization is as Figure 5 shown.
[0132] Therefore, the present invention adopts the above-mentioned robust trajectory optimization method for aerodynamic assisted deorbiting in rarefied atmosphere, uses the Riemann–Stieltjes integral to convert the system of stochastic ordinary differential equations into a system of deterministic ordinary differential equations, establishes extended process constraints and cost functions, and converts the optimization model containing system parameter uncertainties into a standard form. And a point selection strategy is proposed to simplify the problem complexity, and finally a robust trajectory optimization model suitable for the pseudospectral method is constructed. Finally, use the GPOPS tool to solve this standard optimization model, and realize the robust trajectory generation under a dynamic system with multi-dimensional parameter uncertainties.
[0133] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that they can still modify or equivalently replace the technical solutions of the present invention, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present invention.
Claims
1. A robust trajectory optimization method for aerodynamic assisted deorbiting in rarefied atmosphere, characterized in that: It includes the following steps: Step S1: Establish a deterministic dynamics model; Step S2: Determine the optimization objective, which is to maximize the energy change through the planetary influence sphere region; Step S3: Establish deterministic constraint functions, including altitude constraint and heating rate constraint: Altitude constraint: The altitude of the spacecraft is not lower than the allowed minimum altitude: H = r - R0 ≥ H min where r is the Cartesian position of the spacecraft in the inertial system centered at the planet, and R0 is the radius of the planet; Heating rate constraint: The total heating rate of the spacecraft does not exceed the allowed maximum value: where, is the total heating rate, and are the convective and radiative heating rates, respectively, and are calculated as follows: where C c represents a constant that depends on the atmospheric composition, r n represents the radius of the stagnation zone, C r is a constant that depends on the planetary atmosphere, F(v) is a look-up table function that depends on the magnitude of the flight speed and the atmospheric composition, and the exponents b and c are determined by ρ and v; Step S4: Establish a deterministic standard optimization model; Step S5: Based on the RS integral theory, establish an optimization model including uncertainties; Step S6: Based on the orthogonal sampling strategy, establish a standard robust trajectory optimization model; Step S7: Use the pseudospectral method to solve the optimization model in Step S6.
2. The robust trajectory optimization method for aeroassisted deorbiting in rarefied atmosphere according to claim 1, wherein: In Step S1, the dynamic equation for planar atmospheric flight without considering thrust and planetary rotation is: where the orbital state \(x = [r, v]\), \(u\) is the control variable, and the lift coefficient \(C\) L is used as the control variable. The specific calculation formula is as follows: where r = (x,y) and v = (v x ,v y ) are the Cartesian position and velocity states of the spacecraft in the inertial frame centered on the planet, and μ represents the planet mass times the gravitational constant. L and D represent the accelerations due to lift and drag, respectively, and are expressed as: where S is the aerodynamic reference area of the spacecraft, m is the mass of the spacecraft, v n = (-v y , v x ) is the normal vector perpendicular to the velocity vector v, v is the magnitude of the velocity ||v||, and σ is the motion direction index of the entry orbit relative to the planet: prograde and retrograde correspond to 1 and -1 respectively; Drag coefficient C D and lift coefficient C L have the following relationship: where: (L / D) max represents the maximum lift-to-drag ratio, which is (L / D) max at the lift coefficient; Regarding atmospheric modeling, the density ρ is approximated as an exponential function within the atmosphere and is defined by the equation ρ = ρ0exp(-βH), where β represents the reciprocal of the standard altitude of the atmosphere, ρ0 is the reference density at sea level, and H is the altitude of the spacecraft.
3. A robust trajectory optimization method for aeroassisted orbit descent in rarefied atmosphere according to claim 2, characterized in that: In Step S2, the energy of the spacecraft in the sun-planet system is defined as: where V is the magnitude of the heliocentric velocity of the spacecraft, μ s is the gravitational constant times the mass of the sun, r s represents the distance between the spacecraft and the sun, and is expressed by the energy change of the sphere of influence of the planet as: ΔE = E + - E - where E + and E - represent the energy values at the outlet and inlet of the planetary influence sphere, respectively.
4. A method for optimizing a robust trajectory of a thin - atmosphere aerodynamic - assisted orbit - descent, as described in claim 3, characterized in that: In Step S4, construct a deterministic AGA trajectory optimization model: Min.J[x,u] = -|ΔE| H min -H(x(t), C L (t)) ≤ 0。 5. A method for optimizing a robust trajectory of a thin atmosphere aerodynamic assisted orbit descent according to claim 5, characterized in that: In Step S5, the dynamic system with uncertainties is described as: where p represents the M-dimensional uncertain parameter of the dynamic model, and the objective function is defined as the integral of the original cost function with respect to the uncertain parameter p: ∫…∫-|ΔE(p)|dα(p) where ΔE(p) represents the cost function of energy change, which is a function of the values of uncertain parameters, is the known joint probability distribution function CDF of p. The multiplicity of the integral corresponds to the dimension of the uncertain parameters, and the integration region is over the probability distribution domain of p; Inspired by the construction of the objective function, construct the Riemann–Stieltjes optimization RSO model of the dynamic system containing uncertain parameters as follows: Min.J[x(·,·),C L (·)] = ∫ dom(p) -|ΔE(p)|dα(p) ∫ dom(p) (H min -H(x(t),C L (t); p))dα(p) ≤ C2 where C1, C2 denote constants, dom(p) represents the domain of the probability distribution of p, and ∫ dom(p) is a shorthand notation for a multiple integral.
6. A robust trajectory optimization method for aeroassisted orbit descent in rarefied atmosphere according to claim 5, characterized in that: In Step S6, discretize the RSO model obtained in Step S5. First, approximate the integral of the original cost function of the objective function with respect to the uncertain parameter p as a finite sum form, and obtain the optimization objective function in the following way: where i = 1, 2, ..., n, where n represents the set of values of the uncertain parameters with a definite probability density function PDF, and W i represents the weight coefficient after product discretization, is the joint CDF with respect to the uncertain parameter p, and its differential dα(p) is the joint PDF of p, and W i is the value of the joint PDF of the uncertain parameter corresponding to The discrete form of the constraint function is as follows: In the AGA optimization problem, two independent random variables ρ0 and φ = S / m are considered, and the point selection strategy is introduced with a two-dimensional Gaussian distribution. For the two-dimensional uncertainty parameter p = (p1, p2) that satisfies the Gaussian distribution, three points are uniformly selected on the axis of the uncertainty parameter p k with the mean value μ k as the center and the standard deviation σ k as the boundary, where k = 1, 2,..., M, M = 2, to obtain 2M + 1 = 5 sampling points. The weight coefficients of the sampling points are obtained through the joint PDF of the uncertainty parameters, as follows: where μ represents the mean of the uncertainty parameter, Σ represents its covariance matrix, and rewrite the AGA optimization problem model of the dynamic uncertainty as a standard form:
7. A lean atmosphere aerodynamic assisted deorbiting robust trajectory optimization method according to claim 6, characterized in that: In Step S7, use the continuous optimal control problem solver GPOPS-II to solve the AGA optimization problem model of the dynamic uncertainty obtained in Step S6.