A Small-Thrust Rendezvous Trajectory Optimization Method for Lunar-Earth Halo Orbit with Zero Initial Value Guess

Abstract

The present invention discloses a method for optimizing the small-thrust rendezvous trajectory of the Earth-Moon Halo orbit with zero initial value guess, belonging to the field of small-thrust rendezvous trajectory optimization, and solving the problems of inability to guarantee the optimal result, long calculation time, poor robustness and initial value sensitivity; including: discretizing the original control trajectory and state trajectory to obtain an objective function that satisfies the rendezvous mission constraints; solving the intermediate variables of the objective function, and calculating the discrete co-state quantity through the co-state mapping equation; combining the Pontryagin minimum principle, obtaining the optimal control of the tracking spacecraft through the Hamiltonian to construct a two-point boundary value model; using the discrete co-state quantity as the initial guess of the two-point boundary value model, and using the SQP algorithm to optimize the discrete co-state quantity to obtain the output result; integrating the forward dynamics equation and the co-state equation through the output result to obtain the optimal control rate of the tracking spacecraft and the small-thrust rendezvous trajectory; the present invention overcomes the initial value sensitivity problem, shortens the calculation time, ensures the optimal result and the robustness of the calculation process.

Description

Technical Field

[0001] The present invention belongs to the technical field of small-thrust rendezvous trajectory optimization for the Earth-Moon Halo orbit, and particularly relates to a method for optimizing the small-thrust rendezvous trajectory of the Earth-Moon Halo orbit without initial value guessing. Background Technique

[0002] With the exploration of space by humans, the Earth-Moon space has attracted extensive attention from scholars. At present, the typical orbits existing in the Earth-Moon space are considered important transfer stations for human exploration of the Moon and entry into deep space. China has launched the Queqiao relay satellite, which operates on the Halo orbit at the Earth-Moon L2 point. This satellite has participated in the exploration of the far side of the Moon and serves as a relay communication satellite to maintain the communication link between the tracking station and the lunar rover. The Lunar GateWay planned to be built by the United States focuses on the 9:2-L2 southern family near-rectilinear halo orbit (NRHO) near the Moon. It aims to build a space station to support future deep-space exploration missions. In order to maintain the long-term operation of spacecraft, it is necessary to develop rendezvous and docking (RVD) technology on typical orbits in the Earth-Moon space. The rendezvous and docking mission in the Earth-Moon space was first proposed in the Apollo program, which considered a linear model of the circular restricted three-body problem. Subsequent studies have used relative dynamics models, such as the elliptical linear relative motion equation (ELERM) and the circular linear relative motion equation (CLERM) constructed in the local vertical-local horizontal (LVLH) coordinate system centered on the Moon. After that, scholars have proposed strategies for rendezvous using impulses and continuous thrusts.

[0003] In recent years, compared with traditional chemical fuel propulsion, small-thrust propulsion technology with higher specific impulse and smaller thrust has been developed to improve the utilization efficiency of propellants. As shown in Table 1 of typical propulsion technologies and performances, examples of small-thrust propulsion include electric propulsion and solar sail propulsion suitable for deep-space missions. Electric propulsion converts electrical energy into kinetic energy through electrothermal, electrostatic, or electromagnetic means. Solar sails use lightweight, highly reflective films to reflect photons and generate thrust from solar radiation pressure. For rendezvous problems, electric propulsion is more suitable for close-range missions, while solar sails have not been successfully applied due to their lower thrust and more complex sail control (including the control of sails and support mechanisms). These new propulsion technologies correspond to small-thrust control strategies, such as linear quadratic regulator (LQR) control.

[0004] Table 1

[0005]

[0006]

[0007] Compared with the optimization problems constituted by traditional impulsive maneuvers, the low-thrust propulsion mode increases the complexity of the optimization problem due to the higher dimension of the optimization domain. Generally speaking, the low-thrust rendezvous trajectory optimization problem is usually solved by direct methods or indirect methods;

[0008] For example, the publication number: CN112084581A, the invention name: A method and system for optimizing the low-thrust perturbed rendezvous trajectory of a spacecraft; this technical solution calculates the low-thrust rendezvous impulse velocity increment according to the given initial orbital position and target orbital position of the spacecraft, and then substitutes it as the initial value and the optimality criterion of the solution into the indirect method optimization model to obtain the optimal low-thrust transfer trajectory; although this method solves to a certain extent the problem of large dynamic recurrence calculations and low calculation efficiency when the transfer time is long in the traditional indirect method low-thrust rendezvous optimization algorithm, it still faces the problems of non-optimal solution set and sensitivity of the co-state quantity at the initial moment.

[0009] Although some techniques have been proposed in the prior art to reduce the sensitivity of the co-state quantity at the initial moment, such as the homotopy method and the continuation method, most researchers prefer a heuristic hybrid method. However, the most obvious problem is its robustness, that is, different solutions will be generated under the same initial condition settings.

[0010] In summary, the purpose of low-thrust rendezvous trajectory optimization is to quickly solve the optimal solution in the given scenario, and obtain the optimal control rate of the tracking spacecraft and the low-thrust rendezvous trajectory of the Earth-Moon Halo orbit. The direct method does not consider the necessary conditions and cannot guarantee the optimality of the result; the indirect method requires a guessed value of the co-state quantity and has the problem of initial value sensitivity; the heuristic hybrid method has a long calculation time and poor robustness. Summary of the Invention

[0011] To solve the above technical problems, the present invention proposes a method for optimizing the low-thrust rendezvous trajectory of the Earth-Moon Halo orbit without initial value guessing. This method is a method for calculating and optimizing the continuous low-thrust spacecraft / aircraft rendezvous trajectory design problem without initial value guessing. Comparing with the defects of traditional methods, the problem of initial value sensitivity will make it difficult to obtain a convergent solution that meets the constraints in the calculation process, resulting in the spacecraft being unable to complete the rendezvous mission. Using traditional heuristic methods requires a large amount of computing resources (from a few minutes to dozens of minutes), and it is impossible to achieve fast calculation, thus missing a good mission window. The non-fuel optimal control result will cause the spacecraft to lose unnecessary fuel in actual control, thereby reducing the on-orbit life of the spacecraft. The non-robust calculation process will make the algorithm only applicable to specific rendezvous scenarios / tasks, and it is difficult to be competent for the calculation of the optimal orbital control rate of the spacecraft in multi-scenarios / tasks. The purpose of the present invention is to overcome the initial value sensitivity problem of the indirect method, shorten the calculation time, ensure the optimality of the result, and the robustness of the calculation process.

[0012] The object of the present invention is specifically achieved through the following technical solutions:

[0013] The present invention discloses a small-thrust rendezvous trajectory optimization method for a lunar Halo orbit with zero initial value guess, and the method includes:

[0014] Step 1: After discretizing the received original small-thrust rendezvous control trajectory and state trajectory of the lunar Halo orbit, determine the control quantity and state quantity at the discrete Legendre-Gauss-Radau collocation points to obtain an objective function that satisfies the rendezvous mission constraints;

