Earth-Moon orbit optimization design method and system based on four-pulse ephemeris model based on manifold splicing

Abstract

The present invention discloses a method and system for optimizing the design of Earth-Moon orbits under a four-pulse ephemeris model based on manifold splicing, and relates to the technical field of designing probe transfers between special Earth-Moon orbits. The technical key points of the present invention include: processing the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system based on invariant manifolds, and obtaining the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively; constructing a four-pulse dual-circle restrictive four-body model; using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, and using a particle swarm optimization algorithm and a sequential quadratic programming algorithm to solve the four-pulse dual-circle restrictive four-body model to obtain the optimal spliced ​​orbit. The present invention uses two optimization models to design low-energy transfer orbits, which can accelerate iterations, improve computational efficiency and algorithm convergence, and realize rapid search for low-energy transfer orbits in specific launch windows, meeting the needs of space stations for Earth-Moon low-energy orbit transfers.

Description

Technical Field

[0001] The present invention relates to the technical field of design for probe transfer between special Earth-Moon orbits, and in particular to a method and system for optimizing the Earth-Moon orbit under a four-pulse ephemeris model based on manifold splicing. Background Art

[0002] The unique position for observing the space environment and the low-energy transit station for interplanetary exploration have made the application value of the libration point orbit in space missions increasingly apparent. The libration point orbit can provide a landing point for the deployment of lunar exploration relay satellites and support the navigation missions of landers and rovers on the far side of the moon. Deploying a manned space station on the libration point orbit facilitates material transportation and refueling, as well as the implementation of deep space exploration missions such as Mars exploration and asteroid exploration. At the same time, the invariant manifold can be used to achieve low-energy transfer of the spacecraft, effectively improving the spacecraft's ability to carry payloads. my country's manned lunar space station is an important part of the mission to achieve manned lunar landing and long-term human presence on the lunar base, following the three-step lunar exploration plan of "orbiting, landing, and returning."

[0003] Most current research on the manifold splicing problem focuses on splicing the invariant manifolds of the Earth-Moon circular restricted three-body model with the invariant manifolds of the Sun-Earth circular restricted three-body model. The transfer problem is based on the Halo orbit between the Sun-Earth libration point and the Earth-Moon libration point. However, traditional manifold splicing methods for designing low-energy transfer orbits are difficult to target under the ephemeris model, are sensitive to the initial values ​​of the transfer orbit in the three-body model, and are difficult to solve at multiple points. Furthermore, optimization and correction under the ephemeris model are slow and difficult to solve. These limitations significantly hinder the diversity of transfer orbits and make transfer orbit design difficult. Therefore, considering the needs of manned and material transport, the transfer problem from low-Earth orbit to the Halo orbit where the space station is located urgently needs to be solved. Summary of the Invention

[0004] In view of this, the present invention proposes a method and system for optimizing the Earth-Moon orbit based on a four-pulse ephemeris model with manifold splicing, in order to solve the technical problems of the traditional method with long transfer time and inability to achieve low-energy transfer of the spacecraft.

[0005] According to one aspect of the present invention, a method for optimizing the Earth-Moon orbit using a four-pulse ephemeris model based on manifold splicing is proposed, comprising the following steps:

[0006] Step 1: Based on the invariant manifold, the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are processed to obtain the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively; the optimal manifold splicing points of the Halo orbits include the optimal position and optimal speed of the spacecraft;

[0007] Step 2: Construct a four-pulse dual-circle restricted four-body model; use the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, and use the particle swarm optimization algorithm and sequential quadratic programming algorithm to solve the four-pulse dual-circle restricted four-body model to obtain the optimal splicing orbit.

[0008] Furthermore, in step 1, the process of obtaining the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system includes:

[0009] Step 11: Calculate and obtain a Halo orbit initial value array for the Sun-Earth system or the Earth-Moon system; the initial value array includes the positions and velocities of multiple discrete points during the spacecraft's flight;

[0010] Step 12: using a numerical method to solve the initial value array of the Halo orbit to obtain an invariant manifold of the Halo orbit; the invariant manifold includes an unstable manifold and a stable manifold;

[0011] Step 1: Use the point selection method to splice the invariant manifold of the Halo orbit and calculate the optimal manifold splicing point with comprehensive position deviation and velocity deviation.

[0012] Furthermore, the specific process of step one includes:

[0013] Under the circular confined three-body dynamics model, the Halo orbit starts on the xz plane, and its initial value is obtained by calculating the Richardson third-order analytical solution;

[0014] The initial value is subjected to a cycle of dynamic integration to obtain the final state; the error between the final state and the target state is calculated. If the error is greater than the given error, the final state is substituted into the differential correction equation to obtain the initial value of the next iteration, that is, the corrected initial value;

[0015] Repeat the above process until the error between the final state and the target state is less than the given error, then terminate the iteration and obtain the final corrected initial value;

[0016] The final corrected initial value is subjected to a period of dynamic integration to obtain the Halo orbit initial value array.

[0017] Furthermore, the differential correction solution equation in step 1 is expressed as:

[0018]

[0019] Where, is the state transfer matrix, which represents the linear estimation relationship between the initial state and the final state of the system; x0 represents the position of the spacecraft on the X axis; represents the speed of the spacecraft on the Y axis; Indicates the velocity of the spacecraft on the X axis at the final state; Indicates the acceleration of the spacecraft on the X axis at the end state; y f Indicates the position of the spacecraft on the Y axis at the final state; Indicates the velocity of the spacecraft on the Y axis at the final state; Indicates the velocity of the spacecraft on the Z axis at the final state; It represents the acceleration of the spacecraft on the Z axis at the final state; δ represents the correction amount, and the variable after it represents the correction amount of each iteration of the variable.

[0020] Furthermore, the specific processes of steps one and two include:

[0021] The state transfer matrix corresponding to the final state obtained by the dynamic integration of the final corrected initial value after one cycle is expressed as a single-value matrix, and the single-value matrix is ​​solved to obtain the eigenvalues ​​and eigenvectors; the eigenvalues ​​include: λ1>1 and

[0022] The initial value point of the integration of the unstable manifold is obtained by the following calculation: Where p i Represents the column data in the Halo orbit initial value array; ε represents the perturbation factor; D u represents the eigenvector corresponding to the eigenvalue λ1;

[0023] The initial value point of the integration of the stable manifold is obtained by the following calculation: Where D s represents the eigenvector corresponding to the eigenvalue λ2;

