Improved method for fast optimization of spacecraft low earth continuous thrust orbit
By parameterizing Bezier curves and smooth step functions and employing a hierarchical homotopy strategy, the computational complexity and initial value sensitivity of the near-Earth continuous thrust trajectory optimization problem are resolved, achieving efficient and robust trajectory optimization that approximates the theoretically optimal Bang-Bang control.
Patent Information
- Application Number
- CN202611101876.4
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-07-23
- Publication Date
- 2026-11-06
AI Technical Summary
Existing technologies lack a way to systematically integrate shape primitive parameterization and smooth switching approximation to efficiently and robustly solve the near-Earth continuous thrust trajectory optimization problem. They suffer from high computational complexity, sensitivity to initial conditions, and difficulty in approximating the theoretically optimal Bang-Bang structure.
We use Bezier curves to parameterize thrust control variables and orbital geometry variables, introduce a smooth step function, and solve the nonlinear programming problem through a hierarchical homotopy strategy. We then combine the low-dimensional parameterization of Bezier curves and the smooth step function to approximate the Bang-Bang switch structure.
It significantly reduces the dimensionality of the optimization problem, improves the solution speed, reduces sensitivity to initial conditions, approximates the theoretically optimal control structure, and improves computational efficiency and robustness.
Smart Images
Figure CN123049853A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of spacecraft orbital dynamics and trajectory optimization, and particularly to an improved method for rapid optimization of near-Earth continuous thrust orbits for spacecraft. Background Technology
[0002] For spacecraft employing low-thrust electric propulsion systems, the orbital transfer and maintenance problem can be reduced to an optimal control problem constrained by dynamic equations. Traditional solution methods mainly fall into two categories: indirect methods (transforming the problem into a two-point boundary value problem) and direct methods (discretizing a continuous problem into a nonlinear programming problem). Indirect methods are theoretically elegant and can obtain first-order optimality conditions, but they are extremely sensitive to initial values of costate variables and have a small radius of convergence. Direct methods are more robust, but typically require discretizing the state and control variables at numerous time points, resulting in high dimensionality of optimization variables, heavy computational burden, and the solution may lack a clear optimal control structure.
[0003] In fuel optimization problems, theoretical solutions (such as the minimum principle) often indicate that the thrust amplitude exhibits a "bang-bang" structure, meaning the thrust is either maximum or zero. However, forcing this discontinuous switching in direct methods leads to numerical solution difficulties and instability. Shape element methods (such as using Bezier curves, B-splines, etc.) can significantly reduce the dimensionality of the optimization problem by parameterizing the trajectory or control history with a set of smooth basis functions and a small number of control points. Simultaneously, introducing smooth step approximation functions (such as the Sigmoid function) can approximate the Bang-bang control within the optimization framework, maintaining the differentiability of the problem while guiding the solution towards the theoretically optimal structure. Currently, there is a lack of a complete method that systematically integrates shape element parameterization and smooth switching approximation, and provides efficient and robust solutions for near-Earth orbit perturbation environments, overcoming the problems of high computational complexity, sensitivity to initial conditions, and difficulty in approximating the theoretically optimal Bang-bang structure in near-Earth continuous thrust orbit optimization problems. Summary of the Invention
[0004] Purpose of the invention: The technical problem to be solved by the present invention is to provide an improved method for rapid optimization of near-Earth continuous thrust orbits for spacecraft, addressing the shortcomings of existing technologies.
[0005] To address the aforementioned technical problems, this invention discloses an improved method for rapid optimization of near-Earth continuous thrust orbits for spacecraft, comprising the following steps:
[0006] Step 1: Establish a continuous thrust dynamics model of the spacecraft in near-Earth space;
[0007] Step 2: Use Bezier curves to parameterize the thrust control variables or trajectory geometry variables, transforming the infinite-dimensional optimization problem into a finite-dimensional parameter optimization problem;
[0008] Step 3: Introduce a smooth step function in the control parameterization. By adjusting the steepness parameter, the actual thrust amplitude control can be made to approximate the Bang-Bang switch structure.
[0009] Step 4: With the goal of minimizing fuel consumption or maximizing terminal mass, and using the parameterized dynamic equations, initial and terminal state constraints, and path constraints as conditions, construct a nonlinear programming problem;
[0010] Step 5: Use a hierarchical homotopy strategy to solve the problem. Start with a small κ value to obtain a solution to the smooth problem. Use this solution as an initial value and gradually increase the κ value for continuous optimization. Finally, obtain the trajectory and thrust control parameters (i.e., thrust direction vector and throttling rate) that approximate Bang-Bang control.
[0011] In the geocentric inertial frame of reference, the spacecraft's state vector is: ,in For position vectors, It is a velocity vector. For mass, the kinetic equation is:
[0012]
[0013] in, It is the differential of the spacecraft's position vector. It is the differential of the spacecraft's velocity vector. The derivative of the spacecraft's mass, The gravitational constant of Earth, For J2 perturbation acceleration, This is the engine's maximum thrust. For the purpose of comparison, For standard gravitational acceleration, the control variables include the unit vector of the thrust direction. and throttling rate t represents time.
[0014] The parameterization of thrust control variables or orbital geometry variables in step 2 specifically includes: thrust direction parameterization and throttling rate parameterization.
[0015] The thrust direction parameterization includes: representing the thrust direction unit vector using three independent nth-order Bezier curves. A certain projection:
[0016]
[0017] in, These are the control points, i.e., the optimization variables. For n-degree Bernstein basis functions, the actual direction is... , For normalized time.
[0018] The throttling rate parameterization includes: representing the smooth portion of the throttling rate using an m-order Bezier curve.
[0019]
[0020] in, As control points and their values are guaranteed to be within a certain range by constraints. Inside, These are Bernstein basis functions of degree n.
[0021] The smooth step function mentioned in step 3 takes the form of a Sigmoid function:
[0022]
[0023] in, For switching functions, For steepness parameter, when hour, Or 1.
[0024] Throttling rate modified by smooth step function The actual expression is:
[0025] or .
[0026] Step 4, which describes constructing a nonlinear programming problem, specifically includes:
[0027] The optimization variable is: Bezier control points. The initial estimate of the total time it takes for the aircraft to complete its orbital change. ;
[0028] The objective function is: minimize fuel consumption. or minimize ;
[0029] Among them, final quality From the mass flow differential equation The result is obtained by integrating the product of the two products.
[0030] The constraints include: parameterized dynamic equations, boundary conditions, path constraints, and control constraints.
[0031] Step 5 specifically includes:
[0032] Step 5-1, Initial value generation: Based on the initial scenario constraints, the initial and final position and velocity boundary conditions are extracted using the low-thrust analytical solution based on the average orbital elements. The initial values are generated using a third-order Bezier curve, and higher-order mapping is performed through interpolation to obtain the initial Bezier control points.
[0033] Step 5-2, Solving the hierarchical homotopy:
[0034] Step 5-2 specifically includes:
[0035] Step 5-2-1: Set a small initial value for the steepness parameter. Solve the above optimization problem using a nonlinear programming solver to obtain the solution. ;
[0036] Step 5-2-2, with As the initial value, Increase to Solve the new problem and obtain the solution. ;
[0037] Step 5-2-3: Repeat step 5-2-2 until... Reaching the preset limit The final solution is obtained. Control at this time It should be close to the Bang-Bang format.
[0038] Beneficial effects:
[0039] 1. Computationally efficient: By using low-dimensional parameterization of Bezier curves, the number of problem variables is reduced from hundreds or thousands in traditional direct methods to dozens. Combined with analytical / semi-analytical gradient calculation, the solution speed is greatly improved.
[0040] 2. Optimal Structure: By using a smooth step function and homotopy method, the optimization solution is guided to naturally approximate the theoretically optimal Bang-Bang control structure, thereby obtaining better or closer fuel consumption than the theoretical limit.
[0041] 3. Strong robustness: The convex hull and smoothness of the shape primitives, as well as the homotopy method's solution path "from smooth to sharp", reduce the sensitivity to the initial optimization value and improve the robustness of the algorithm's convergence.
[0042] 4. Convenient constraint handling: The endpoint interpolation properties of Bezier curves facilitate the handling of fixed boundary conditions, and their smoothness also facilitates the handling of path constraints. Attached Figure Description
[0043] Figure 1 This is a flowchart illustrating the overall implementation process of the method of the present invention.
[0044] Figure 2 Flowchart for smoothing step function application and homotopy parameter tuning.
[0045] Figure 3 A flowchart illustrating a complete implementation of task input to result analysis.
[0046] Figure 4 The time history of the baseline operating condition state quantity in the embodiment is shown.
[0047] Figure 5 The time history of the control quantity under the baseline operating condition in the embodiment is shown.
[0048] Figure 6 The example shows the three-dimensional transfer trajectory from 400 km to 800 km.
[0049] Figure 7 The example illustrates the evolution and convergence process of homotopy steepness.
[0050] Figure 8 This example compares the state variables and control variables at different target heights. Detailed Implementation
[0051] The specific scenario of this embodiment is: designing a transfer trajectory for a specific continuous low-thrust spacecraft, so that it can enter the target trajectory from the current trajectory through continuous propulsion, and calculating the change in the magnitude of thrust required during the transfer process.
[0052] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments;
[0053] S1: Establish a continuous thrust dynamics model for near-Earth orbit;
[0054] In the geocentric inertial frame of reference, the spacecraft's state vector is: ,in For position vectors, It is a velocity vector. For mass, the kinetic equation is:
[0055]
[0056] in, It is the spacecraft's position vector. It is the differential of the spacecraft's velocity vector. The derivative of the spacecraft's mass, The gravitational constant of Earth, The perturbation acceleration is J2 (the second spherical harmonic coefficient of Earth's gravitational field). This is the engine's maximum thrust. For the purpose of comparison, For standard gravitational acceleration, the control variables include the unit vector of the thrust direction. and throttling rate t represents time.
[0057] 2. Shape primitive parameterization (taking control parameterization as an example)
[0058] Normalized time As a parameter, where Flight time (can be fixed or optimized).
[0059] Thrust direction parameterization: The angular variation of the thrust direction in spherical coordinates can be represented using an nth-order Bezier curve, or the three-dimensional unit vector can be parameterized and then normalized. For example, the thrust direction unit vector can be represented by three independent nth-order Bezier curves. A certain projection:
[0060]
[0061] in, These are control points (optimization variables). For n-degree Bernstein basis functions, the actual direction is... .
[0062] Throttling rate parameterization: The smooth portion of the throttling rate is represented by an m-order Bezier curve.
[0063]
[0064] in, As control points, and their values are guaranteed to be within a certain range through constraints. Inside;
[0065] 3. Smooth switching approximate embedding
[0066] Define a smooth step function, for example, using the Sigmoid function form:
[0067]
[0068] in, This is a switching function, often related to costates in fuel optimization problems. Here, it can be initially defined as some control function (such as the one described later). as well as ). This is the steepness parameter. When... hour, Or 1 (hard switch). Adjust the throttling rate. The actual expression is modified as follows:
[0069] or
[0070] The latter form allows for switching functions With smooth shape Together, they determine the timing and shape of the switching. In optimization, Starting from a small value (such as 1), at this time Smoothness makes optimization problems easier to converge; then gradually increase... (e.g., up to 50 or greater), making A Bang-Bang structure that approximates 0 or 1.
[0071] 4. Formulation of nonlinear programming problems (i.e., optimization problems)
[0072] The optimization variable is: Bezier control points. Total possible time ;
[0073] The objective function is: minimize fuel consumption. or minimize ;
[0074] Among them, final quality From the mass flow differential equation The result is obtained by integrating the product of the two products.
[0075] The constraints include:
[0076] Dynamic equations (substituted with parameterized control).
[0077] Boundary conditions: Given; The target orbit requirements are met.
[0078] Path constraints: such as height The distance from the Earth's center is obtained by taking the 2-norm of the position vector, which is the magnitude of the vector, and then subtracting the Earth's radius. , can then represent height;
[0079] Control constraints: The parameters have been naturally satisfied or the boundary constraints have been imposed.
[0080] 5. Fast solution process
[0081] a) Initial value generation: Based on the initial scenario constraints, the initial and final position and velocity boundary conditions are extracted using the low-thrust analytical solution based on the average orbital elements. The initial values are generated using a third-order Bezier curve, and higher-order mapping is performed through interpolation to obtain the initial Bezier control points.
[0082] b) Solving for layered homotopy:
[0083] 1. Set a small initial value for the steepness parameter. Solve the above optimization problem using a nonlinear programming solver (such as sequential quadratic programming (SQP), interior-point method, nonlinear least squares, etc.) to obtain the solution. .
[0084] 2. with As the initial value, Increase to Solve the new problem and obtain the solution. .
[0085] 3. Repeat step 2 until... Reaching the preset limit The final solution is obtained. At this time, the throttling rate It should be close to the Bang-Bang format.
[0086] The final solution It refers to the set of optimized variables of the transformed nonlinear programming problem (which can also be understood as the solution to the optimization problem). Each set of solutions can uniquely determine the thrust direction vector and the throttling rate.
[0087] By leveraging the property that the derivative and integral of a Bezier curve are still Bezier curves, and using a fast recursive algorithm for their values (de Casteljau algorithm), state integrals and gradients can be efficiently calculated, accelerating optimization iterations.
[0088] Example:
[0089] Mission: To elevate the spacecraft from a 400km circular orbit to an 800km circular orbit, maintaining the same orbital inclination. Consider central gravity and J2 perturbation. Thrusts. The initial mass is 1000 kg.
[0090] Simulation parameter settings: To ensure the reproducibility of this embodiment, the initial orbital inclination angle is set to 30°, and the right ascension of the ascending node, the argument of perigee, and the initial true perigee angle are all set to 0°; the flight time (tf) for target altitude conditions of 600km, 800km, and 1000km is set to 6d, 12d, and 18d, respectively (the unit d is days). Other parameters are shown in Table 1.
[0091] Table 1. Simulation Parameter Settings
[0092] parameter numerical values illustrate Initial orbit 400km circular orbit, i=30° Ω=ω=ν=0° Target orbit / flight time 600 / 800 / 1000km; 6 / 12 / 18d Circular track, average inclination angle remains constant. Spacecraft and thrusters <![CDATA[m0=1000kg;Tmax=0.3N;Isp=300s]]> Continuous thrust, mass changes with propellant consumption Earth Model <![CDATA[μ=3.986004418e+14m³ / s²;J2=1.08262668e-03]]> Central gravity and J2 perturbation Control Parameterization Direction: 5th-order Bezier; Switching function: 4th-order Bezier RTN direction normalization, Sigmoid output throttling rate Homotopic steepness κ=1, 3, 8, 20, 50 The solution from the previous stage is used as the initial value for the next stage. Integral and Optimization DOP853; SLSQP <![CDATA[Final propagation relative tolerance 2×10⁻¹ 0 >
[0093] The dynamics are integrated using a seven-dimensional state in a geocentric inertial frame, and the thrust direction is expressed using six three-dimensional Bezier control points. For this coplanar orbital ascent embodiment, the normal component is fixed at 0, the tangential component is normalized, and radial control points are used to correct the terminal eccentricity. The throttling rate is obtained from five switching control points via a Sigmoid function. After each homotopy stage, the overall bias of the switching function is recalibrated to eliminate semi-major axis drift caused by steepness changes. Each homotopy stage only accepts control solutions converged by SLSQP and passed physical constraint acceptance. When inter-stage refinement is required, intermediate steepness is automatically inserted between the previous steepness and the current target steepness at the geometric midpoint, with a maximum of three subdivisions, thus ensuring continuous transmission of the control sequence and stable arrival at the final steepness.
[0094] Implementation:
[0095] a) Initial value generation: Based on the initial and final boundary conditions, a third-order Bezier curve function is used to obtain the curve parameters by substituting the boundary conditions. The curve parameters are then mapped to an eighth-order Bezier curve using interpolation to obtain the position curve (that is, the control points of the curve can be obtained).
[0096] b) Layered homotopy solution: Thrust direction vector Parameterization was performed using a 5th-order Bezier curve (6 three-dimensional control points, 18 variables in total). The throttling rate switching logic adopted... ,in Parameterized using a 4th-order Bezier curve (5 variables). Starting from 1, the number is gradually increased to 50 in 5 stages, using the SQP solver;
[0097] Result: The final throttling rate obtained (Throttle rate characterizes the magnitude of thrust.) The history shows a clear two-stage "full thrust - coasting - full thrust" quasi-Bang-Bang structure. Compared to not using... Functions, using only Bezier smoothing to represent the throttling rate Compared to other methods, fuel consumption is reduced by approximately 5%. The total optimization computation time (including the entire homotopy process) is reduced by approximately 60% compared to the traditional high-dimensional direct collocation method.
[0098] The specific simulation results are as follows: all three sets of operating conditions use the same dynamic model, control parameterization method, homotopy sequence, and termination error caliber. The final results are shown in Table 2, with the state variables, control variables, and three-dimensional trajectories shown in [Table 2]. Figures 4 to 8 .
[0099] Table 2. Simulation results for three sets of working conditions
[0100] Operating conditions tf / d Propellant / kg Hotopy calculation / s |Δa| / km ef |Δi average| / (°) 400→600 km 6 37.706 129.5 0.015 5.215e-04 0.00044 400→800 km 12 71.888 304.2 0.031 8.014e-04 0.00098 400→1000 km 18 102.749 313.8 0.041 1.158e-04 0.00007
[0101] The final mass under the baseline operating conditions from 400km to 800km was 928.112kg, propellant consumption was 71.888kg, and thrust duty cycle was 67.996%. The terminal semi-major axis error (Δa) was 0.031km, the terminal eccentricity (ef) was 8.014e-04, the average inclination change of the first and last orbits (|Δi average| / ) was 0.00098°, and the lowest orbital altitude was 385.170km. The total calculation time for five-stage homotopy was 304.2s; this time is the actual measured value of this workstation, reflecting the complete calculation scale from the initial steepness to the final control sequence.
[0102] like Figure 5 As shown, the throttling rate under the baseline operating condition successively exceeds 0.5 around t=4.533d and t=8.378d, forming the first thrust segment, the coasting segment, and the second thrust segment. The thrust magnitude is calculated as T(t)=0.3δ(t)N: the output is approximately 0.3N when the throttling rate is close to 1, and the thrust is shut off when it is close to 0. The optimized RTN (Radial-Tangential-Normal) direction satisfies uT≥0.99999979, |uR|≤6.427e-04, and uN=0.0, that is, the thrust is mainly applied along the track tangentially, and the terminal eccentricity is corrected by a small radial component.
[0103] During control execution, at any given time t, the Bezier basis function is calculated based on the normalized time to obtain the RTN direction unit vector and throttle rate. Then, using the current position and velocity, a transformation matrix from the RTN to the geocentric inertial frame is constructed to obtain the inertial frame thrust vector. The throttle rate serves as the thruster switching and power command, while the inertial frame thrust direction serves as the attitude pointing command. Together, they form the thrust control sequence at that moment and drive the spacecraft to complete the orbital transfer.
[0104] Figure 7 The results show that the throttling rate curve is smooth when κ=1; as κ increases to 3, 8, 20 and 50, the transition zone on both sides of the sliding section gradually narrows, and the final throttling rate is a minimum of 1.621e-08 and a maximum of 1.0, thus approximating the Bang-Bang switch structure. Figure 8 Furthermore, it is shown that at target altitudes of 600km, 800km, and 1000km, the method can satisfy the terminal error constraint within a given flight time.
[0105] This invention provides an improved method for rapid optimization of near-Earth continuous thrust orbits for spacecraft. Many methods and approaches exist for implementing this technical solution; the above description is merely a preferred embodiment of the invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of this invention, and these improvements and modifications should also be considered within the scope of protection of this invention. All components not explicitly stated in this embodiment can be implemented using existing technologies.
Claims
1. An improved method for fast optimization of spacecraft low earth continuous thrust orbit, characterized in that, Includes the following steps: Step 1: Establish a continuous thrust dynamics model of the spacecraft in near-Earth space; Step 2: Use Bezier curves to parameterize the thrust control variables or trajectory geometry variables, transforming the infinite-dimensional optimization problem into a finite-dimensional parameter optimization problem; Step 3: Introduce a smooth step function in the control parameterization. By adjusting the steepness parameter, the actual thrust amplitude control can be made to approximate the Bang-Bang switch structure. Step 4: With the goal of minimizing fuel consumption or maximizing terminal mass, and using the parameterized dynamic equations, initial and terminal state constraints, and path constraints as conditions, construct a nonlinear programming problem; Step 5: Use a hierarchical homotopy strategy to solve the problem. Start with a small κ value to obtain a solution to the smooth problem. Use this solution as an initial value and gradually increase the κ value for continuous optimization. Finally, obtain the trajectory and thrust control parameters that approximate Bang-Bang control, namely the thrust direction vector and throttling rate.
2. The improved method for rapid optimization of near-Earth continuous thrust orbits for spacecraft according to claim 1, characterized in that, In the geocentric inertial frame of reference, the spacecraft's state vector is: ,in For position vectors, It is a velocity vector. For mass, the kinetic equation is: in, It is the differential of the spacecraft's position vector. It is the differential of the spacecraft's velocity vector. The derivative of the spacecraft's mass, The gravitational constant of Earth, For J2 perturbation acceleration, This is the engine's maximum thrust. For the purpose of comparison, For standard gravitational acceleration, the control variables include the unit vector of the thrust direction. and throttling rate t represents time.
3. The improved method for rapid optimization of near-Earth continuous thrust orbits for spacecraft according to claim 1, characterized in that, The parameterization of thrust control variables or orbital geometry variables in step 2 specifically includes: thrust direction parameterization and throttling rate parameterization.
4. The improved method for rapid optimization of near-Earth continuous thrust orbits for spacecraft according to claim 1, characterized in that, The thrust direction parameterization includes representing the thrust direction unit vector using three independent nth-order Bezier curves. A certain projection: in, These are the control points, i.e., the optimization variables. For n-degree Bernstein basis functions, the actual direction is... , For normalized time.
5. The improved method for rapid optimization of a spacecraft's near-Earth continuous thrust trajectory according to claim 4, characterized in that, The throttling rate parameterization includes: representing the smooth portion of the throttling rate using an m-order Bezier curve. in, As control points and their values are guaranteed to be within a certain range by constraints. Inside, These are Bernstein basis functions of degree n.
6. The improved method for rapid optimization of a spacecraft's near-Earth continuous thrust trajectory according to claim 4, characterized in that, The smooth step function mentioned in step 3 takes the form of a Sigmoid function: in, For switching functions, For steepness parameter, when hour, Or 1.
7. The improved method for rapid optimization of a spacecraft's near-Earth continuous thrust trajectory according to claim 6, characterized in that, Throttling rate modified by smooth step function The actual expression is: or .
8. The improved method for rapid optimization of a spacecraft's near-Earth continuous thrust trajectory according to claim 6, characterized in that, Step 4, which describes constructing a nonlinear programming problem, specifically includes: The optimization variable is: Bezier control points. The initial estimate of the total time it takes for the aircraft to complete its orbital change. ; The objective function is: minimize fuel consumption. or minimize ; Among them, final quality From the mass flow differential equation The result is obtained by integrating the product of the two products. The constraints include: parameterized dynamic equations, boundary conditions, path constraints, and control constraints.
9. The improved method for rapid optimization of a spacecraft's near-Earth continuous thrust trajectory according to claim 6, characterized in that, Step 5 specifically includes: Step 5-1, Initial value generation: Based on the initial scenario constraints, the initial and final position and velocity boundary conditions are extracted using the low-thrust analytical solution based on the average orbital elements. The initial values are generated using a third-order Bezier curve, and higher-order mapping is performed through interpolation to obtain the initial Bezier control points. Step 5-2: Solve the layered homotopy problem.
10. An improved method for rapid optimization of a spacecraft's near-Earth continuous thrust trajectory according to claim 9, characterized in that, Step 5-2 specifically includes: Step 5-2-1: Set a small initial value for the steepness parameter. Solve the above optimization problem using a nonlinear programming solver to obtain the solution. ; Step 5-2-2, with As the initial value, Increase to Solve the new problem and obtain the solution. ; Step 5-2-3: Repeat step 5-2-2 until... Reaching the preset limit The final solution is obtained. Control at this time It is close to the Bang-Bang form.