[0015] Step 2: Use the open-source nonlinear programming solver IPOPT to solve the intermediate variables of the objective function;

[0016] Step 3: Take the intermediate variables as inputs and calculate the discrete co-state quantities through the co-state mapping equation;

[0017] Step 4: Combining with the Pontryagin minimum principle, obtain the optimal control of the chasing spacecraft through the Hamiltonian, and form a two-point boundary value model from the optimal control of the chasing spacecraft, the dynamic equation of the chasing spacecraft, the co-state equation, and the terminal constraints;

[0018] Step 5: Take the discrete co-state quantities as the initial guesses of the co-state quantities at the initial moment of the two-point boundary value model, use the SQP algorithm to optimize the co-state quantities at the initial moment, obtain the output results of the two-point boundary value model, and obtain the minimum value of the objective function through the output results;

[0019] Step 6: Use the output results to integrate the dynamic equation and co-state equation of the chasing spacecraft forward to obtain the optimal control rate of the chasing spacecraft and the small-thrust rendezvous trajectory of the lunar Halo orbit.

[0020] In Step 1, the Legendre-Gauss-Radau collocation points refer to the roots of the P N (τ) + P N-1 (τ) polynomial;

[0021] In the formula, P N (τ) is the Nth-order Legendre polynomial, and P N-1 (τ) is the (N - 1)th-order Legendre polynomial, where P N (τ) is calculated by the following formula:

[0022]

[0023] Replace N in P N (τ) with N - 1 to obtain the expression of P N-1 (τ), where τ is the independent variable, d is the differential symbol, and N is the specified order.

[0024] In step 1, the objective functions that satisfy the rendezvous mission constraints include: the objective function of time-optimal control and the objective function of fuel-optimal control; where,

[0025] (1) The objective function of time-optimal control is: Among them, time-optimal control refers to satisfying the control rate of the chaser spacecraft in the rendezvous mission, making the time required for rendezvous the shortest;

[0026] (2) The objective function of fuel-optimal control is: Among them, fuel-optimal control refers to satisfying the control rate of the chaser spacecraft in the rendezvous mission, making the fuel required for rendezvous the least;

[0027] In the formula, J1 is the time required for rendezvous, t0 is the initial time, t f is the terminal time, J2 is the fuel required for rendezvous, and U is the control vector of the chaser spacecraft.

[0028] In step 2, the intermediate variables include: Gauss quadrature weights, Lagrange multiplier matrix, and differential matrix.

[0029] In step 3, taking the intermediate variables as inputs, the discrete co-states calculated through the co-state mapping equation include: the discrete co-states of time-optimal control and the discrete co-states of fuel-optimal control; where, the co-state mapping equation is:

[0030]

[0031] In the formula, the discrete co-states of time-optimal control and the discrete co-states of fuel-optimal control calculated through the co-state mapping equation are composed of the discrete co-states at the Legendre-Gauss-Radau collocation points and the discrete co-states at the endpoints ; is the (M + 1)-th column of the transposed matrix D T of the differential matrix D, Λ k represents the k-th row of the Lagrange multiplier matrix Λ, ω k is the value of the Gauss quadrature weight ω k at the k-th Legendre-Gauss-Radau collocation point, k is the number of the collocation point, and M is the total number of collocation points.

[0032] In step 4, combining the Pontryagin minimum principle, the optimal control of the chaser spacecraft obtained through the Hamiltonian includes: the time-optimal control U1 of the chaser spacecraft obtained through the Hamiltonian H1 of time-optimal control and the fuel-optimal control U2 of the chaser spacecraft obtained through the Hamiltonian H2 of fuel-optimal control; where,

[0033] (1) The Hamiltonian H1 of time-optimal control is:

[0034]