[0024] After obtaining the initial value points of the integration of the two manifolds, the circular restricted three-body dynamics model is used for recursion to obtain the invariant manifold array of the Halo orbit; the invariant manifold array includes an unstable manifold array and a stable manifold array.

[0025] Furthermore, the specific process of steps 1 and 3 includes:

[0026] Calculate the minimum position deviation Δr between the position in the unstable manifold array and the position in the stable manifold array i ;

[0027] Calculate the minimum velocity deviation Δv between the velocity in the unstable manifold array and the velocity in the stable manifold array i ;

[0028] Minimum position deviation Δr i Deviation from minimum speed Δv i Normalize and perform weighted summation to obtain the comprehensive deviation

[0029] right Sort and select the point corresponding to the smallest comprehensive deviation, which is the optimal manifold splicing point of the comprehensive position deviation and velocity deviation.

[0030] Furthermore, the four-pulse dual-circle restricted four-body model in step 4 includes an upper optimization model and a lower optimization model; wherein the upper optimization model is expressed as:

[0031]

[0032] The lower layer optimization model is expressed as:

[0033]

[0034] In the formula, min means minimization; and They represent the optimal velocities at the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system, respectively; ε s represents the disturbance factor corresponding to the Sun-Earth system, ε m represents the disturbance factor corresponding to the Earth-Moon system; Indicates the preset disturbance factor threshold; and They respectively represent the optimal positions of the optimal manifold splicing points in the Halo orbit corresponding to the Sun-Earth system and the optimal manifold splicing points in the Halo orbit corresponding to the Earth-Moon system; t1 represents the time it takes for the spacecraft in the Halo orbit of the Sun-Earth system to run from the final corrected initial value to the optimal manifold splicing point, and t2 represents the time it takes for the spacecraft in the Halo orbit of the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point.

[0035] Furthermore, in step 2, the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the particle swarm optimization algorithm and the sequential quadratic programming algorithm are used to solve the four-pulse based double-circle restricted four-body model to obtain the optimal splicing orbit, which includes:

[0036] Step 21: Using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, the particle swarm optimization algorithm is used to solve the upper optimization model, and the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are output; wherein, during the iterative optimization process of the particle swarm optimization algorithm, a perturbation factor threshold is set. At the end of each iteration, it is determined whether the locally optimal perturbation factor obtained each time is less than or equal to the perturbation factor threshold. If it is greater, the optimization result is penalized, that is, its performance index is increased;

[0037] Step 22: Take the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system, as well as the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, use the sequential quadratic programming algorithm to solve the lower-level optimization model, and output the time it takes for the spacecraft in the Halo orbits of the Sun-Earth system and the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point.

[0038] Furthermore, the method further includes step three: after obtaining the optimal spliced ​​orbit, correcting the position and state of the data points in the optimal spliced ​​orbit based on the four-pulse ephemeris model to obtain the optimal spliced ​​orbit under the four-pulse ephemeris model.

[0039] According to another aspect of the present invention, a system for optimizing the Earth-Moon orbit using a four-pulse ephemeris model based on manifold splicing is proposed, comprising:

[0040] An optimal manifold splicing point calculation module is configured to process the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system based on invariant manifolds, and obtain the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively; the optimal manifold splicing points of the Halo orbits include the optimal position and optimal velocity of the spacecraft; the specific process includes: calculating and obtaining an array of initial values ​​of the Halo orbits of the Sun-Earth system or the Earth-Moon system; the initial value array includes the positions and velocities of multiple discrete points during the flight of the spacecraft; solving the initial value array of the Halo orbit using a numerical method to obtain the invariant manifold of the Halo orbit; the invariant manifold includes an unstable manifold and a stable manifold; splicing the invariant manifolds of the Halo orbit using a point selection method, and calculating the optimal manifold splicing points of the comprehensive position deviation and velocity deviation;

[0041] The optimal splicing orbit acquisition module is configured to construct a four-pulse dual-circle restricted four-body model; using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, using the particle swarm optimization algorithm and the sequential quadratic programming algorithm to solve the four-pulse dual-circle restricted four-body model to obtain the optimal splicing orbit; specifically comprising: using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, using the particle swarm optimization algorithm to solve the upper optimization model, and outputting the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system; wherein, in the particle swarm optimization algorithm, the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the particle swarm optimization algorithm is used to solve the upper optimization model, and output the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used. During the iterative optimization process of the method, a perturbation factor threshold is set. At the end of each iteration, it is judged whether the local optimal perturbation factor obtained each time is less than or equal to the perturbation factor threshold. If it is greater, the optimization result is penalized, that is, its performance index is increased. The optimal manifold splicing point of the Halo orbit corresponding to the Sun-Earth system and the Earth-Moon system, as well as the optimal perturbation factor of the Halo orbit corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the sequential quadratic programming algorithm is used to solve the lower-level optimization model, and the time it takes for the spacecraft in the Halo orbit of the Sun-Earth system and the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point is output.

[0042] The beneficial technical effects of the present invention are:

[0043] The low-energy transfer orbit design method under a high-precision ephemeris model proposed in this invention uses two optimization models to perform low-energy transfer orbit design, which can accelerate iterations to a certain extent and improve computational efficiency. At the same time, the initial value has a low impact on the optimization results, improving the convergence of the algorithm and enabling rapid search for low-energy transfer orbits within a specific launch window, meeting the needs of space stations for Earth-Moon low-energy orbit transfers. The universal optimization method and model provided by this invention can efficiently design low-energy transfer orbits using manifolds under a three-body model. This method can not only serve as a design method for cargo transfer orbits for future lunar base construction, but also serve as a reference for the deployment of future lunar space stations, promoting development and technological progress in these fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] The present invention can be better understood by referring to the description given below in conjunction with the accompanying drawings, which together with the following detailed description are included in this specification and form a part of this specification, and are used to further illustrate the preferred embodiments of the present invention and explain the principles and advantages of the present invention.

[0045] Figure 1 This is a flow chart of the Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to an embodiment of the present invention.

[0046] Figure 2 Schematic diagram of the transfer track design concept in an embodiment of the present invention.

[0047] Figure 3 Schematic diagram of the periodic Halo orbit generated in an embodiment of the present invention.

[0048] Figure 4 This is an example diagram of a special orbital invariant manifold in an embodiment of the present invention.

[0049] Figure 5 Schematic diagram of error deviation of the optimal manifold splicing point in an embodiment of the present invention.

[0050] Figure 6 This is a flow chart for solving a four-pulse based two-circle restricted four-body model in an embodiment of the present invention.

[0051] Figure 7 Schematic diagram of the transfer orbit from the Sun-Earth Halo orbit to the Earth-Moon Halo orbit in an embodiment of the present invention.

[0052] Figure 8 Schematic diagram of the transfer orbit under the double-circle restricted four-body model based on four pulses in an embodiment of the present invention.

[0053] Figure 9 Schematic diagram of low-energy transfer orbit under the high-precision ephemeris model in an embodiment of the present invention. DETAILED DESCRIPTION

[0054] In order to enable those skilled in the art to better understand the present invention, exemplary embodiments or examples of the present invention will be described below with reference to the accompanying drawings. Obviously, the described embodiments or examples are only some of the embodiments or examples of the present invention, and not all of them. Based on the embodiments or examples of the present invention, all other embodiments or examples obtained by those skilled in the art without creative work should fall within the scope of protection of the present invention.

[0055] To solve the problem of transferring from LEO to Earth-Moon Halo orbit under high-precision ephemeris model, multi-point shooting is usually used to correct the initial transfer orbit converted to high-precision ephemeris model. [1] In order to save energy, the initial transfer orbit usually adopts the design method of invariant manifold [2]However, the simple invariant manifold of the Earth-Moon system is difficult to approach large celestial bodies, and the local optimization design method and multi-point shooting have strict requirements on the initial values, low convergence, and it is difficult to select the number of points for multi-point shooting. Therefore, the present invention proposes a method and system for optimizing the Earth-Moon orbit under a four-pulse ephemeris model based on manifold splicing. The initial values ​​of the transfer orbit are designed by splicing the Halo orbits of the Sun-Earth system and the Earth-Moon system. Then, the proposed optimization model is solved by combining the hybrid optimization method of the particle swarm algorithm and the sequential quadratic programming algorithm (SQP) to further optimize the initial values, and then obtain the transfer orbit under the double-circle restricted four-body model. The optimized orbit has high convergence and high computational efficiency under the multi-point shooting method, and the initial value guess of this method is simple, which can approach the transfer orbit with the minimum fuel. It has potential application value in the future design of the Earth-Moon transfer orbit and the future utilization and development of the moon.

[0056] The embodiment of the present invention provides a method for optimizing the Earth-Moon orbit based on a four-pulse ephemeris model with manifold splicing. Figure 1 As shown, the method includes the following steps:

[0057] Step 1: Based on the invariant manifold, the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are processed to obtain the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively; the optimal manifold splicing points of the Halo orbits include the optimal position and optimal speed of the spacecraft; specifically, the following steps are performed:

[0058] Step 11: Calculate and obtain a Halo orbit initial value array for the Sun-Earth system or the Earth-Moon system; the initial value array includes the positions and velocities of multiple discrete points during the spacecraft's flight;

[0059] Step 12: using a numerical method to solve the initial value array of the Halo orbit to obtain an invariant manifold of the Halo orbit; the invariant manifold includes an unstable manifold and a stable manifold;

[0060] Step 13: Use the point selection method to splice the invariant manifold of the Halo orbit and calculate the optimal manifold splicing point with comprehensive position deviation and velocity deviation;

[0061] Step 2: Construct a four-pulse dual-circle restricted four-body model; using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, solve the four-pulse dual-circle restricted four-body model using the particle swarm optimization algorithm and the sequential quadratic programming algorithm to obtain the optimal splicing orbit; specifically, the following steps are performed:

[0062] Step 21: Using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, the particle swarm optimization algorithm is used to solve the upper optimization model, and the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are output; wherein, during the iterative optimization process of the particle swarm optimization algorithm, a perturbation factor threshold is set. At the end of each iteration, it is determined whether the locally optimal perturbation factor obtained each time is less than or equal to the perturbation factor threshold. If it is greater, the optimization result is penalized, that is, its performance index is increased;

[0063] Step 22: Take the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system, as well as the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, use the sequential quadratic programming algorithm to solve the lower-level optimization model, and output the time it takes for the spacecraft in the Halo orbits of the Sun-Earth system and the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point.

[0064] The method begins with step 1. In step 1, based on the invariant manifold, the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are processed to obtain the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively.

[0065] According to an embodiment of the present invention, Figure 2 As shown, Figure 2 The transfer orbit design concept of the present invention is shown. First, in step 11, the Halo orbit initial value array of the Sun-Earth system or the Earth-Moon system is calculated and obtained. The specific process includes: under the circular restricted three-body dynamics model, the Halo orbit starts from the xz plane, and its initial value is obtained by calculating the Richardson third-order analytical solution; the initial value is subjected to a cycle of dynamic integration to obtain the final state; the error between the final state and the target state is calculated. If the error is greater than the given error, the final state is substituted into the differential correction solution equation to calculate the initial value of the next iteration, that is, the corrected initial value; the above process is repeated until the error between the final state and the target state is less than the given error, and the iteration is terminated to obtain the final corrected initial value; the final corrected initial value is subjected to a cycle of dynamic integration to obtain the Halo orbit initial value array.

[0066] Specifically, in the circular restricted three-body dynamics model coordinate system, the center of mass motion of the spacecraft is:

[0067]

[0068] Where: x, y, z are the three-axis coordinates of the spacecraft in the normalized coordinate system; is the velocity of the spacecraft in the three-axis coordinate system in the normalized coordinate system; is the acceleration of the three-axis coordinates of the spacecraft in the normalized coordinate system; μ is the mass parameter of the system; r1 is the distance from the spacecraft to the earth; r2 is the distance from the spacecraft to the moon.

[0069] The Halo orbit is a libration point orbit that exists in the circular restricted three-body dynamics model. According to the third-order approximate analytical solution given by Richardson, the Halo orbit in the circular restricted three-body dynamics model is determined by three parameters: amplitude parameter, north-south family parameter and phase parameter. The amplitude parameter is generally represented by the normal amplitude, the x-axis amplitude A x and the y-direction amplitude A y Available through A z The phase parameter φ determines the specific phase of the initial point on the Halo orbit. The Richardson third-order analytical solution is:

[0070]

[0071] Due to the symmetry of the three-body problem, Halo orbits are divided into two families, the northern and southern families, corresponding to the cases where n is 1 and 3 in the third-order analytical solution. For the L1 point, n = 1 corresponds to the northern family of Halo orbits; for the L2 point, n = 3 corresponds to the northern family of Halo orbits. The northern family of Halo orbits spends most of their time above the plane of motion of the main celestial body, thus providing good coverage of the northern hemisphere of the main celestial body. Similarly, the southern family of Halo orbits provides good coverage of the southern hemisphere of the main celestial body. Therefore, in space mission orbit design, the appropriate Halo orbit can be selected based on actual conditions, such as if the landing point is in the southern (northern) hemisphere.