[0035] where Λ = [Λ r , Λ v , λ m represents the co-state quantity, Λ r is the position co-state quantity, Λ v is the velocity co-state quantity, and λ m is the mass co-state quantity; v is the velocity of the tracking spacecraft, g(r) is the term related to the position r of the tracking spacecraft in the dynamic equation of the tracking spacecraft, h(v) is the term related to the velocity v of the tracking spacecraft in the dynamic equation of the tracking spacecraft, U is the control vector of the tracking spacecraft, T max is the maximum available thrust, m is the mass of the tracking spacecraft, u is the magnitude of the control vector of the tracking spacecraft, I sp is the specific impulse, and g0 is the acceleration due to gravity at sea level;

[0036] According to the Pontryagin minimum principle, the time-optimal control of the tracking spacecraft is obtained through the Hamiltonian H1 of the time-optimal control as follows: where λ v is the magnitude of the velocity co-state quantity;

[0037] (2) The Hamiltonian H2 of the fuel-optimal control is:

[0038]

[0039] According to the Pontryagin minimum principle, the fuel-optimal control of the tracking spacecraft is obtained through the Hamiltonian H2 of the fuel-optimal control as follows: where u2 ∈ [0, T max , representing the magnitude of the fuel-optimal control vector of the tracking spacecraft.

[0040] In step four, the two-point boundary value model includes: the two-point boundary value model of the time-optimal control and the two-point boundary value model of the fuel-optimal control; where,

[0041] The two-point boundary value model of the time-optimal control consists of the dynamic equation of the tracking spacecraft, the time-optimal control of the tracking spacecraft, the co-state equation of the time-optimal control, and the terminal constraint of the time-optimal control;

[0042] The two-point boundary value model of the fuel-optimal control consists of the dynamic equation of the tracking spacecraft, the fuel-optimal control of the tracking spacecraft, the co-state equation of the fuel-optimal control, and the terminal constraint of the fuel-optimal control; where the dynamic equation of the tracking spacecraft is:

[0043]

[0044] where, To track the derivative of the state variables of the spacecraft, f(X,U) is the dynamic equation that the state variables of the tracking spacecraft obey. X = [r, v] is the state variables of the tracking spacecraft, r is the position of the tracking spacecraft, and v is the velocity of the tracking spacecraft. is the derivative of the position of the tracking spacecraft, is the derivative of the velocity of the tracking spacecraft, is the derivative of the mass of the tracking spacecraft, and u = |U| represents the magnitude of the control vector U of the tracking spacecraft;

[0045] (1) The co-state equation of the time-optimal control is:

[0046]

[0047] where is the derivative of the co-state variable of the time-optimal control, and H1 is the Hamiltonian of the time-optimal control;

[0048] The terminal constraint of the time-optimal control is:

[0049]

[0050] where φ1([X,Λ1] T ,t0,t f ) is the terminal constraint value of the time-optimal control, Λ1 is the co-state variable of the time-optimal control of the tracking spacecraft, t0 is the initial time, and t f is the terminal time, r(t f ) is the position of the tracking spacecraft at the terminal time, r f is the position of the target spacecraft at the terminal time, v(t f ) is the velocity of the tracking spacecraft at the terminal time, v f is the velocity of the target spacecraft at the terminal time, and λ m (t f ) is the co-state variable of the mass of the tracking spacecraft at the terminal time, and H1(t f ) is the Hamiltonian at the terminal time, and 08 is an eight-dimensional zero vector;

[0051] (2) The co-state equation of the fuel-optimal control is:

[0052]

[0053] where is the derivative of the co-state variable of the fuel-optimal control, and H2 is the Hamiltonian of the fuel-optimal control;

[0054] The terminal constraint of the fuel-optimal control is:

[0055]

[0056] where φ2([X,Λ2]T , t0, t f ), Λ2 is the co - state quantity of the fuel - optimal control of the tracking spacecraft, and 07 is a seven - dimensional zero vector.

[0057] In step 5, the SQP algorithm is used to optimize the co - state quantity at the initial moment, and the output result of the two - point boundary - value model is obtained. The minimum value of the objective function obtained from the output result includes: the minimum value of the objective function of the time - optimal control and the minimum value of the objective function of the fuel - optimal control; where,

[0058] (1) The calculation method of the minimum value of the objective function of the time - optimal control includes:

[0059] Use the SQP algorithm to optimize the discrete co - state quantity of the time - optimal control; stop the calculation when the control rate calculated by the two - point boundary - value model of the time - optimal control satisfies the terminal constraint of the time - optimal control, and output the co - state quantity at the initial moment of the time - optimal control and the terminal moment of the time - optimal control. The minimum value of the objective function of the time - optimal control is obtained from the co - state quantity at the initial moment of the time - optimal control and the terminal moment of the time - optimal control;

[0060] (2) The calculation method of the minimum value of the objective function of the fuel - optimal control includes:

[0061] Construct the objective function J of the energy - optimal control that satisfies the rendezvous mission constraints;

[0062] Use the open - source nonlinear programming solver IPOPT to solve the intermediate variables of the objective function of the energy - optimal control;

[0063] Take the intermediate variables as inputs and calculate the discrete co - state quantity of the energy - optimal control through the co - state mapping equation;

[0064] Combined with the Pontryagin minimum principle, obtain the energy - optimal control U3 of the tracking spacecraft through the Hamiltonian H3 of the energy - optimal control;

[0065] Take the discrete co - state quantity of the energy - optimal control as the initial guess of the co - state quantity at the initial moment of the two - point boundary - value model of the fuel - optimal control. Use the SQP algorithm to optimize the co - state quantity at the initial moment. Based on the homotopy iteration process, obtain the Hamiltonian of the fuel - optimal control and the fuel - optimal control of the tracking spacecraft from the Hamiltonian of the energy - optimal control and the energy - optimal control of the tracking spacecraft; stop the calculation when the control rate calculated by the two - point boundary - value model of the fuel - optimal control satisfies the terminal constraint of the fuel - optimal control, and output the co - state quantity at the initial moment of the fuel - optimal control. The minimum value of the objective function of the fuel - optimal control is obtained from the co - state quantity at the initial moment of the fuel - optimal control; where,

[0066] The objective function of the energy - optimal control that satisfies the rendezvous mission constraints is:

[0067] The Hamiltonian H3 of the energy-optimal control is as follows:

[0068] Energy-optimal control of the chaser spacecraft where u3 ∈ [0, T max , representing the magnitude of the energy-optimal control vector of the chaser spacecraft.

[0069] In step five, the iterative calculation method of the homotopy iteration process is as follows:

[0070]

[0071] In the formula, is the homotopy parameter for each step; i = 10, 9, …, 1, corresponding to the homotopy parameter ∈ = 1 to ∈ = 0; where ∈ = 1 corresponds to the Hamiltonian of the energy-optimal control, ∈ = 0 corresponds to the Hamiltonian of the fuel-optimal control, and H4 represents the Hamiltonian in the homotopy iteration process;

[0072] For each iteration step, the initial co-state quantity is optimized by the SQP algorithm, stops after satisfying the terminal constraints of the fuel-optimal control, and enters the next iteration step; when the homotopy parameter ∈ = 0, the calculated Hamiltonian of the energy-optimal control is the same as the Hamiltonian of the fuel-optimal control, and the Hamiltonian of the fuel-optimal control is obtained.

[0073] In step six, the integration process of the forward integration uses a Runge-Kutta 7 / 8-order variable-step integrator for forward integration.

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

[0075] 1. After discretizing the received original small-thrust rendezvous control trajectory and state trajectory of the Earth-Moon Halo orbit, the control quantity and state quantity at the discrete Legendre-Gauss-Radau collocation points are determined, and a technical solution for the objective function that satisfies the rendezvous mission constraints is obtained. The conversion from the original continuous small-thrust rendezvous problem to the nonlinear programming process expands the convergence domain of the original problem and reduces the difficulty of solving the continuous small-thrust rendezvous problem.

[0076] 2. The open-source nonlinear programming solver IPOPT is used to solve the intermediate variables of the objective function, and the intermediate variables are used as inputs. The discrete co-state quantity is calculated through the co-state mapping equation, and the discrete co-state quantity is used as the initial guess of the co-state quantity at the initial moment of the two-point boundary value model; a reasonable initial value guess is provided for the co-state quantity at the initial moment of the two-point boundary value model. Compared with the traditional indirect method of selecting random co-state quantities, the initial value sensitivity problem is overcome, and the calculation speed is greatly improved compared with the heuristic method.

[0077] 3. The Pontryagin minimum principle is applied to ensure that the calculation result is the optimal solution.

[0078] 4. In the calculation of the homotopy iteration process, homotopy parameters with gradually decreasing modulus and gradually decreasing step size are used, which improves the robustness of the calculation process.

[0079] 5. The optimal control rate of the tracking spacecraft can be obtained by integrating the minimum value of the objective function forward.

[0080] 6. The technical solutions disclosed in the present invention use different convergence tolerances respectively, which greatly improves the calculation speed compared with the traditional indirect method and heuristic algorithm with a single convergence tolerance. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] The present invention will be further described in detail below with reference to the drawings and embodiments.

[0082] Figure 1 It is a schematic diagram of the earth-moon rotating coordinate system provided by an embodiment of the present invention.

[0083] Figure 2 It is a schematic diagram of the rendezvous trajectory on the Halo orbit near L1 provided by Example 1 of the present invention.

[0084] Figure 3 It is a schematic diagram of the rendezvous trajectory on the Halo orbit near L2 provided by Example 1 of the present invention.

[0085] Figure 4 It is a schematic diagram of the rendezvous trajectory on the near-rectilinear halo orbit of the L2 southern family provided by Example 1 of the present invention.

[0086] Figure 5 It is a schematic diagram of the Halo orbit near L2 and the near-rectilinear halo orbit of the L2 southern family provided by Example 3 of the present invention.

[0087] Figure 6 It is a schematic diagram of the running time of three algorithms in the MT problem provided by Example 4 of the present invention.

[0088] Figure 7 It is a schematic diagram of the running time of three algorithms in the ME problem provided by Example 4 of the present invention.

[0089] Figure 8 It is a schematic diagram of the running time of three algorithms in the MF problem provided by Example 4 of the present invention.

[0090] Figure 9 It is a schematic diagram of the fuel consumption of three algorithms for calculating the MT problem provided by Example 4 of the present invention.

[0091] Figure 10 It is a schematic diagram of the fuel consumption of three algorithms for calculating the ME problem provided by Example 4 of the present invention.

[0092] Figure 11 It is a schematic diagram of fuel consumption for calculating the MF problem by three algorithms provided in Example 4 of the present invention. Detailed implementation manners

[0093] An embodiment of the present invention provides a method for optimizing the small-thrust rendezvous trajectory of a lunar Halo orbit without initial value guess, and the method includes:

[0094] Step 1: After discretizing the received original small-thrust rendezvous control trajectory and state trajectory of the lunar Halo orbit, determine the control quantity and state quantity at the discrete Legendre-Gauss-Radau collocation points, and obtain an objective function that satisfies the rendezvous mission constraints;

[0095] Step 2: Use the open-source nonlinear programming solver IPOPT to solve the intermediate variables of the objective function;

[0096] Among them, the open-source nonlinear programming solver IPOPT can be downloaded from the corresponding website, and the preferred download address is: https: / / coin-or.github.io / Ipopt / .

[0097] Step 3: Take the intermediate variables as inputs, and calculate the discrete adjoint variables through the adjoint mapping equation;

[0098] Step 4: Combine the Pontryagin minimum principle, obtain the optimal control of the chasing spacecraft through the Hamiltonian, and form a two-point boundary value model from the optimal control of the chasing spacecraft, the dynamic equation of the chasing spacecraft, the adjoint equation, and the terminal constraints;

[0099] The Pontryagin minimum principle is a principle for determining the optimal control rate, and the specific expression is: the value of the optimal control rate should make the system Hamiltonian minimum. Among them, the system Hamiltonian is determined by the dynamic equation and the objective function.

[0100] Step 5: Take the discrete adjoint variables as the initial guess of the adjoint variables at the initial time of the two-point boundary value model, use the SQP algorithm to optimize the adjoint variables at the initial time, obtain the output result of the two-point boundary value model, and obtain the minimum value of the objective function through the output result;

[0101] Step 6: Use the output result to integrate the dynamic equation and adjoint equation of the chasing spacecraft forward to obtain the optimal control rate of the chasing spacecraft and the small-thrust rendezvous trajectory of the lunar Halo orbit.

[0102] In Step 1, the Legendre-Gauss-Radau collocation points refer to the roots of the P N (τ) + P N-1 polynomial;

[0103] In the formula, P N(τ) is the Nth order Legendre polynomial, P N-1 (τ) is the N-1 order Legendre polynomial, where P N (τ) is calculated by the following formula:

[0104]

[0105] P N Replace N in (τ) with N-1, and we get P N-1 (τ), where τ is the independent variable, d is the differential symbol, and N is the specified order.

[0106] In step 1, the objective function that satisfies the rendezvous task constraints includes: the objective function of time optimal control and the objective function of fuel optimal control; wherein,

[0107] (1) The objective function of time optimal control is: Among them, time optimal control refers to satisfying the tracking spacecraft control rate in the rendezvous mission so that the time required for rendezvous is the shortest;

[0108] (2) The objective function of fuel optimal control is: Among them, optimal fuel control refers to satisfying the tracking spacecraft control rate in the rendezvous mission so that the fuel required for rendezvous is minimized;

[0109] Where J1 is the time required for the intersection, t0 is the initial time, t f is the terminal time, J2 is the fuel required for rendezvous, and U is the tracking spacecraft control vector.

[0110] In step 2, the intermediate variables include: Gaussian quadrature weights, Lagrange multiplier matrix and differential matrix.

[0111] In step 3, the intermediate variables are used as inputs, and the discrete co-state quantities calculated by the co-state mapping equation include: the discrete co-state quantity of time optimal control and the discrete co-state quantity of fuel optimal control; wherein the co-state mapping equation is:

[0112]

[0113] In the formula, the discrete co-state quantity of time optimal control and the discrete co-state quantity of fuel optimal control calculated by the co-state mapping equation are given by the discrete co-state quantity at the Legendre-Gauss-Radau collocation point and discrete co-state quantities at the endpoints composition, is the transposed matrix D of the differential matrix D T The M+1th column, Λ k represents the kth row of the Lagrange multiplier matrix Λ, ω k is the Gaussian quadrature weight ωk The value at the k-th Legendre-Gauss-Radau collocation point, where k is the number of the collocation point and M is the total number of collocation points.

[0114] In step 4, combining the Pontryagin's minimum principle, the optimal control of the chaser spacecraft obtained through the Hamiltonian includes: the time-optimal control U1 of the chaser spacecraft obtained through the Hamiltonian H1 of the time-optimal control and the fuel-optimal control U2 of the chaser spacecraft obtained through the Hamiltonian H2 of the fuel-optimal control; where

[0115] (1) The Hamiltonian H1 of the time-optimal control is:

[0116]

[0117] In the formula, Λ = [Λ r , Λ v , λ m represents the co-state variables, Λ r is the position co-state variable, Λ v is the velocity co-state variable, λ m is the mass co-state variable; v is the velocity of the chaser spacecraft, g(r) is the term related to the position r of the chaser spacecraft in the dynamic equation of the chaser spacecraft, h(v) is the term related to the velocity v of the chaser spacecraft in the dynamic equation of the chaser spacecraft, U is the control vector of the chaser spacecraft, T max is the maximum available thrust, m is the mass of the chaser spacecraft, u is the magnitude of the control vector of the chaser spacecraft, I sp is the specific impulse, and g0 is the gravitational acceleration at sea level. r = [r x , r y , r z , v = [v x , v y , v z ;

[0118]

[0119] Among them, and respectively represent the distances between the chaser spacecraft and the Earth and the Moon, and μ = 0.01215 is the normalized value of the lunar mass. Among them, the definition of the Earth-Moon rotating coordinate system is: the x-axis of this coordinate system points from the Earth to the Moon, the z-axis is along the angular momentum direction of the lunar orbit, the y-axis is determined according to the orthogonal rule, and the coordinate system center is the Earth-Moon barycenter position. The schematic diagram of the Earth-Moon rotating coordinate system is as shown in Figure 1 .

[0120] As shown in Figure 1As shown, the motion of the tracking spacecraft in the Earth-Moon space can be approximated as the Circular Restricted Three-Body Problem (CR3BP). The Earth and the Moon move in a circular plane around their common center of mass. The motion of the tracking spacecraft is affected by the combined gravitational forces of the two, and the mass of the tracking spacecraft is negligible.

[0121] According to the Pontryagin minimum principle, the time-optimal control of the tracking spacecraft is obtained through the Hamiltonian H1 of time-optimal control as follows: where λ b is the modulus of the velocity co-state quantity;

[0122] (2) The Hamiltonian H2 of fuel-optimal control is:

[0123]

[0124] According to the Pontryagin minimum principle, the fuel-optimal control of the tracking spacecraft is obtained through the Hamiltonian H2 of fuel-optimal control as follows: where u2 ∈ [0, T max , representing the modulus of the fuel-optimal control vector of the tracking spacecraft.

[0125] In step four, the two-point boundary value model includes: the two-point boundary value model of time-optimal control and the two-point boundary value model of fuel-optimal control. The known parameters of the model are the state quantities at the initial time, the state quantities at the terminal time, and the parameter to be solved is the co-state quantity at the initial time; the state quantities at the initial time are composed of the state quantities of the tracking spacecraft and the target spacecraft at the start time of the rendezvous, and the state quantities at the terminal time are composed of the state quantities of the tracking spacecraft and the target spacecraft at the end time of the rendezvous. Due to the rendezvous constraints, the state quantities of the tracking spacecraft and the target spacecraft are equal at the terminal time. The present invention uses r f = [x, y, z] and to represent the state quantity of the target spacecraft at the terminal time; where,

[0126] The two-point boundary value model of time-optimal control is composed of the dynamic equation of the tracking spacecraft, the time-optimal control of the tracking spacecraft, the co-state equation of time-optimal control, and the terminal constraint of time-optimal control;

[0127] The two-point boundary value model of fuel-optimal control is composed of the dynamic equation of the tracking spacecraft, the fuel-optimal control of the tracking spacecraft, the co-state equation of fuel-optimal control, and the terminal constraint of fuel-optimal control; where, the dynamic equation of the tracking spacecraft is:

[0128]

[0129] where, To track the derivative of the state variables of the spacecraft, f(X,U) is the dynamic equation that the state variables of the tracking spacecraft obey. X = [r,v] are the state variables of the tracking spacecraft, r is the position of the tracking spacecraft, and v is the velocity of the tracking spacecraft. is the derivative of the position of the tracking spacecraft, is the derivative of the velocity of the tracking spacecraft, is the derivative of the mass of the tracking spacecraft, and u = |U| represents the magnitude of the control vector U of the tracking spacecraft;

[0130] (1) The co-state equation of the time-optimal control is:

[0131]

[0132] where, is the derivative of the co-state quantity of the time-optimal control, and H1 is the Hamiltonian of the time-optimal control;

[0133] The terminal constraint of the time-optimal control is:

[0134]

[0135] where, φ1([X,Λ1] T ,t0,t f ) is the terminal constraint value of the time-optimal control, Λ1 is the co-state quantity of the time-optimal control of the tracking spacecraft, t0 is the initial time, t f is the terminal time, r(t f ) is the position of the tracking spacecraft at the terminal time, r f is the position of the target spacecraft at the terminal time, v(t f ) is the velocity of the tracking spacecraft at the terminal time, v f is the velocity of the target spacecraft at the terminal time, λ m (t f ) is the mass co-state quantity of the tracking spacecraft at the terminal time, H1(t f ) is the Hamiltonian at the terminal time, and 08 is an eight-dimensional zero vector;

[0136] (2) The co-state equation of the fuel-optimal control is:

[0137]

[0138] where, is the derivative of the co-state quantity of the fuel-optimal control, and H2 is the Hamiltonian of the fuel-optimal control; X = [r,v] are the state variables of the tracking spacecraft;

[0139] (2) The co-state equation of the fuel-optimal control is:

[0140]

[0141] In the formula, is the derivative of the co-state variable of the fuel optimal control, and H2 is the Hamiltonian of the fuel optimal control;

[0142] The terminal constraint of the fuel optimal control is:

[0143]

[0144] In the formula, φ2([X,Λ2] T ,t0,t f ) is the terminal constraint value of the fuel optimal control, Λ2 is the co-state variable of the fuel optimal control of the tracking spacecraft, and 07 is a seven-dimensional zero vector.

[0145] In step five, the SQP algorithm is used to optimize the co-state variable at the initial time to obtain the output result of the two-point boundary value model. The minimum value of the objective function obtained from the output result includes: the minimum value of the objective function of the time optimal control and the minimum value of the objective function of the fuel optimal control; among them,

[0146] (1) In the time optimal control, the terminal time is free. Therefore, the variables to be solved in the two-point boundary value model of the time optimal control are [Λ, t f , a total of 8 dimensions, where t f is the terminal time; the calculation method of the minimum value of the objective function of the time optimal control includes:

[0147] Use the SQP algorithm to optimize the discrete co-state variable of the time optimal control; stop the calculation when the control rate calculated by the two-point boundary value model of the time optimal control satisfies the terminal constraint of the time optimal control, and output the co-state variable at the initial time and the terminal time of the time optimal control. The minimum value of the objective function of the time optimal control is obtained from the co-state variable at the initial time and the terminal time of the time optimal control;

[0148] (2) In the fuel optimal control, the terminal time is fixed. Therefore, the variable to be solved in the two-point boundary value model of the fuel optimal control is Λ of 7 dimensions. Since it is relatively difficult to directly solve the fuel optimal control, the energy optimal (ME) control is first constructed and solved, and then the fuel optimal control is solved through a homotopy iteration process. The calculation method of the minimum value of the objective function of the fuel optimal control includes:

[0149] Construct the objective function of the energy optimal (ME) control that satisfies the rendezvous mission constraints:

[0150] Use the open-source nonlinear programming solver IPOPT to solve the intermediate variable of the objective function J of the energy optimal control;

[0151] Taking the intermediate variables: Gauss quadrature weights, differential matrices, and Lagrange multipliers as inputs, calculate the discrete co - states of the energy - optimal control through the co - state mapping equation;

[0152] Combined with the Pontryagin minimum principle, through the Hamiltonian H3 of the energy - optimal control: Obtain the energy - optimal control of the chaser spacecraft where, u3 ∈ [0, T max , representing the modulus of the energy - optimal control vector of the chaser spacecraft;

[0153] Taking the discrete co - states of the energy - optimal control as the initial guess of the co - states at the initial time of the two - point boundary - value model of the fuel - optimal control, use the SQP algorithm to optimize the co - states at the initial time. Based on the homotopy iteration process, stop the calculation when the control rate calculated by the two - point boundary - value model of the fuel - optimal control satisfies the terminal constraint of the fuel - optimal control, output the co - states at the initial time of the fuel - optimal control, and obtain the minimum value of the objective function of the fuel - optimal control from the co - states at the initial time of the fuel - optimal control.

[0154] In step five, the iterative calculation method of the homotopy iteration process is:

[0155]

[0156] In the formula, is the homotopy parameter for each step; i = 10, 9, …, 1, corresponding to the homotopy parameter ∈ = 1 to ∈ = 0; where, ∈ = 1 corresponds to the Hamiltonian of the energy - optimal control, ∈ = 0 corresponds to the Hamiltonian of the fuel - optimal control, and H4 represents the Hamiltonian in the homotopy iteration process;

[0157] For each iteration step, optimize the co - states Λ at the initial time by the SQP algorithm, stop after satisfying the terminal constraint of the fuel - optimal control, and enter the next iteration step; when the homotopy parameter ∈ = 0, the calculated Hamiltonian of the energy - optimal control is the same as the Hamiltonian of the fuel - optimal control, and obtain the Hamiltonian of the fuel - optimal control. The corresponding result at this time is the fuel - optimal control, that is, the output result of the two - point boundary - value model of the fuel - optimal control, and this result makes the corresponding objective function of the fuel - optimal control take the minimum value. This iterative process corresponds to the transition from the energy - optimal control to the fuel - optimal control and is called the homotopy iteration process.

[0158] In step six, the output results of the two - point boundary - value model of the time - optimal control are the co - states at the initial time of the time - optimal control and the terminal time of the time - optimal control, and the output result of the two - point boundary - value model of the fuel - optimal control is the co - states at the initial time of the fuel - optimal control; the integration process of the forward integration uses a Runge - Kutta 7 / 8 - order variable - step - size integrator for forward integration.

[0159] To enable those skilled in the art to understand the technical solution of the present invention, five specific examples are provided below to elaborate in detail on the technical solution ZIGA disclosed by the present invention:

[0160] ZIGA is applied to the continuous low-thrust rendezvous problem of the Earth-Moon Halo orbit

[0161] Three typical low-thrust rendezvous problems are solved using the zero-initial-value guessing algorithm, and the results are compared with those of the traditional two-impulse maneuver rendezvous (referred to as the three-body Lambert problem (TBLP)). The initial orbits of the target spacecraft and the chaser spacecraft are selected from three typical Earth-Moon orbits.

[0162] Example 1: Halo orbit near L1

[0163] The Halo orbit around L1 is usually selected for scientific missions. In the rendezvous mission, the target and chaser spacecraft are located on the L1 southern family Halo orbit, with an amplitude of 5000 km along the z-axis and a period of 11.99 days. Assume that the chaser spacecraft is at the 0-degree (r y = 0) phase at the start, while the target is at the 12.5-degree phase. In the ME and MF problems, assume the rendezvous time is 0.25 of a Halo orbit period. In step two of ZIGA, the convergence tolerance is set to 10 -5 for quickly solving the MT and ME problems, and in step five, the tolerance is set to 10 -15 for calculating high-precision results. Additionally, T max = 30 mN, I sp = 3000 s, and m = 1500 kg.

[0164] In Figure 2 、 Figure 3 and Figure 4 , I represents the initial position, F represents the final position, O represents the rendezvous orbit, the subscript C represents the chaser spacecraft, and T represents the target spacecraft. The words in parentheses represent the corresponding rendezvous problems. Since the final time is fixed, the final positions of the ME problem, MF problem, and TBLP coincide. However, the rendezvous position of the MT problem is ahead of their phases on the orbit, which means less time is consumed for the rendezvous.

[0165] Example 2: Halo orbit near L2

[0166] Queqiao Satellite is the first communication relay satellite placed in the Earth-Moon L2-Halo orbit, connecting the Earth and the far side of the Moon. In this example, the L2 northern family Halo orbit with an amplitude of 5000 km along the z-axis and a period of 14.91 days is used. The rendezvous time is 0.25 of an orbit period. Assume that the chaser spacecraft is at the apogee (r y = 0), and the phase of the target spacecraft is 12.5 degrees ahead of the chaser spacecraft.

[0167] To solve the rendezvous problem in Halo orbits, it is usually necessary to select a suitable value for λ(t0) and t f In the MT problem, it is well known that the engine will always use the maximum available thrust. The introduction of this prior information greatly reduces the dimension of the NLP problem, thus promoting its rapid convergence with arbitrary precision. However, this comes at the cost of an increase in CPU running time. In the ME problem, low-order numerical integration methods can be used to calculate the dynamic equations. However, during the RPM calculation process, the thrust curve must be approximated by interpolation, which results in a slightly increased solution time for the ME problem compared to the MT problem. The MF problem is the most difficult to solve and requires determining the appropriate rendezvous time and reasonable initial values. To accelerate the solution of the MF problem, a larger tolerance is used in step two of ZIGA, which enables the rapid acquisition of the co-state mapping results and initiates the subsequent steps of the algorithm. Considering the rapid change of the right side of the dynamic equation near the switching point of the optimal bang-bang control, a high-order variable-step numerical integration method is usually required in step five. Considering the terminal relative distance requirement for rendezvous and docking, the tolerance in step five is usually smaller than that in step two, resulting in a longer running time for the MF problem than for the MT and ME problems. If the rendezvous time is less than the solution of the MT problem, or if random initial-time co-state quantities are adopted, the result is non-convergence.

[0168] Example 3: L2 South Family Near-Rectilinear Halo Orbit

[0169] The near-rectilinear halo orbit (NRHO) is a special subset of the north-south branch halo orbits around the L1 and L2 points. LunarGateWay utilizes the L2 South Family Near-Rectilinear Halo Orbit because of its excellent performance, including low maintenance costs, easy access to the lunar surface and return to Earth, and convenience for round trips to deep space. To asymptotically generate the L2 South Family Near-Rectilinear Halo Orbit from the L2 Lyapunov orbit, another state correction function is constructed in this example, which assumes δx = 0.

[0170]

[0171] Figure 5 Shows the process of constructing the L2 Lyapunov orbit to the L2 South Family Near-Rectilinear Halo Orbit by alternately using these two correction functions. The ratio in parentheses represents the ratio of the orbital period to the lunar rotation period. Table 2 shows the comparison of the small-thrust rendezvous results and the impulsive rendezvous results.

[0172] Table 2

[0173]

[0174]

[0175] The present invention uses the same 9:2 NRHO as Lunar GateWay, with a perilune of 3250 km and an orbital period of 6.562 days. Assume that the chasing spacecraft is at the apolune at the beginning, and the target spacecraft has a phase lead of 12.5 degrees, and the rendezvous is set after 0.25 orbital periods. Table 2 shows that compared with the three-body Lambert problem (TBLP) (also known as the double-pulse rendezvous), using low-thrust propulsion can reduce fuel consumption by 42.36% to 84.62%. Different control strategies of low-thrust propulsion also affect the final mass, and the final mass is usually negatively correlated with the rendezvous time.

[0176] Example 4: Comparison with Other Algorithms

[0177] The low-thrust trajectory optimization problem usually encounters two main challenges: the sensitivity of the co-state variables at the initial time and the optimality of the results. To solve the sensitivity problem, the co-state variables at the initial time are usually estimated first, and then the indirect shooting method (ISM) is applied to obtain the optimal solution. Two common algorithms for estimating the co-state variables at the initial time for solving the low-thrust trajectory optimization problem are: particle swarm optimization (PSO) and ISM based on random initial co-state variables.

[0178] The present invention uses an AMD Ryzen5-5600H CPU@3.3GHz and 16GB RAM for simulation, and conducts it independently 10 times. The CPU running time is calculated using the tic-toc function in MATLAB software, as Figure 6 、 Figure 7 、 Figure 8 shown. The fuel consumption results are as Figure 9 、 Figure 10 、 Figure 11 shown.

[0179] The PSO algorithm consumes a large amount of computing resources and has a low convergence speed, especially in the fuel-optimal two-point boundary value problem. This is due to the continuous generation of new and different information by the particles and the frequent information exchange. In addition, due to the new information generated near the last generation, the PSO algorithm may converge to a suboptimal solution. In contrast, the ISM algorithm produces optimal results by using Newton iteration to optimize along the gradient direction. However, the co-state variables at the initial time are highly sensitive, which may make it difficult for the algorithm to converge and require more trial and error for different random initial guesses. In the most sensitive fuel-optimal two-point boundary value problem, compared with the PSO algorithm, the ISM algorithm may require more computing time.

[0180] The comparison of the three algorithms shows that although the PSO algorithm overcomes the sensitivity problem, it has the problems of long running time and non-optimal solutions. The ISM algorithm performs best in the MT problem, but has convergence problems in the ME and MF problems. The ZIGA method produces optimal results in the MT and MF problems, with a slightly longer running time than ISM in the MT problem and much shorter in the MF problem. In addition, the result of the ME problem is an intermediate step in solving the MF problem and does not affect the final optimality. For the MT problem, ZIGA converges in about 2 seconds, the ME problem converges in 0.5 seconds, and the MF problem converges in 15 seconds. Therefore, ZIGA is suitable for calculating the optimal solution of the small-thrust rendezvous problem on the typical Earth-Moon Halo orbit, overcoming the sensitivity problem of the adjoint variables at the initial moment and having no convergence problem.

[0181] The beneficial effects of the embodiments of the present invention are as follows:

[0182] 1. After discretizing the received original small-thrust rendezvous control trajectory and state trajectory of the Earth-Moon Halo orbit, the control quantities and state quantities at the discrete Legendre-Gauss-Radau collocation points are determined, and the technical solution of the objective function that satisfies the rendezvous mission constraints is obtained. The conversion from the original continuous small-thrust rendezvous problem to the nonlinear programming process expands the convergence domain of the original problem and reduces the difficulty of solving the continuous small-thrust rendezvous problem.

[0183] 2. The open-source nonlinear programming solver IPOPT is used to solve the intermediate variables of the objective function. Taking the intermediate variables as inputs, the discrete adjoint variables are calculated through the adjoint mapping equation, and the discrete adjoint variables are used as the initial guess of the adjoint variables at the initial moment of the two-point boundary value model; it provides a reasonable initial value guess for the adjoint variables at the initial moment of the two-point boundary value model, overcomes the initial value sensitivity problem compared with the traditional indirect method of selecting random adjoint variables, and greatly improves the calculation speed compared with the heuristic method.

[0184] 3. The Pontryagin minimum principle is applied to ensure that the calculation result is the optimal solution.

[0185] 4. In the calculation of the homotopy iteration process, homotopy parameters with gradually decreasing modulus and gradually decreasing step size are used, which improves the robustness of the calculation process.

[0186] 5. The optimal control rate of the tracking spacecraft can be obtained by integrating the minimum value of the objective function forward.

[0187] 6. The technical solutions disclosed in the present invention use different convergence tolerances respectively, and greatly improve the calculation speed compared with the traditional indirect method and heuristic algorithm with a single convergence tolerance.

[0188] As described above, it is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should all be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims described above.

Claims

1. A small-thrust rendezvous trajectory optimization method for the Earth-Moon Halo orbit with zero initial value guess, characterized in that, The method includes: Step 1: After discretizing the received original small-thrust rendezvous control trajectory and state trajectory of the Earth-Moon Halo orbit, determine the control quantities and state quantities at the discrete Legendre-Gauss-Radau collocation points to obtain an objective function that satisfies the rendezvous mission constraints; Step 2: Use the open-source nonlinear programming solver IPOPT to solve the intermediate variables of the objective function; Step 3: Take the intermediate variables as inputs and calculate the discrete co-state quantities through the co-state mapping equation; Step 4: Combining with the Pontryagin minimum principle, obtain the optimal control of the chasing spacecraft through the Hamiltonian. The two-point boundary value model is composed of the optimal control of the chasing spacecraft, the dynamic equation of the chasing spacecraft, the co-state equation, and the terminal constraints; Step 5: Take the discrete co-state quantities as the initial guess of the co-state quantities at the initial time of the two-point boundary value model, use the SQP algorithm to optimize the co-state quantities at the initial time, obtain the output result of the two-point boundary value model, and obtain the minimum value of the objective function through the output result; Step 6: Use the output result to integrate the dynamic equation and co-state equation of the chasing spacecraft forward to obtain the optimal control rate of the chasing spacecraft and the small-thrust rendezvous trajectory of the Earth-Moon Halo orbit; In Step 4, combining with the Pontryagin minimum principle, obtaining the optimal control of the chasing spacecraft through the Hamiltonian includes: obtaining the time-optimal control U1 of the chasing spacecraft through the Hamiltonian H1 of time-optimal control and obtaining the fuel-optimal control U2 of the chasing spacecraft through the Hamiltonian H2 of fuel-optimal control; where (1) The Hamiltonian H1 of time-optimal control is: where Λ = [Λ r , Λ v , λ m represents the costate variables, Λ r is the position costate variable, Λ v is the velocity costate variable, λ m is the mass costate variable; v is the velocity of the chaser spacecraft, g(r) is the term related to the position r of the chaser spacecraft in the dynamic equation of the chaser spacecraft, h(v) is the term related to the velocity v of the chaser spacecraft in the dynamic equation of the chaser spacecraft, U is the control vector of the chaser spacecraft, T max is the maximum available thrust, m is the mass of the chaser spacecraft, u is the magnitude of the control vector of the chaser spacecraft, I sp is the specific impulse, and g0 is the acceleration due to gravity at sea level; According to the Pontryagin minimum principle, the time-optimal control of the tracking spacecraft is obtained through the Hamiltonian H1 of time-optimal control as follows: where λ v is the modulus of the velocity co-state quantity; (2) The Hamiltonian H2 of fuel-optimal control is: According to the Pontryagin minimum principle, the fuel-optimal control of the chaser spacecraft is obtained through the Hamiltonian H2 of the fuel-optimal control as follows: where u2 ∈ [0, T max , representing the magnitude of the fuel-optimal control vector of the chaser spacecraft.

2. The method according to claim 1, characterized in that, In Step 1, the Legendre-Gauss-Radau collocation points refer to the roots of the P N (τ) + P N-1 (τ) polynomials; where, P N (τ) is the Nth-order Legendre polynomial, and P N-1 (τ) is the (N - 1)th-order Legendre polynomial, where P N (τ) is calculated by the following formula: Replace N in P N (τ) with N - 1 to obtain the expression of P N-1 (τ), where τ is the independent variable, d is the differential symbol, and N is the specified order.

3. The method according to claim 1, wherein In Step 1, the objective function that satisfies the rendezvous mission constraints includes: the objective function of time-optimal control and the objective function of fuel-optimal control; where (1) The objective function of time-optimal control is as follows: Among them, time-optimal control means satisfying the control rate of the tracking spacecraft in the rendezvous mission to minimize the time required for rendezvous; (2) The objective function of fuel optimal control is as follows: Among them, fuel optimal control means satisfying the control rate of the chaser spacecraft in the rendezvous mission to minimize the fuel required for rendezvous; Where, J1 is the time required for rendezvous, t0 is the initial time, t f is the terminal time, J2 is the fuel required for rendezvous, and U is the control vector of the chaser spacecraft.

4. The method according to claim 3, wherein In Step 2, the intermediate variables include: Gaussian quadrature weights, Lagrange multiplier matrices, and differential matrices.

5. The method according to claim 4, characterized in that, In Step 3, taking the intermediate variables as inputs, the discrete co-state quantities calculated through the co-state mapping equation include: the discrete co-state quantities of time-optimal control and the discrete co-state quantities of fuel-optimal control; where the co-state mapping equation is: In the formula, the discrete co-state variables of the time-optimal control and the discrete co-state variables of the fuel-optimal control calculated by the co-state mapping equation are composed of the discrete co-state variables at the Legendre-Gauss-Radau collocation points and the discrete co-state variables at the endpoints . is the transpose matrix D T of the differential matrix D, the (M + 1)-th column of Λ k represents the k-th row of the Lagrange multiplier matrix Λ, ω k is the value of the Gaussian quadrature weight ω k at the k-th Legendre-Gauss-Radau collocation point, where k is the number of the collocation points and M is the total number of collocation points.

6. The method according to claim 1, wherein In Step 4, the two-point boundary value model includes: the two-point boundary value model of time-optimal control and the two-point boundary value model of fuel-optimal control; where The two-point boundary value model of time-optimal control is composed of the dynamic equation of the chasing spacecraft, the time-optimal control of the chasing spacecraft, the co-state equation of time-optimal control, and the terminal constraints of time-optimal control; The two-point boundary value model of fuel-optimal control is composed of the dynamic equation of the chasing spacecraft, the fuel-optimal control of the chasing spacecraft, the co-state equation of fuel-optimal control, and the terminal constraints of fuel-optimal control; where the dynamic equation of the chasing spacecraft is: In the formula, is the derivative of the state variables of the tracking spacecraft, f(X,U) is the dynamic equation obeyed by the state variables of the tracking spacecraft, X = [r, v] is the state variables of the tracking spacecraft, r is the position of the tracking spacecraft, and v is the velocity of the tracking spacecraft. is the derivative of the position of the tracking spacecraft. is the derivative of the velocity of the tracking spacecraft. is the derivative of the mass of the tracking spacecraft, and u = |U| represents the modulus of the control vector U of the tracking spacecraft. (1) The co-state equation of time-optimal control is: In the formula, is the derivative of the co-state quantity of the time-optimal control, and H1 is the Hamiltonian of the time-optimal control; The terminal constraint of time-optimal control is: where φ1([X,Λ1] T ,t0,t f ) is the terminal constraint value of the time-optimal control, Λ1 is the co-state quantity of the time-optimal control of the chaser spacecraft, t0 is the initial time, t f is the terminal time, r(t f ) is the position of the chaser spacecraft at the terminal time, r f is the position of the target spacecraft at the terminal time, v(t f ) is the velocity of the chaser spacecraft at the terminal time, v f is the velocity of the target spacecraft at the terminal time, λ m (t f ) is the mass co-state quantity of the chaser spacecraft at the terminal time, H1(t f ) is the Hamiltonian at the terminal time, 08 is the eight-dimensional zero vector; (2) The co-state equation of fuel-optimal control is: In the formula, is the derivative of the co-state quantity of the optimal fuel control, and H2 is the Hamiltonian of the optimal fuel control; The terminal constraint of fuel-optimal control is: where φ2([X,Λ2] T ,t0,t f ) is the terminal constraint value of the optimal fuel control, Λ2 is the co-state quantity of the optimal fuel control of the tracking spacecraft, and 07 is a seven-dimensional zero vector.

7. The method according to claim 6, characterized in that, In step 5, the SQP algorithm is used to optimize the initial co-state variables, and the output results of the two-point boundary value model are obtained. The minimum values of the objective functions obtained from the output results include: the minimum value of the objective function for time-optimal control and the minimum value of the objective function for fuel-optimal control. Among them, (1) The calculation method for the minimum value of the objective function for time-optimal control includes: Using the SQP algorithm to optimize the discrete co-state variables for time-optimal control; stop the calculation when the control rate calculated by the two-point boundary value model for time-optimal control satisfies the terminal constraints for time-optimal control, and output the initial co-state variables at the initial time for time-optimal control and the terminal time for time-optimal control. The minimum value of the objective function for time-optimal control is obtained from the initial co-state variables at the initial time for time-optimal control and the terminal time for time-optimal control; (2) The calculation method for the minimum value of the objective function for fuel-optimal control includes: Construct the objective function J for energy-optimal control that satisfies the rendezvous mission constraints; Use the open-source nonlinear programming solver IPOPT to solve the intermediate variables of the objective function for energy-optimal control; Take the intermediate variables as inputs and calculate the discrete co-state variables for energy-optimal control through the co-state mapping equation; Combined with the Pontryagin minimum principle, obtain the energy-optimal control U3 of the tracking spacecraft through the Hamiltonian H3 of energy-optimal control; Take the discrete co-state variables for energy-optimal control as the initial guess of the initial co-state variables of the two-point boundary value model for fuel-optimal control, use the SQP algorithm to optimize the initial co-state variables, and based on the homotopy iteration process, obtain the Hamiltonian for fuel-optimal control and the fuel-optimal control of the tracking spacecraft from the Hamiltonian of energy-optimal control and the energy-optimal control of the tracking spacecraft; stop the calculation when the control rate calculated by the two-point boundary value model for fuel-optimal control satisfies the terminal constraints for fuel-optimal control, and output the initial co-state variables at the initial time for fuel-optimal control. The minimum value of the objective function for fuel-optimal control is obtained from the initial co-state variables at the initial time for fuel-optimal control. Among them, The objective function of the energy-optimal control that meets the rendezvous mission constraints is as follows: The Hamiltonian H3 for energy-optimal control is as follows: Energy-optimal control of a tracking spacecraft where u3 ∈ [0, T max , representing the magnitude of the energy-optimal control vector of the tracking spacecraft.

8. The method according to claim 7, wherein In step 5, the iterative calculation method for the homotopy iteration process is: wherein, is the homotopy parameter for each step; i = 10, 9, …, 1, corresponding to the homotopy parameter ∈ = 1 to ∈ = 0; wherein, ∈ = 1 corresponds to the Hamiltonian of the energy-optimal control, ∈ = 0 corresponds to the Hamiltonian of the fuel-optimal control, and H4 represents the Hamiltonian in the homotopy iteration process; At each iteration step, the SQP algorithm is used to optimize the initial co-state variables, stop after satisfying the terminal constraints for fuel-optimal control, and enter the next iteration step; when the homotopy parameter ∈ = 0, the Hamiltonian for energy-optimal control calculated is the same as the Hamiltonian for fuel-optimal control, and the Hamiltonian for fuel-optimal control is obtained.

9. The method according to claim 1, wherein In step 6, the integration process of forward integration uses a Runge-Kutta 7 / 8-order variable-step integrator for forward integration.

Images