[0072] In the rendezvous coordinate system of the circularly restricted three-body dynamics model, the Halo orbit is symmetric about the xz plane. Therefore, when the orbit crosses the xz plane, the velocity is perpendicular to the xz plane, that is, the velocities in both the x and z directions are 0. Based on this characteristic of the Halo orbit, a differential correction solution is given.

[0073] Assume that the Halo orbit starts from the xz plane, and the initial state X0 has the following form:

[0074]

[0075] Take X0 as the initial value for the orbital integration. When the orbit returns to the xz plane again, the integration time is considered to be half a period t f , assuming that the final state after half a period of integration is:

[0076]

[0077] Assumptions is the initial correction value of the initial state quantity X0, after half a cycle t f The change at the end is:

[0078]

[0079] Expanded to:

[0080]

[0081] The integral termination condition is to integrate to the xz surface, which can be written as the expression:

[0082]

[0083] The differential correction solution equation is as follows:

[0084]

[0085] In the formula, Φ is the state transfer matrix, which is a 6×6 square matrix. ij Where i is the corresponding row number of the state transfer matrix, and j is the corresponding column number of the state transfer matrix. is the state transfer matrix, which represents the linear estimation relationship between the initial state and the final state of the system; x0 represents the position of the spacecraft on the X axis; represents the speed of the spacecraft on the Y axis; Indicates the velocity of the spacecraft on the X axis at the final state; Indicates the acceleration of the spacecraft on the X axis at the end state; y f Indicates the position of the spacecraft on the Y axis at the final state; Indicates the velocity of the spacecraft on the Y axis at the final state; Indicates the velocity of the spacecraft on the Z axis at the final state; It represents the acceleration of the spacecraft on the Z axis at the final state; δ represents the correction amount, and the variable after it represents the correction amount of each iteration of the variable.

[0086] In summary, the initial value obtained by Richardson's third-order analytical solution is iterated using differential correction to solve the equation. The iterative process is the Newton iteration process until the iterative variable is less than the error. The specific process is: the initial value obtained by Richardson's third-order analytical solution is the Halo orbit initial value X0. After the dynamic integration to the integration termination condition, the final state is X f ; Calculate the final state X f If the error is greater than the given error, the iteration continues until the final state error is less than the given error, and the final corrected initial value X'0 is output. f |≤ε1, ε1 represents a small value, usually 1×10 -7 ; The constraint variables are and The control variables are x0, z0 and By adjusting x0, z0 and Make and If the value is less than ε2, the orbit is considered to be periodic. ε2 is expressed as a small value, usually 1×10 -7 The initial value of the third-order analytical solution requires the z-axis amplitude parameter of Halo. The z-axis amplitude parameter is related to z0, so fix z0, modify x0 and Figure 3 The following figure shows an example of using differential correction to solve the equation for iteration. The final corrected initial value X'0 is shown in Table 1.

[0087] Table 1 Examples of initial values ​​of Halo orbits corresponding to the Earth-Moon system

[0088]

[0089] Finally, the final corrected initial value is subjected to a period of dynamic integration to obtain the Halo orbit initial value array.

[0090] Then, in steps one and two, the Halo orbit initial value array is solved using a numerical method to obtain the invariant manifold of the Halo orbit; the invariant manifold includes an unstable manifold and a stable manifold. The specific process includes: expressing the state transfer matrix corresponding to the final state obtained by the dynamic integration of the final corrected initial value after one cycle as a single-value matrix, and solving the single-value matrix to obtain the eigenvalues ​​and eigenvectors; the eigenvalues ​​include: λ1>1 and The initial value point of the integration of the unstable manifold is obtained by the following calculation: Where p i represents the column data in the Halo orbit initial value array; ε represents the perturbation factor; D u represents the eigenvector corresponding to the eigenvalue λ1; the initial value point of the integration of the stable manifold is obtained by the following calculation: Where D s Represents the eigenvector corresponding to the eigenvalue λ2; after obtaining the initial value points of the integration of the two manifolds, the circular restricted three-body dynamics model is used for recursion to obtain the invariant manifold array of the Halo orbit; the invariant manifold array includes an unstable manifold array and a stable manifold array.

[0091] When a spacecraft is performing orbit transfer, it needs to apply a velocity pulse to itself to change its velocity state and enter the transfer orbit. The embodiment of the present invention uses an invariant manifold and a splicing design of invariant manifolds to optimize the four-pulse transfer orbit, where the four pulses correspond to a pulse from the Earth parking orbit to the Sun-Earth stable manifold, a pulse from the Sun-Earth stable manifold to the Sun-Earth unstable manifold, a pulse from the Sun-Earth unstable manifold to the Earth-Moon stable manifold, and a pulse from the Earth-Moon stable manifold to the lunar orbit.

[0092] The set of Halo initial value arrays obtained after a period of dynamic integration in each step is represented as P1. The dimension of P1 is 6×N, where N is the number of discrete points for the Halo orbit to be discretely selected. That is, the selection of a Halo orbit in a period is represented by N points. The calculation of the invariant manifold of the Halo orbit requires the initial value of the integral of the invariant manifold, and obtaining the initial value of the integral of the invariant manifold requires the use of a single-valued matrix. The single-valued matrix is ​​defined as: The state transition matrix running along the Halo orbit for one cycle is called the single-valued matrix M, that is:

[0093]

[0094] Where Φ(t0;T) is the state transition matrix after running for one cycle T.

[0095] The next step is to calculate the eigenvalues ​​and eigenvectors of the single-value matrix. The eigenvalues ​​of the single-value matrix have the following form. For the same Halo orbit, different points have the same single-value matrix: λ1>1; λ3=λ4=1; |λ5|=|λ6=1|.

[0096] The calculation of the unstable manifold requires the eigenvector corresponding to the eigenvalue λ1, so any point p in the Halo initial value array obtained in step 1 i , the initial value of the integral of the unstable manifold corresponding to i∈(0,1,2…,N) can be obtained by the following formula:

[0097]

[0098] The calculation of the stable manifold requires the eigenvector corresponding to the eigenvalue λ2, so any point p in the Halo initial value array obtained in step 1 is i The initial value of the integral of the corresponding stable manifold can be obtained by the following formula:

[0099]

[0100] Among them, ε is the perturbation factor, which is usually designed through experience; the ± signs correspond to the two branches of the manifold, one flowing to the large celestial body and the other flowing to the small celestial body.

[0101] After obtaining the initial value point of the integration of the manifold, the invariant manifold array of the Halo orbit can be obtained by recursion according to the circular restricted three-body dynamics model (this process is the existing technology); the dimension of the array is N×6×N1, where the array data of the manifold obtained at any point on the Halo orbit is 6×N1, and N1 is the number of points discretized on each manifold. An example of the invariant manifold of the Halo orbit is as follows Figure 4 shown.

[0102] Then, in steps 1 and 3, the invariant manifolds of the Halo orbit are spliced ​​using the point selection method to calculate the optimal manifold splicing point with comprehensive position deviation and velocity deviation. The specific process includes: calculating the minimum position deviation Δr between the position in the unstable manifold array and the position in the stable manifold array i ; Calculate the minimum speed deviation Δv between the speed in the unstable manifold array and the speed in the stable manifold array i ; For the minimum position deviation Δr i Deviation from minimum speed Δv i Normalize and perform weighted summation to obtain the comprehensive deviation right Sort and select the point corresponding to the smallest comprehensive deviation, which is the optimal manifold splicing point of the comprehensive position deviation and velocity deviation.

[0103] According to steps one and two, the invariant manifolds of the Halo orbit of the Sun-Earth system and the Halo orbit of the Earth-Moon system are generated. The difference between the generation of the Earth-Moon Halo orbit and the Sun-Earth Halo orbit is that the normalized quantities such as distance, time and speed are converted into normalized quantities of the Sun-Earth system, and the mass ratio in dynamics is converted into the mass ratio of the Sun-Earth system to obtain the Sun-Earth Halo orbit. The invariant manifold of the Halo orbit of the Sun-Earth system obtained is spliced ​​with the Halo orbit of the Earth-Moon system. A very important point in the manifold splicing process is to select the splicing point of the two manifolds. Although the manifolds on both sides end at the Poincare section, the Poincare section is a two-dimensional plane under the circular restrictive three-body model of the Sun and the Earth, and the yz plane with x=0 is selected here, the manifolds on both sides of the Poincare section still cannot guarantee that the position deviation and velocity deviation can meet the low-energy orbit transfer requirements. Therefore, a point selection method is given below to ensure that the obtained manifolds can be effectively spliced.

[0104] The Halo orbit has been discretized in the above, and the discrete points can be expressed as:

[0105] X=[X1,X2,...,X N ]

[0106] Perform perturbation integration on each starting discrete point, that is, solve the unstable manifold corresponding to each point. The end of the Poincare section is the integration termination condition, and the state matrix of each manifold on the Poincare section is represented as follows:

[0107] X u =[X u1 ,X u2 ,...,X uN ]

[0108] Here, each column represents a point state on the unstable manifold:

[0109]

[0110] Similarly, for the target libration point orbit, the reverse integration after perturbation can be performed on each target discrete point to solve the stable manifold corresponding to each point. The integration termination condition is also set to the Poincare section, and the state matrix representation of each manifold on the Poincare section is obtained as follows:

[0111] X s =[X s1 ,X s2 ,...,X sN ]

[0112] Here, each column represents a point state on a stable manifold:

[0113]

[0114] The minimum position and velocity deviations are calculated based on the manifold cross-section state of the starting track, and the expressions are as follows:

[0115]

[0116] Where r ui is the position of the end of the unstable manifold at any point on the Poincare section, r sj To stabilize the position of the end of the manifold at any point on the Poincare section, v ui is the velocity of any point on the Poincare section at the end of the unstable manifold, v sj is the velocity of any point on the Poincare section at the end of the stable manifold.

[0117] Normalize the deviations to the maximum value, so that the position and velocity deviations are unified to the range [0,1]:

[0118]

[0119] The position deviation and speed deviation are weighted and summed to obtain the comprehensive deviation to balance the influence of the two deviations on the state at this point. The mathematical expression is as follows:

[0120]

[0121] right Sort and select the point corresponding to the smallest comprehensive deviation, which is the optimal manifold splicing point obtained by combining position deviation and velocity deviation. The error deviation example is as follows Figure 5 shown.

[0122] Then, step 2 is executed. In step 2, a four-pulse dual-circle restricted four-body model is constructed. The optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values. The particle swarm optimization algorithm and sequential quadratic programming algorithm are used to solve the four-pulse dual-circle restricted four-body model to obtain the optimal splicing orbit. Specifically, the following steps are performed:

[0123] Step 21: Using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, the particle swarm optimization algorithm is used to solve the upper optimization model, and the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are output; wherein, during the iterative optimization process of the particle swarm optimization algorithm, a perturbation factor threshold is set. At the end of each iteration, it is determined whether the locally optimal perturbation factor obtained each time is less than or equal to the perturbation factor threshold. If it is greater, the optimization result is penalized, that is, its performance index is increased;

[0124] Step 22: Take the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system, as well as the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, use the sequential quadratic programming algorithm to solve the lower-level optimization model, and output the time it takes for the spacecraft in the Halo orbits of the Sun-Earth system and the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point.

[0125] According to an embodiment of the present invention, the optimal manifold splicing point of the comprehensive position deviation and velocity deviation is obtained by step 1, and the time t0 of the Earth-Moon Halo orbit and the Sun-Earth Halo orbit from the starting point to the corresponding optimal splicing point can be obtained. 1 and t0 2 The result and the optimal manifold splicing point are used as the initial optimization input. The upper layer uses the particle swarm optimization algorithm for global optimization, and the optimization variables are the perturbation factors corresponding to the Earth-Moon manifold and the Sun-Earth manifold. The lower layer uses the sequential quadratic programming algorithm (SQP algorithm) for local optimization, and the perturbation factor output by the upper layer is used as one of the inputs of the lower model.

[0126] The specific process is: select the initial value of the differentially corrected Halo orbit as the starting point in the Sun-Earth Halo orbit and the Earth-Moon Halo orbit respectively. and The final state point reached by the Halo orbit after running t1 and t2 times from the starting point is and and is the starting point of the Sun-Earth manifold and the Earth-Moon manifold. The position of the Sun-Earth manifold and the Earth-Moon manifold when they reach the Poincare section is and Speed ​​status is and The manifold of the Earth-Moon system determines the conversion matrix C(x) of the Earth-Moon three-body system to the Sun-Earth three-body system according to the given epoch starting time. The optimization variables of this problem are the time t1 of running in the Sun-Earth Halo orbit and the time t2 of running in the Earth-Moon Halo orbit. The optimization goal is to minimize the velocity error Δv3 on the Poincare section. The constraints of the optimization model are the bounded constraints of no position error and deviation at the Poincare section. The optimal splicing orbit can be obtained by solving this optimization problem. The solution result of the lower model is affected by the output result of the upper model. When the solution value of the manifold deviation is unreasonable, the result of the lower model does not converge. A large index value is assigned to it as a penalty, and the next iterative solution is entered at the same time. The flowchart of the optimization model is as follows: Figure 6 shown.

[0127] The mathematical description of the optimization model (i.e., the four-pulse two-circle restricted four-body model) is as follows:

[0128] The upper optimization model is:

[0129]

[0130] The lower optimization model is:

[0131]

[0132] In the formula, min means minimization; and They represent the optimal velocities at the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system, respectively; ε s represents the disturbance factor corresponding to the Sun-Earth system, ε m represents the disturbance factor corresponding to the Earth-Moon system; Indicates the preset disturbance factor threshold; and They respectively represent the optimal positions of the optimal manifold splicing points in the Halo orbit corresponding to the Sun-Earth system and the optimal manifold splicing points in the Halo orbit corresponding to the Earth-Moon system; t1 represents the time it takes for the spacecraft in the Halo orbit of the Sun-Earth system to run from the final corrected initial value to the optimal manifold splicing point, and t2 represents the time it takes for the spacecraft in the Halo orbit of the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point.

[0133] Among them, the position and velocity update formula of the particle swarm optimization algorithm is:

[0134]

[0135] Where: is the velocity component of the particle at the kth iteration; is the position component of the particle at the kth iteration; ω is the inertia factor; c1 and c2 are acceleration factors; is the individual extreme value of particle i after the first k iterations; is the population extreme value of the entire particle population after the first k iterations; R is a random number. The acceleration factor update strategy is as follows: three different autonomous populations are designed, and the acceleration factor update strategies of the three populations are shown in Table 2.

[0136] Update strategy for Table 2c1 and c2

[0137]

[0138] Select the number of particles, the maximum number of iterations, the maximum inertia weight coefficient, the minimum inertia weight coefficient, The double-layer optimization algorithm is used to optimize the transfer orbit from the Sun-Earth Halo orbit to the Earth-Moon Halo orbit, as shown in Figure 7 As shown in the figure, the result of the transfer orbit is a 6×N array.

[0139] The point on the Sun-Earth Halo orbit corresponding to the initial point of the stable manifold of the Sun-Earth system in the obtained optimal spliced ​​orbit is taken as the starting point, and the unstable manifold of this point is calculated. The transfer orbit from the Sun-Earth Halo orbit to the LEO orbit is obtained by differential correction to satisfy the parking orbit height constraint. The obtained optimal spliced ​​orbit is spliced ​​with the corrected unstable manifold, as shown in Figure 8 The obtained transfer orbits are all 7×N arrays, where the 7th dimension is the time corresponding to the transfer, and the first 6 dimensions are the position and velocity.

[0140] Furthermore, the method further includes step three: after obtaining the optimal spliced ​​orbit, correcting the position and state of the data points in the optimal spliced ​​orbit based on the four-pulse ephemeris model to obtain the optimal spliced ​​orbit under the four-pulse ephemeris model.

[0141] According to an embodiment of the present invention, the optimized initial value is a 7×N array of transfer orbits for the dual-circular restrictive four-body system obtained in step 2. An initial epoch is selected, the transfer orbits are transformed to the Earth-centered inertial coordinate system, the time is converted to the corresponding ephemeris time, and the array of transfer orbits is discretized. For example, if the number of target points is set to k = 20, N is divided into 20 equal parts, and target points are selected one by one. The position, velocity, and ephemeris time corresponding to each target point can be obtained.

[0142] The specific process is as follows:

[0143] (1) Continuous position correction: The state quantity of k points given by the current iteration value and the integration time of each segment are used for orbit recursion. At this time, the speed of the starting point of this segment is changed according to the position correction equation to make the end point of this segment coincide with the position of the starting point of the next segment. The control variable is the speed component of the starting point of this segment. The constraint variable is the position r of the end point if =[x if ,y if ,z if ] T , the termination condition is that the integration time is fixed. The position correction equation can be obtained:

[0144]

[0145] Repeat this correction operation for each segment of the trajectory until all segments of the trajectory are positioned continuously.

[0146] (2) Speed ​​Continuity Correction: After the first level differential correction, i.e. position continuous correction, the trajectory is continuous in position, but there will be a speed jump at each segment connection point (existing at every maneuvering point). In the second level differential correction, the speed will be corrected to make the trajectory continuous in speed. The corresponding speed differential correction relationship is:

[0147]

[0148] After completing the continuous velocity correction, since the position of each segment end point is also changed, continuous position correction is required. This process is repeated until the position error and velocity error are within the allowable range. It can be considered that a transfer trajectory with continuous position and velocity is obtained.

[0149] Because the transfer trajectory needs to meet engineering constraints and the starting and ending points need to be fixed during the correction process, the calculation efficiency of the transfer trajectory is slow and the convergence speed decreases when the number of correction iterations is large. Therefore, a penalty function is set for the velocity termination condition of the corrected trajectory to obtain an inaccurate convergence solution and then substitute it into the optimization model using the sequential quadratic programming algorithm for solution. The optimization variables of the optimization model are the velocity v4 required to leave the starting trajectory and the velocity v5 required to enter the target trajectory. The optimization target is the sum of the squares of the velocity increments required to leave the starting trajectory and enter the target trajectory. The constraint variables are the constraints of the trajectory dynamics and the position constraints of the starting and target points. The mathematical description of the sequential quadratic programming algorithm is as follows:

[0150]

[0151] Where f(x) is the dynamics of the orbit, x c is the starting point of the starting track, x tThe gravitational bodies added to the ephemeris model are: Earth, Sun, Moon, Venus, Mars, and Saturn. By comprehensively correcting the three-segment transfer orbit, namely the stable manifold and unstable manifold of the Sun and Earth, and the stable manifold of the Earth and Moon, the sequential quadratic programming algorithm is used to solve the optimization problem. The optimization variables are the three-dimensional components of the three-segment orbit insertion pulses. The low-energy transfer orbit under the high-precision ephemeris model is obtained, that is, the Earth-Moon orbit under the four-pulse ephemeris model based on manifold splicing. The orbital data is a 6×N3 array, as shown in Figure 9 shown.

[0152] Another embodiment of the present invention provides a system for optimizing the Earth-Moon orbit using a four-pulse ephemeris model based on manifold splicing, the system comprising:

[0153] An optimal manifold splicing point calculation module is configured to process the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system based on invariant manifolds, and obtain the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively; the optimal manifold splicing points of the Halo orbits include the optimal position and optimal velocity of the spacecraft; the specific process includes: calculating and obtaining an array of initial values ​​of the Halo orbits of the Sun-Earth system or the Earth-Moon system; the initial value array includes the positions and velocities of multiple discrete points during the flight of the spacecraft; solving the initial value array of the Halo orbit using a numerical method to obtain the invariant manifold of the Halo orbit; the invariant manifold includes an unstable manifold and a stable manifold; splicing the invariant manifolds of the Halo orbit using a point selection method, and calculating the optimal manifold splicing points of the comprehensive position deviation and velocity deviation;

[0154] The optimal splicing orbit acquisition module is configured to construct a four-pulse dual-circle restricted four-body model; using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, using the particle swarm optimization algorithm and the sequential quadratic programming algorithm to solve the four-pulse dual-circle restricted four-body model to obtain the optimal splicing orbit; specifically comprising: using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, using the particle swarm optimization algorithm to solve the upper optimization model, and outputting the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system; wherein, in the particle swarm optimization algorithm, the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the particle swarm optimization algorithm is used to solve the upper optimization model, and output the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used. During the iterative optimization process of the method, a perturbation factor threshold is set. At the end of each iteration, it is judged whether the local optimal perturbation factor obtained each time is less than or equal to the perturbation factor threshold. If it is greater, the optimization result is penalized, that is, its performance index is increased. The optimal manifold splicing point of the Halo orbit corresponding to the Sun-Earth system and the Earth-Moon system, as well as the optimal perturbation factor of the Halo orbit corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the sequential quadratic programming algorithm is used to solve the lower-level optimization model, and the time it takes for the spacecraft in the Halo orbit of the Sun-Earth system and the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point is output.

[0155] The functions of the Earth-Moon orbit optimization design system under the four-pulse ephemeris model based on manifold splicing in an embodiment of the present invention can be described by the aforementioned Earth-Moon orbit optimization design method under the four-pulse ephemeris model based on manifold splicing. Therefore, for the parts not described in detail in the system embodiment, please refer to the above method embodiment and will not be repeated here.

[0156] Although the present invention has been described with respect to a limited number of embodiments, those skilled in the art, having benefit of the foregoing description, will appreciate that other embodiments are contemplated within the scope of the invention thus described. This disclosure is intended to be illustrative rather than restrictive of the scope of the invention, which is defined by the appended claims.

[0157] The documents cited in the present invention are as follows:

[0158] [1] Lei Hanlun. Libration points, invariant manifolds and low-energy orbits[D]. Nanjing: Nanjing University, 2015.

[0159] [2] Li Mingtao. Design and optimization of energy-saving orbits for collinear libration missions[D]. Beijing: Chinese Academy of Sciences, 2010.

Claims

1. A method for optimizing the Earth-Moon orbit based on a four-pulse ephemeris model with manifold splicing, characterized by: The following steps are involved: Step 1: Based on the invariant manifold, the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are processed to obtain the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively; the optimal manifold splicing points of the Halo orbits include the optimal position and optimal speed of the spacecraft; Step 2: Construct a four-pulse dual-circle restricted four-body model; using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, use the particle swarm optimization algorithm and the sequential quadratic programming algorithm to solve the four-pulse dual-circle restricted four-body model to obtain the optimal splicing orbit; Step 3: After obtaining the optimal spliced ​​orbit, the positions and states of the data points in the optimal spliced ​​orbit are corrected based on the four-pulse ephemeris model to obtain the optimal spliced ​​orbit under the four-pulse ephemeris model; Select an initial epoch, transform the transfer orbit into the geocentric inertial coordinate system, convert the time into the corresponding ephemeris time, discretize the array of the transfer orbit, select the number of target points, select the target points in turn, and obtain the position, velocity and ephemeris time corresponding to each target point; Step three specifically includes: (1) Continuous position correction: The trajectory is recursively calculated using the state quantities of the k points given by the current iteration value and the integration time of each segment. At this time, the velocity of the starting point of this segment is changed according to the position correction equation to make the end point of this segment coincide with the position of the starting point of the next segment; The control variable is the velocity component at the starting point of this segment The constraint variables are the positions of the end points The termination condition is that the integration time is fixed, and the position correction equation is obtained: Repeat this correction operation for each segment of the trajectory until all segments of the trajectory are positioned continuously; (2) Speed ​​Continuity Correction: After the first-level differential correction, i.e., position continuous correction, the trajectory is continuous in position, but there will be a speed jump at each segment connection point; in the second-level differential correction, the speed will be corrected to make the trajectory continuous in speed. The corresponding speed differential correction relationship is: δv=Mδr After completing the continuous velocity correction, since the position of each segment end point is also changed, continuous position correction is required. This process is repeated until the position error and velocity error are within the allowable range, and a transfer trajectory with continuous position and velocity is obtained.

2. The Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to claim 1 is characterized in that: In step 1, the process of obtaining the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system includes: Step 11: Calculate and obtain a Halo orbit initial value array for the Sun-Earth system or the Earth-Moon system; the initial value array includes the positions and velocities of multiple discrete points during the spacecraft's flight; Step 12: using a numerical method to solve the initial value array of the Halo orbit to obtain an invariant manifold of the Halo orbit; the invariant manifold includes an unstable manifold and a stable manifold; Step 1: Use the point selection method to splice the invariant manifold of the Halo orbit and calculate the optimal manifold splicing point with comprehensive position deviation and velocity deviation.

3. The Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to claim 2 is characterized in that: The specific process of step one includes: Under the circular confined three-body dynamics model, the Halo orbit starts on the xz plane, and its initial value is obtained by calculating the Richardson third-order analytical solution; The initial value is subjected to a cycle of dynamic integration to obtain the final state; the error between the final state and the target state is calculated. If the error is greater than the given error, the final state is substituted into the differential correction equation to obtain the initial value of the next iteration, that is, the corrected initial value; Repeat the above process until the error between the final state and the target state is less than the given error, then terminate the iteration and obtain the final corrected initial value; The final corrected initial value is subjected to a period of dynamic integration to obtain the Halo orbit initial value array.

4. The Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to claim 3 is characterized in that: The differential correction solution equation in step 1 is expressed as: Where, is the state transfer matrix, which represents the linear estimation relationship between the initial state and the final state of the system; x0 represents the position of the spacecraft on the X axis; represents the spacecraft's velocity on the Y axis; Indicates the velocity of the spacecraft on the X axis at the final state; Indicates the acceleration of the spacecraft on the X axis at the end state; y f Indicates the position of the spacecraft on the Y axis at the final state; Indicates the velocity of the spacecraft on the Y axis at the final state; Indicates the velocity of the spacecraft on the Z axis at the final state; Indicates the acceleration of the spacecraft on the Z axis at the final state; δ represents the correction amount, and the variable after it represents the correction amount of the variable in each iteration.

5. The Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to claim 4 is characterized in that: The specific process of steps one and two includes: The state transfer matrix corresponding to the final state obtained by the dynamic integration of the final corrected initial value after one cycle is expressed as a single-value matrix, and the single-value matrix is ​​solved to obtain the eigenvalues ​​and eigenvectors; the eigenvalues ​​include: λ1>1 and The initial value point of the integration of the unstable manifold is obtained by the following calculation: Where p i represents the column data in the Halo orbit initial value array; ε represents the perturbation factor; D u represents the eigenvector corresponding to the eigenvalue λ1; The initial value point of the integration of the stable manifold is obtained by the following calculation: Where D s represents the eigenvector corresponding to the eigenvalue λ2; After obtaining the initial value points of the integration of the two manifolds, the circular restricted three-body dynamics model is used for recursion to obtain the invariant manifold array of the Halo orbit; the invariant manifold array includes an unstable manifold array and a stable manifold array.

6. The Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to claim 5 is characterized in that: The specific process of steps 1 and 3 includes: Calculate the minimum position deviation Δr between the position in the unstable manifold array and the position in the stable manifold array i ; Calculate the minimum velocity deviation Δv between the velocity in the unstable manifold array and the velocity in the stable manifold array i ; Minimum position deviation Δr i Deviation from minimum speed Δv i Normalize and perform weighted summation to obtain the comprehensive deviation right Sort and select the point corresponding to the smallest comprehensive deviation, which is the optimal manifold splicing point of the comprehensive position deviation and velocity deviation.

7. The Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to claim 6 is characterized in that: The four-pulse dual-circle restricted four-body model in step 2 includes an upper optimization model and a lower optimization model; wherein the upper optimization model is expressed as: The lower layer optimization model is expressed as: In the formula, min means minimization; and They represent the optimal velocities at the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system, respectively; ε s represents the disturbance factor corresponding to the Sun-Earth system, ε m represents the disturbance factor corresponding to the Earth-Moon system; Indicates the preset disturbance factor threshold; and They respectively represent the optimal positions of the optimal manifold splicing points in the Halo orbit corresponding to the Sun-Earth system and the optimal manifold splicing points in the Halo orbit corresponding to the Earth-Moon system; t1 represents the time it takes for the spacecraft in the Halo orbit of the Sun-Earth system to run from the final corrected initial value to the optimal manifold splicing point, and t2 represents the time it takes for the spacecraft in the Halo orbit of the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point.

8. The Earth-Moon orbit optimization design method based on the four-pulse ephemeris model based on manifold splicing according to claim 7 is characterized in that: In step 2, the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the particle swarm optimization algorithm and sequential quadratic programming algorithm are used to solve the four-pulse based double-circle restricted four-body model to obtain the optimal splicing orbit. Step 21: Using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, the particle swarm optimization algorithm is used to solve the upper optimization model, and the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are output; wherein, during the iterative optimization process of the particle swarm optimization algorithm, a perturbation factor threshold is set. At the end of each iteration, it is determined whether the locally optimal perturbation factor obtained each time is less than or equal to the perturbation factor threshold. If it is greater, the optimization result is penalized, that is, its performance index is increased; Step 22: Take the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system, as well as the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, use the sequential quadratic programming algorithm to solve the lower-level optimization model, and output the time it takes for the spacecraft in the Halo orbits of the Sun-Earth system and the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point.

9. A system for optimizing the Earth-Moon orbit using a four-pulse ephemeris model based on manifold splicing, for implementing the method for optimizing the Earth-Moon orbit using a four-pulse ephemeris model based on manifold splicing according to any one of claims 1 to 8, characterized in that: include: An optimal manifold splicing point calculation module is configured to process the initial values ​​of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system based on invariant manifolds, and obtain the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system respectively; the optimal manifold splicing points of the Halo orbits include the optimal position and optimal velocity of the spacecraft; the specific process includes: calculating and obtaining an array of initial values ​​of the Halo orbits of the Sun-Earth system or the Earth-Moon system; the initial value array includes the positions and velocities of multiple discrete points during the flight of the spacecraft; solving the initial value array of the Halo orbit using a numerical method to obtain the invariant manifold of the Halo orbit; the invariant manifold includes an unstable manifold and a stable manifold; splicing the invariant manifolds of the Halo orbit using a point selection method, and calculating the optimal manifold splicing points of the comprehensive position deviation and velocity deviation; The optimal splicing orbit acquisition module is configured to construct a four-pulse dual-circle restricted four-body model; using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, using the particle swarm optimization algorithm and the sequential quadratic programming algorithm to solve the four-pulse dual-circle restricted four-body model to obtain the optimal splicing orbit; specifically comprising: using the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system as initial values, using the particle swarm optimization algorithm to solve the upper optimization model, and outputting the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system; wherein, in the particle swarm optimization algorithm, the optimal manifold splicing points of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the particle swarm optimization algorithm is used to solve the upper optimization model, and output the optimal perturbation factors of the Halo orbits corresponding to the Sun-Earth system and the Earth-Moon system are used. During the iterative optimization process of the method, a perturbation factor threshold is set. At the end of each iteration, it is judged whether the local optimal perturbation factor obtained each time is less than or equal to the perturbation factor threshold. If it is greater, the optimization result is penalized, that is, its performance index is increased. The optimal manifold splicing point of the Halo orbit corresponding to the Sun-Earth system and the Earth-Moon system, as well as the optimal perturbation factor of the Halo orbit corresponding to the Sun-Earth system and the Earth-Moon system are used as initial values, and the sequential quadratic programming algorithm is used to solve the lower-level optimization model, and the time it takes for the spacecraft in the Halo orbit of the Sun-Earth system and the Earth-Moon system to run from the final corrected initial value to the optimal manifold splicing point is output.

Images