Design and control method for relative motion configuration of earth-moon libration point spacecraft formation

Abstract

The application discloses a method for designing and controlling relative motion configuration of spacecraft formation in the Earth-Moon libration point orbit, and belongs to the technical field of spacecraft orbit dynamics and control. The method firstly constructs a double-level dynamics model composed of a circular restricted three-body problem model and a high-precision ephemeris model, and defines the orbit amplitude difference and phase difference as configuration parameters; then, according to different combinations of the orbit amplitude difference and phase difference, three types of natural relative motion configurations are set, periodic orbits are quickly generated under the circular restricted three-body problem model, and multiple turns are corrected by using the high-precision ephemeris model and two-stage multiple shooting; for the controlled tasks of fixed position, plane and space encirclement, a hierarchical constraint model is established, and a two-stage optimization strategy is adopted to solve the control pulse. The application realizes efficient design and stable control of the formation configuration in a complex dynamics environment, and is suitable for various libration point orbit tasks.

Description

Technical Field

[0001] This invention relates to the field of spacecraft orbital dynamics and control technology, specifically to the design and control method of relative motion configuration of spacecraft formations in the Earth-Moon translational point orbit. Background Technology

[0002] The Earth-Moon libration orbit is a special region of gravitational equilibrium within the Earth-Moon system, offering advantages such as low fuel consumption and long-term loiter capability. It has become a core orbital choice for missions such as deep space exploration, relay communication, and lunar polar observation. When spacecraft formations fly in libration orbits, multi-satellite coordination enables three-dimensional observation and data relay enhancement, significantly improving mission efficiency. Therefore, the design and control of the relative motion configuration of spacecraft formations in libration orbits is currently a research hotspot and engineering challenge in the field of aerospace dynamics.

[0003] Currently, research on translational point orbit spacecraft formations mainly focuses on natural and controlled configurations, but the following key issues remain: Preservation and solution of configuration under high-precision models are difficult: Existing relative motion and control designs mostly rely on periodic orbits and linearized structures in simplified models. However, under high-precision ephemeris models, multibody gravity and environmental perturbations introduce non-autonomous effects, which destroy the strict periodic orbit structure. This leads to enhanced coupling between relative geometric constraints and control timing, and the model is extremely sensitive to initial values ​​and disturbances, which significantly increases the difficulty and computational cost of solving multi-loop formation configurations.

[0004] The natural relative motion configuration law is not fully characterized: the natural configuration is induced by the amplitude difference and phase difference of the initial state of the spacecraft. Existing work mainly focuses on the analysis of a single orbital family or specific operating conditions. There is a lack of unified configuration maps and evolution laws for different translational point orbital families and different offset operating conditions. As a result, engineering selection and safety margin assessment still mainly rely on scattered calculations and empirical judgments, lacking systematic guidance.

[0005] Numerical stability issues in strongly nonlinear regions: In strongly nonlinear regions, especially in the near-lunar segment of the near-straight halo orbit NRHO, traditional numerical construction methods are prone to iterative non-convergence or numerical divergence. When the altitude of the near-lunar end of the orbit is low or the sensitivity of state transition increases, the convergence domain of splice point selection and differential correction shrinks significantly, reducing the reliability of batch generation of multi-cycle reference orbits and corresponding configurations, and limiting the coverage of the configuration database.

[0006] The contradiction between globality and optimality in solving controlled configurations: For controlled configurations such as fixed positions and forced planar encirclement, it is necessary to simultaneously determine the pulse time and pulse amplitude and satisfy multiple constraints such as relative distance, relative velocity and task cycle. Under high-precision models, the computational cost of directly performing global optimization is huge; while relying solely on surrogate models can improve efficiency, the accuracy in the neighborhood of the optimal solution is often insufficient, resulting in a complex solution process that is highly sensitive to algorithm parameters and initial guesses, making it difficult to balance global exploration and local optimization.

[0007] Therefore, there is an urgent need for a spacecraft formation relative motion configuration design and control method that can bridge the gap between simplified and high-precision models, systematically reveal the laws of natural configuration, possess high numerical stability, and take into account both global and local optimization. Summary of the Invention

[0008] To address the problems existing in the background art, this invention proposes a method for designing and controlling the relative motion configuration of spacecraft formations in a lunar translational point orbit, comprising the following steps: S1. For spacecraft formations in the Earth-Moon translation point orbit, the orbital amplitude difference and phase difference are defined as parameters describing the relative motion configuration of the formation. A two-level dynamic model is constructed, which includes a circular restricted three-body problem model and a high-precision ephemeris model. For different translation points and orbital families, a series of periodic orbits with different amplitudes are generated as the reference orbits for formation design. S2. Based on the aforementioned reference orbit, three types of natural relative motion configurations are defined according to different combinations of orbit amplitude differences and phase differences; S3. Rapidly generate multi-cycle periodic orbits and natural relative motion configurations under the circular restricted three-body problem model; S4. The results of the circular restricted three-body problem model are corrected for multiple orbits under a high-precision ephemeris model; S5. Establish a constrained optimization model of the controlled configuration based on the preset formation task requirements, and solve for the velocity increment required to maintain the configuration.

[0009] Specifically, in step S1, let Let represent the state vector of the spacecraft. Then, in the Earth-Moon rotating coordinate system, the normalized equations of motion for the circularly restricted three-body problem are expressed as: ; in, Represents the potential function Regarding respectively Partial derivatives, potential function Defined as: ; in, The Earth-Moon gravitational constant is... These represent the distances between the spacecraft and the two primary celestial bodies: , .

[0010] Specifically, the process of establishing the high-precision ephemeris model in step S1 is as follows: The motion equations of the spacecraft are established in the lunar-centered inertial coordinate system. These equations include a circularly restricted three-body gravitational term and additional perturbation terms. The additional perturbation terms specifically include solar gravitational perturbation, solar radiation pressure perturbation, and Earth's non-spherical gravitational perturbation. The perturbation generated by the solar gravitational perturbation is expressed as: ; in, The gravitational constant of the Sun. This represents the position vector from the Moon to the Sun. It is the position vector of the spacecraft in the lunar center J2000 inertial coordinate system; The calculation process for the solar radiation pressure perturbation term is as follows: (1) Calculate the solar radiation pressure based on the relative positions of the spacecraft, the sun, and the moon: ; in The speed of light in a vacuum. Total solar irradiance: ; in Solar brightness power, This represents the distance from the spacecraft to the sun. (2) Determine whether the spacecraft is in the lunar shadow region. If the conditions for a solar eclipse are met... If the spacecraft is in the shadow region, the light pressure perturbation is zero; otherwise, the light pressure perturbation is calculated. ; in, The distance from the center of the sun to the spacecraft. The distance from the center of the sun to the point of tangency on the moon. The angle between the line connecting the center of the sun to the spacecraft and the line connecting the center of the sun to the center of the moon. The angle between the line connecting the center of the Sun and the center of the Moon and the line connecting the center of the Sun and the point of tangency of the Moon. The cross-sectional area of ​​the spacecraft. This is the radiation pressure coefficient. The value is between 0 and 2, where 0 represents no reflection, 1 represents a black body, and 2 represents complete reflection. m For spacecraft mass; For the non-spherical gravitational perturbation term of the Earth, the perturbation acceleration of the spacecraft by the Earth's non-spherical gravity is calculated using the spherical harmonic function expansion of the Earth's gravitational field potential function. The spherical harmonic function expansion is expressed as: make It is extracted into term J2 to characterize the effect of Earth's non-spherical gravitational perturbation; The perturbation acceleration of the spacecraft caused by Earth's non-spherical gravity is then calculated using the following steps: a. Calculate the Earth's gravitational potential function with respect to radial distance. Earth's latitude Earth's longitude partial derivatives , , ; b. Express the partial derivatives in the spherical coordinate system using the coordinate components of the rectangular coordinate system: ; ; ; in, ; c. Multiply the partial derivatives by the gradient components and then synthesize them to obtain the gravitational acceleration of the spacecraft in the body-fixed coordinate system: ; d. Then calculate the differential acceleration: ; in, Let be the rotation matrix from the lunar inertial coordinate system to the Earth-fixed coordinate system. Let be the position vector of the spacecraft in the lunar center J2000 inertial coordinate system. This is the position vector pointing from the Moon's center of mass to the Earth's center of mass.

[0011] Specifically, the three types of natural relative motion configurations mentioned in step S2 are as follows: (1) Same amplitude but different phase condition: The master spacecraft and the slave spacecraft operate on the same translational periodic orbit, with a difference in amplitude. Furthermore, the two spacecraft have a non-zero initial phase difference. ; (2) Same phase, different amplitude condition: The master spacecraft and the slave spacecraft operate on two adjacent translational point periodic orbits in the family, and the two spacecraft have the same initial phase angle. Meanwhile, there is a non-zero orbital amplitude difference between the two spacecraft. ; (3) Different amplitude and different phase conditions: The master spacecraft and the slave spacecraft operate on two adjacent translational point periodic orbits in the family, and the two spacecraft have a non-zero initial phase difference. It also has a non-zero orbital amplitude difference. .

[0012] Specifically, step S3 includes: S31. Under the dynamic model of the circular restricted three-body problem, the initial states of the reference orbit and the target orbit are obtained by loading the reference orbit data and interpolating. The reference orbit is the flight orbit of the master spacecraft, and the target orbit is the flight orbit of the slave spacecraft. The corresponding initial multi-cycle periodic orbits are generated by numerical integration. S32. The initial multi-cycle periodic orbit is corrected and spliced ​​into a continuous orbit through a two-stage multiple-targeting method. S33. Calculate the relative state of the master and slave spacecraft in the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system. The relative state is obtained by subtracting the states of the master and slave spacecraft in the same time period. Output the maximum / minimum relative distance and the maximum relative coordinate component, and then generate a three-dimensional configuration diagram to characterize the natural relative motion configuration and geometric features. The maximum / minimum relative distance is the norm of the relative state.

[0013] Specifically, the two-stage multiple-targeting method described in step S32 includes: S321. Based on the preset orbital period and the number of calculated revolutions, discretize the entire orbit on the time axis and set a series of target points. Each firing point includes a position vector. Velocity vector and the corresponding flight time ; S322, First Stage Speed ​​Correction: Fix the position of all firing points and flight time For each firing point except the last firing point The initial velocity correction at this point is solved using a single-shot differential corrector. This makes the endpoint from that point and the next target point... Location The orbits overlap, thereby eliminating the spatial discontinuity in the trajectory at the firing point, and the initial velocity correction amount. The calculation formula is: ; in, This is the velocity-position sensitivity submatrix in the state transition matrix. This refers to the end position error; S323, Second Stage Spacetime Joint Correction: After completing the first stage velocity correction, by... Location of each firing point and time Find the partial derivative and construct a formula that includes the adjustment amount of all firing point positions. and time adjustment amount Global Jacobian matrix M The sum of the velocity discontinuities at each firing point of the trajectory With minimization as the objective, the target adjustment amount is calculated using the least squares method, i.e. = , in, for: ; Then, the position and time of each firing point are updated using the solved correction values, and the process returns to step S322 until the position and velocity discontinuities at all firing points converge to the preset tolerance range.

[0014] Specifically, step S4 is as follows: for the multi-cycle periodic orbit generated under the circular restricted three-body problem model in S3, the segmented state is further corrected using a two-stage multi-targeting method under a high-precision dynamic model containing perturbations, and then transformed to the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system. The absolute configuration diagrams of the two spacecraft are drawn. Finally, by interpolating the synchronous orbit time axis, the three-dimensional relative position between the orbit pair is calculated, the maximum relative coordinate component is extracted as the relative motion configuration performance index, and the relative configuration diagrams of the two spacecraft are drawn.

[0015] Specifically, in the process of using the two-stage multiple target method under the high-precision dynamic model, when the time of the periodic orbit is segmented according to the orbital period and the number of calculated orbits to obtain the initial value of the segment, for the near-straight halo orbit NRHO, sensitive region segment points located near the lunar perihelion are removed. The sensitive region segment points are specifically time nodes that have strong nonlinear characteristics in the dynamic environment and are prone to numerical iteration divergence in the differential correction process. Subsequently, the remaining segment point states are converted to the lunar inertial coordinate system as the initial guess values ​​for the two-stage multiple target firing.

[0016] Specifically, the formation task described in step S5 includes: (1) Fixed position formation: The spacecraft maintains a fixed relative position with respect to the host spacecraft in a selected coordinate system, and the relative position deviation is less than a preset distance threshold; (2) Planar forced orbit formation: The spacecraft moves around the host spacecraft along a closed polygonal path or a quasi-closed path in a specific plane of the selected coordinate system. The vertices of the polygon are the control points. The spacecraft orbits the host spacecraft by activating maneuvering pulses at the control points. (3) Three-dimensional spatial formation: the spacecraft moves around the host spacecraft in three-dimensional space and satisfies the ratio constraint between the maximum and minimum relative distances; Specifically, the solution process for the fixed-position formation controlled configuration is as follows: The reference orbit data is loaded and interpolated to obtain the initial state of the reference orbit. The initial state of the target orbit is constructed by adding a specific initial offset to the reference orbit. The corresponding initial multi-cycle periodic orbits are generated by numerical integration. Then, the velocity increment matching the configuration offset is obtained through a one-stage multiple target method, thereby obtaining the fixed position configuration diagram. The solution process for the planar forced orbiting formation controlled configuration and the space three-dimensional orbiting formation controlled configuration is as follows: The reference orbit data is loaded and interpolated to obtain the initial state of the reference orbit. An optimized offset is obtained through a two-stage optimization algorithm. The optimized offset is added to the reference orbit to obtain the standard orbit. Then, a one-stage multiple target method is used to obtain the velocity increment that matches the configuration, thereby obtaining the corresponding configuration diagram. The two-stage optimization algorithm, designed for forced planar encirclement formations and controlled configurations of three-dimensional spatial encirclement formations, specifically includes the following steps: (1) Determine design variables based on formation plane type: If the formation is restricted to a fixed coordinate plane, including the xy plane, xz plane, and yz plane, the design variables are: Corner position parameters , The angular position parameter controls the orientational distribution of the polygon vertices on the planar circumference, which is the number of configuration points. If the formation is in a spatial plane, the design variables are as follows: In addition to the angular position parameters, two additional planar orientation parameters are introduced. , used to determine the normal direction of the spatial formation plane; (2) Generate the target relative position of the spacecraft with respect to the host spacecraft based on design variables: For a fixed coordinate plane, the target's relative position is determined by the maximum relative distance of the formation. Determined by the angular position parameter, if it is in the xy plane, the target's relative position is specifically expressed as: ; If it is the xz / yz plane, then set the corresponding coordinate axis components to 0; For a spatial plane, the orientation parameter is used. Construct two orthogonal basis vectors in the plane Then, the relative position of the target is generated by combining the angular position parameters, that is: ; (3) Design the objective function with the goal of minimizing the cumulative velocity increment norm required for spacecraft to maintain formation: ; in To design the set of variables, For the first k The speed increment of the secondary maneuver; The constraints include: , A preset maneuverability threshold is set for the spacecraft; the relative position deviation is less than the maximum relative distance of the formation. ; (4) Solve through two-stage optimization: The first stage involves global optimization using a proxy. Within the range of design variable values, global sampling is performed. For each sample point, the following steps are executed: a. Generate the relative position of the target at each maneuver moment based on the design variables; b. Perform orbit propagation and position targeting under a high-precision ephemeris model to solve for the velocity increment required for each maneuver; c. Calculate the cumulative velocity increment and constraint violation, and iteratively update the proxy approximation model to guide the search to converge to a feasible region with a cost lower than the preset value, outputting candidate design variable solutions. The second stage involves local refinement of the sequential quadratic programming algorithm: using the candidate design variable solutions output in the first stage as initial values, the algorithm is used to optimize local constraints. A feasibility-first mechanism is enabled to suppress constraint violations during the iteration process, while further reducing the cumulative velocity increment. Finally, the diagonal position parameters and the plane orientation parameters in the spatial plane case are continuously corrected within the same design variable space until the convergence criterion is met, and the final design variable solution is output. Then, the target relative position at each maneuver moment is obtained through the final design variable solution as the optimization bias.

[0017] In summary, the beneficial technical effects of the present invention are as follows: 1. Achieved integrated closed-loop design of formation configuration under complex dynamic environments: This invention integrates the construction of a multi-cycle periodic orbit library, the classification and generation of natural relative motion configurations, the modeling of controlled configuration constraints, two-stage optimization and multi-target solution into a unified process, enabling comparable and reusable results output under the same computational framework for different orbital families, different configuration conditions and different task constraints, thereby significantly reducing the cost of scheme design and iteration and improving the efficiency of engineering implementation.

[0018] 2. By bridging the gap between simplified and high-precision models, the accuracy and engineering feasibility of configuration design are ensured: This invention achieves rapid configuration scanning and regularity extraction under a circular restricted three-body problem model, and introduces perturbation factors such as solar gravity, solar radiation pressure, and Earth's non-spherical gravity into a high-precision ephemeris model for multi-cycle propagation and correction. This forms a continuous migration path from theoretical analysis to the engineering environment, making configuration selection and control cost assessment closer to real mission conditions, and improving the credibility and applicability of the conclusions.

[0019] 3. A two-stage, multi-target differential correction strategy is proposed, which significantly improves the numerical stability in the strongly nonlinear region: In response to the problem that near-linear halo orbits such as NRHO have strong nonlinear characteristics in the near-lunar segment and are prone to iterative divergence, this invention effectively expands the convergence domain of differential correction by removing sensitive region segmentation points and adopting a two-stage correction strategy, reduces the sensitivity to initial values, and ensures the reliability of generating multi-cycle reference orbits and corresponding configurations, thereby supporting the construction of a larger configuration database.

[0020] 4. A hierarchical constraint and two-stage optimization solution mechanism for controlled configurations is established, taking into account both global exploration and local optimization: For various controlled tasks such as fixed position, forced planar encirclement and three-dimensional spatial encirclement, this invention embeds constraints into the optimization model and adopts a two-stage strategy that combines global initial screening with sequential secondary programming for local refinement. This solves the problem of high computational cost and easy getting trapped in local optima under high-precision models, and can efficiently obtain control pulse sequences that satisfy complex constraints and have a relatively low cost for velocity increment. Attached Figure Description

[0021] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is the solar eclipse condition judgment diagram in this invention; Figure 3 This is a schematic diagram of the first-stage target correction in this invention; Figure 4 This is a simulation result diagram of the two-stage multiple target firing method in the embodiment of the present invention; Figure 5 This is a simulation result diagram of extrapolating one circle using GMAT in an embodiment of the present invention; Figure 6 This is a comparison chart of the simulation results of the two-stage multiple-targeting method in this embodiment of the invention and the STK simulation results; Figure 7 This is a position error diagram comparing the simulation results of the two-stage multiple target firing method in this embodiment of the invention with the STK simulation results; Figure 8 This is a speed error diagram comparing the simulation results of the two-stage multiple-targeting method in this embodiment of the invention with the STK simulation results; Figure 9 This is a schematic diagram of fixed-position formation in this invention; Figure 10 This is a schematic diagram of the planar forced encirclement formation in this invention. Detailed Implementation

[0022] To make the technical means, creative features, objectives and effects of this invention clearer and easier to understand, the invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0023] Example like Figure 1 As shown in the embodiment of the present invention, the method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit specifically includes the following steps: S1. For spacecraft formations in the Earth-Moon libration point orbit, the formation adopts a master-slave framework with a master spacecraft and a slave spacecraft. The orbital amplitude difference and phase difference are defined as parameters describing the relative motion configuration of the formation. The relative distance, relative velocity, relative azimuth, and velocity increment are defined as parameters as controlled configuration constraints and output indicators, so that different orbital families, different operating conditions, and different control strategies are comparable under the same indicator system. For different translation points and orbit families, a series of periodic orbits with different amplitudes are generated as reference orbits for subsequent formation design. The specific data of the reference orbits generated in this part are existing data and are only used as a rough reference. Subsequent steps will be based on this reference for correction and optimization. Construct a two-level dynamic model comprising a circular restricted three-body problem model and a high-precision ephemeris model: (1) Circular restricted three-body problem model: Let Let represent the state vector of the spacecraft. Then, in the Earth-Moon rotating coordinate system, the normalized equations of motion for the circularly restricted three-body problem are expressed as: ; in, Represents the potential function Regarding respectively Partial derivatives, potential function Defined as: ; in, The Earth-Moon gravitational constant is... These represent the distances between the spacecraft and the two primary celestial bodies: , .

[0024] (2) High-precision ephemeris model: The motion equations of the spacecraft are established in the lunar inertial coordinate system, which include circular restricted three-body gravitational terms and additional perturbation terms. The additional perturbation terms specifically include solar gravitational perturbation, solar radiation pressure perturbation and Earth's non-spherical gravitational perturbation. The perturbation caused by solar gravitational perturbation is represented as follows: ; in, The gravitational constant of the Sun. This represents the position vector from the Moon to the Sun. It is the position vector of the spacecraft in the lunar center J2000 inertial coordinate system; Solar radiation pressure perturbation term: Solar radiation exerts stress on spacecraft; to model this stress, the following stress model is used: ; ; in, It is the total solar irradiance. It is the speed of light in a vacuum. It is solar brightness power. It is the distance to the sun. By considering a spherical satellite and only specular reflection, we can obtain: ; in, Its direction is opposite to that of the sun. It is the cross-sectional area of ​​the satellite, assumed to be circular. It is the radiation pressure coefficient. The value is between 0 and 2, where 0 represents no reflection, 1 represents a black body, and 2 represents complete reflection.

[0025] For solar radiation pressure perturbation, a solar eclipse condition needs to be met, meaning the spacecraft is not visible to the sun, implying it is behind the moon. In this case, the solar radiation pressure perturbation is zero. Figure 2 As shown, the conditions for lunar occultation of a spacecraft are: ; Right now ; Based on geometric considerations, we obtain: ; Solve : ; This can then be expressed as : ; ; Therefore, the conditions for a solar eclipse can be derived from geometric considerations.

[0026] Earth's non-spherical gravitational perturbation terms: First, the potential function of Earth's non-spherical gravitational field is given. The specific expression is: The Earth is decomposed into several small volume elements, each so small that it can be considered a point mass, with its mass denoted as . Its geocentric position vector is denoted as Therefore, there is ; This integration is performed over the entire Earth, of which... It is an extraterrestrial experimental particle Small volume elements inside the Earth The distance is ; ; here That is, the test particle With volume element The angle subtended by the Earth's center, by It has the following expansion: ; in Right now of Legendre polynomials have ; Substituting the above expression into the bitwise function yields... ; For the test particle and volume element geocentric position vector and In the Earth-fixed coordinate system, the relationship between rectangular coordinates and spherical coordinates is as follows: ; ; in and Given the geocentric latitude and longitude of the test particle and volume element in the Earth-fixed coordinate system, the spherical trigonometric relations can be used to give... ; According to the theory of spherical harmonics, Legendre polynomials It can also be displayed in the following forms: ; in yes The association of Legendre polynomials, when It then degenerates into a general Legendre polynomial. ,Will Substituting the expansion of the bitwise function, we get: If remember ; ; The above equation then becomes ; Two constants are introduced here. and The former refers to the total mass of the Earth, as indicated by the symbols used previously. The latter is the equatorial radius of the reference ellipsoid. Regarding the reference ellipsoid, the above coefficients... These gravitational potential coefficients cannot be directly given by an integral expression because there is currently no specific form to represent the shape and density distribution of the Earth. In fact, these coefficients are obtained through a combined adjustment of satellite geodesy and ground geodesy. Furthermore, along with these gravitational potential coefficients, a corresponding Earth reference ellipsoid, including the geocentric gravitational constant, is also provided. Equatorial radius and flatness The center of this reference ellipsoid (considered the "geocenter") is taken as the origin of the associated Earth-fixed coordinate system, and the internationally accepted orientation is CIO. Axial direction.

[0027] The above equation is the general form of the Earth's gravitational potential function expression. The first major term corresponds to a spherical celestial body, while the subsequent terms are non-spherical correction terms. The magnitude of the spherical correction reflects the irregularity of the shape and the non-uniformity of the internal mass distribution, all of which are determined by two sets of coefficients. and The numerical value reflects this.

[0028] As can be seen from the above Earth gravitational potential function, the non-spherical part of the Earth contains two terms with completely different properties: one corresponding to... ,at this time, , One type of term is unrelated to longitude; while the other type corresponds to... Dependent on longitude .

[0029] Distinguish between these two types of terms and write them in the following form: ; in It is the gravitational constant, and ; ; Therefore, the Earth's gravitational potential function can be written as: ; The terms mentioned above that are unrelated to longitude are called zone harmonic terms, while those that are related to longitude are called field harmonic terms; correspondingly (or This is called the band harmonic coefficient. and ( $ is called the field harmonic coefficient.

[0030] This invention makes , The J2 perturbation of Earth is captured, among which This represents the J2 term, the second-order spherical harmonic term. J2 is the main term of the Earth's non-spherical gravitational field, which is related to the Earth's equatorial bulge and typically affects the orbital plane and inclination. This term corresponds to the Earth's shape, i.e., the Earth's equatorial radius is greater than its polar radius. This indicates that the term is independent of longitude; The perturbation acceleration of the spacecraft caused by Earth's non-spherical gravity can then be solved using the following steps: a. Calculate the Earth's gravitational potential function with respect to radial distance. Earth's latitude Earth's longitude partial derivatives , , ; b. Express the partial derivatives in the spherical coordinate system using the coordinate components of the rectangular coordinate system: ; ; ; in, ; c. Multiply the partial derivatives by the gradient components and then synthesize them to obtain the gravitational acceleration of the spacecraft in the body-fixed coordinate system: ; d. Then calculate the differential acceleration: ; in, Let be the rotation matrix from the lunar inertial coordinate system to the Earth-fixed coordinate system. Let be the position vector of the spacecraft in the lunar center J2000 inertial coordinate system. This is the position vector pointing from the Moon's center of mass to the Earth's center of mass.

[0031] S2. Based on the reference orbit, three types of natural relative motion configurations are defined according to different combinations of orbital amplitude and phase differences. Due to the instability of translational periodic orbits, both the reference star (master spacecraft) and the orbiting star (slave spacecraft) must be operating in periodic orbits for their relative motion to form a natural configuration. The natural orbital configuration consists of two parts of relative motion: one is the phase difference between the reference star and the orbiting star. The relative motion is caused by two factors: firstly, the difference in z-axis amplitude between the reference star and the orbiting star; and secondly, the difference in z-axis amplitude between the reference star and the orbiting star. The resulting relative motion will , As a formation configuration parameter, the natural configuration characteristics of orbital formations are analyzed, and natural orbital formations are divided into three categories: (1) Same amplitude but different phase condition: The master spacecraft and the slave spacecraft operate on the same translational periodic orbit, with a difference in amplitude. Furthermore, the two spacecraft have a non-zero initial phase difference. ; (2) Same phase, different amplitude condition: The master spacecraft and the slave spacecraft operate on two adjacent translational point periodic orbits in the family, and the two spacecraft have the same initial phase angle. Meanwhile, there is a non-zero orbital amplitude difference between the two spacecraft. ; (3) Different amplitude and different phase conditions: The master spacecraft and the slave spacecraft operate on two adjacent translational point periodic orbits in the family, and the two spacecraft have a non-zero initial phase difference. It also has a non-zero orbital amplitude difference. .

[0032] For different orbits, its amplitude The definition is as follows: Halo orbit: the out-of-plane amplitude of the halo orbit. Defined as the maximum deviation of the track from the xy plane in three-dimensional space, that is, the maximum distance from any point on the track to the xy plane.

[0033] NRHO (Near-Linear Halo Orbit): Amplitude of NRHO Defined as the perpendicular distance from the perilunar point to the xy plane.

[0034] Vertical orbit: the vertical orbit's out-of-plane amplitude. Defined as the maximum vertical distance of the orbit relative to the xy plane; it reflects the intensity of the significant vertical figure-eight oscillation of this type of orbit.

[0035] Short-period orbit: Characteristic amplitude of a short-period orbit Defined as the distance from the translation point to the point where the orbit intersects the line connecting the translation point and the Earth.

[0036] The phase of an orbit is specifically defined as follows: when a point in time is selected, the phase is indicated by the proportion of the period that has elapsed since that selected point in time. Based on the different amplitudes and phases, two reference orbits and a target orbit can be obtained, and the relative motion orbit can be further obtained. The maximum relative coordinate component of the relative motion coordinate can be up to the kilometer level.

[0037] S3. Rapidly generate multi-cycle periodic orbits and natural relative motion configurations under a circularly restricted three-body problem model, specifically including: S31. Under the dynamic model of the circular restricted three-body problem, the initial states of the reference orbit and the target orbit are obtained by loading the reference orbit data and interpolating. The reference orbit is the flight orbit of the master spacecraft, and the target orbit is the flight orbit of the slave spacecraft. The corresponding initial multi-cycle periodic orbits are generated by numerical integration. S32. The initial multi-cycle orbit is corrected and spliced ​​into a continuous orbit through a two-stage multiple-targeting method. In the two-stage multiple-targeting method, the entire orbit is first discretized on the time axis according to the preset orbit period and the number of cycles, and a series of target points are set. Each firing point includes a position vector. Velocity vector and the corresponding flight time The second level (Level 2) assumes that all firing points are uniformly distributed in time. Level 1 corrects the velocity at each firing point, while Level 2 adjusts the position and time of the firing points. If the firing points are within the convergence region of the differential corrector, the process can quickly converge to a continuous trajectory. The following is a detailed description of the two-level iterative process: (1) First-stage speed correction (Level 1): Along each firing point on the track (except the last one), the velocity is adjusted using a single-shot differential corrector so that the end point of each trajectory segment aligns with its next firing point. After this layer is completed, the track is spatially continuous, but a velocity adjustment is still required at each firing point. Thrust correction is needed to maintain the trajectory.

[0038] like Figure 3 As shown, the spacecraft starts from its initial state. Depart, along the nominal track Running, in Includes initial position vector and velocity vector The goal is to make the orbit at a certain moment Reaching the specified final state This status includes the target location. and target speed solid line trajectory The initial trajectory is represented by the dashed line. The corrected trajectory arrives at the target position exactly at the target time (as shown by the bullseye in the figure).

[0039] Single target shooting utilizes state transition matrix To solve for the change in initial velocity To eliminate end position error The state transition matrix describes the evolution of the state perturbation over time, and its linearized expression is: or = ; because and Unconstrained, the above equation can be simplified, and the required initial velocity correction can be obtained as: ; (2) Second-stage spatiotemporal joint correction (Level 2): In this layer, the position and timing of all firing points (including the last point) are adjusted, and the least squares method is used to reduce the total... After this step, the trajectory becomes discontinuous at each firing point, but its overall thrust overhead is smaller, which is beneficial for the convergence of the next Level 1 round. The method assumes that the trajectory consists of at least three firing points, and the goal is to reduce the total thrust overhead. Expenses.

[0040] Specifically, after completing the first stage of speed correction, by... Location of each firing point and time Find the partial derivative and construct a formula that includes the adjustment amount of all firing point positions. and time adjustment amount Global Jacobian matrix M The sum of the velocity discontinuities at each firing point of the trajectory With minimization as the objective, the target adjustment amount is calculated using the least squares method, i.e. = , in, for: ; Then, the position and time of each firing point are updated using the solved correction values, and the execution steps Level 1 are returned until the position and velocity discontinuities at all firing points converge to the preset tolerance range.

[0041] S33. Calculate the relative state of the master and slave spacecraft in the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system. The relative state is obtained by subtracting the states of the master and slave spacecraft in the same time period. Output the maximum / minimum relative distance and the maximum relative coordinate component, and then generate a three-dimensional configuration diagram to characterize the natural relative motion configuration and geometric features. The maximum / minimum relative distance is the norm of the relative state.

[0042] S4. The results of the circular restricted three-body problem model are corrected for multiple orbits under a high-precision ephemeris model. Specifically, for the multiple periodic orbits generated under the circular restricted three-body problem model in S3, the segmented states are further corrected using a two-stage multi-targeting method under a high-precision dynamic model that includes perturbations. The orbits are then converted to the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system. The absolute configuration diagrams of the two spacecraft are drawn. Finally, by interpolating the synchronous orbit time axis, the three-dimensional relative position between the orbit pairs is calculated. The maximum relative coordinate component is extracted as the relative motion configuration performance index, and the relative configuration diagrams of the two spacecraft are drawn.

[0043] Preferably, in the process of using the two-stage multiple target method under the high-precision dynamic model, when the time of the periodic orbit is segmented according to the orbital period and the number of calculated orbits to obtain the initial value of the segment, for the near-straight halo orbit NRHO, the sensitive region segment points located near the lunar perihelion are removed. The sensitive region segment points are specifically the time nodes that have strong nonlinear characteristics in the dynamic environment and are prone to numerical iteration divergence in the differential correction process. Subsequently, the state of the remaining segment points is transformed to the lunar center inertial coordinate system as the initial guess value of the two-stage multiple target.

[0044] Taking the southward NRHO orbit of the L2 translation point in a high-precision model with an amplitude of 3000km as an example, where the comprehensive parameter of light pressure multiplied by the surface-to-mass ratio is taken as 0.005, the method is compared with STK to illustrate its correctness: like Figure 4 The figure shows the simulation results obtained by using the two-stage multiple firing method of this invention. Figure 5 The image shows the simulation results of extrapolating one lap using GMAT, as shown below. Figure 6 The figure shown is a comparison between the simulation results of this method and the simulation results of STK. Figure 7 , Figure 8The figure shows the position error and velocity error compared with the simulation results of this method and the STK simulation results. It can be seen from the figure that the trajectories of the two are basically the same. When STK considers more perturbations than this method, the final error figure shows that the position error is within 10km and the velocity error is within 4.5e-4.

[0045] S5. Establish a constrained optimization model for the controlled configuration based on the preset formation task requirements, and solve for the velocity increment required to maintain the configuration: (1) Fixed-position formation: such as Figure 9 The diagram shows a fixed-position formation, where a spacecraft maintains a fixed relative position with respect to the host spacecraft in a selected coordinate system, and the relative position deviation is less than a preset distance threshold. Specifically, the spacecraft maintains a fixed relative position with respect to the host spacecraft in a designated coordinate system. The deviation range is divided into two cases: less than 10% and 20% of the relative distance. The relatively fixed position must include the positive and negative directions of the three coordinate axes XYZ. The reference coordinate system includes the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system. The relative distance needs to consider cases such as 10 km, 100 km, 200 km, 500 km, and 1000 km.

[0046] The solution process for the fixed-position formation controlled configuration is as follows: The reference orbit data is loaded and interpolated to obtain the initial state of the reference orbit. The initial state of the target orbit is constructed by adding a specific initial offset to the reference orbit. The corresponding initial multi-cycle periodic orbits are generated by numerical integration. Then, the velocity increment matching the configuration offset is obtained through a one-stage multiple target method, thereby obtaining the fixed position configuration diagram. (2) Planar forced encirclement formation: such as Figure 10 The diagram shows a forced planar orbit formation. The slave spacecraft moves around the master spacecraft along a closed polygonal path or a quasi-closed path within a specific plane of a selected coordinate system. The vertices of the polygon are control points. The slave spacecraft achieves orbiting of the master spacecraft by generating maneuvering pulses at the control points. Specifically, the reference coordinate system includes the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system. The configuration needs to include triangles, quadrilaterals, hexagons, and octagons, and is not necessarily a regular polygon. The control points are at the vertices of the polygons. The trajectory does not need to be a straight line. The plane needs to include three coordinate planes and one spatial plane. The control quantity of the spatial plane configuration needs to be smaller than that of the coordinate plane configuration. The orbital period includes three cases: 1.0, 0.5, and 0.25 orbital periods. The maximum distance of the polygon vertices from the reference star (main spacecraft) needs to consider cases such as 10 km, 100 km, 200 km, 500 km, and 1000 km. The minimum relative distance needs to be greater than 10% of the maximum relative distance.

[0047] (3) Three-dimensional spatial formation: the spacecraft moves around the host spacecraft in three-dimensional space and satisfies the ratio constraint between the maximum and minimum relative distances; Specifically, the scenario requires the orbiting satellite (slave spacecraft) to perform orbital control of the reference satellite (master spacecraft) within its orbital period in space. The design of the velocity increment and optimal orbital path around the reference satellite is to use pulse maneuvering and continuous thrust control. The velocity increment must be superior to planar orbital configurations of the same type, amplitude, and orbital period. The orbital period must include three cases: 1.0, 0.5, and 0.25 orbital periods. The maximum distance from the reference satellite must consider cases such as 10 km, 100 km, 200 km, 500 km, and 1000 km. The minimum relative distance must be greater than 10% of the maximum relative distance.

[0048] The solution process for the planar forced orbiting formation controlled configuration and the space three-dimensional orbiting formation controlled configuration is as follows: The reference orbit data is loaded and interpolated to obtain the initial state of the reference orbit. An optimized offset is obtained through a two-stage optimization algorithm. The optimized offset is added to the reference orbit to obtain the standard orbit. Then, a one-stage multiple target method is used to obtain the velocity increment that matches the configuration, thereby obtaining the corresponding configuration diagram. The two-stage optimization algorithm is designed for forced planar encirclement formations and controlled configurations of three-dimensional encirclement formations in space. Specifically, it includes the following steps: (1) Determine design variables based on formation plane type: If the formation is restricted to a fixed coordinate plane, including the xy plane, xz plane, and yz plane, the design variables are: Corner position parameters , The angular position parameter controls the orientational distribution of the polygon vertices on the planar circumference, which is the number of configuration points. If the formation is in a spatial plane, the design variables are as follows: In addition to the angular position parameters, two additional planar orientation parameters are introduced. , used to determine the normal direction of the spatial formation plane; (2) Generate the target relative position of the spacecraft with respect to the host spacecraft based on design variables: For a fixed coordinate plane, the target's relative position is determined by the maximum relative distance of the formation. Determined by the angular position parameter, if it is in the xy plane, the target's relative position is specifically expressed as: ; If it is the xz / yz plane, then set the corresponding coordinate axis components to 0; For a spatial plane, the orientation parameter is used. Construct two orthogonal basis vectors in the plane Then, the relative position of the target is generated by combining the angular position parameters, that is: ; (3) Design the objective function with the goal of minimizing the cumulative velocity increment norm required for spacecraft to maintain formation: ; in To design the set of variables, For the first k The speed increment of the secondary maneuver; The constraints include: , A preset maneuverability threshold is set for the spacecraft; the relative position deviation is less than the maximum relative distance of the formation. ; (4) Solve through two-stage optimization: The first stage involves global optimization using a proxy. Within the range of design variable values, global sampling is performed. For each sample point, the following steps are executed: a. Generate the relative position of the target at each maneuver moment based on the design variables; b. Perform orbit propagation and position targeting under a high-precision ephemeris model to solve for the velocity increment required for each maneuver; c. Calculate the cumulative velocity increment and constraint violation, and iteratively update the proxy approximation model to guide the search to converge to a feasible region with a cost lower than the preset value, outputting candidate design variable solutions. The second stage involves local refinement of the sequential quadratic programming algorithm: using the candidate design variable solutions output in the first stage as initial values, the algorithm is used to optimize local constraints. A feasibility-first mechanism is enabled to suppress constraint violations during the iteration process, while further reducing the cumulative velocity increment. Finally, the diagonal position parameters and the plane orientation parameters in the spatial plane case are continuously corrected within the same design variable space until the convergence criterion is met, and the final design variable solution is output. Then, the target relative position at each maneuver moment is obtained through the final design variable solution as the optimization bias.

[0049] Therefore, the method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit provided by this invention rapidly generates multi-cycle orbits and natural configurations by constructing a two-level dynamic framework of a circular restricted three-body problem and a high-precision ephemeris model. It overcomes the numerical convergence bottleneck in the strongly nonlinear orbit region by employing a two-stage multi-target differential correction strategy and establishes a hierarchical constraint model and a two-stage optimization mechanism to solve the control impulses of the controlled configuration. This achieves an integrated closed-loop design process of "basic environment construction - natural configuration classification and generation - high-precision multi-cycle correction - controlled configuration optimization solution". It solves the pain points of traditional translational orbit formation design, which suffers from low design efficiency, poor engineering adaptability, insufficient configuration reliability, and high control solution cost due to the inability to connect simplified and high-precision models, lack of systematic representation of natural configuration laws, easy divergence in iteration in strongly nonlinear regions, and difficulty in balancing global exploration and local optimization in controlled configurations. It provides an efficient, accurate, and systematic solution for the design and control of spacecraft formation configurations in the complex multibody gravitational environment of translational orbit space missions.

[0050] The embodiments of the present invention have been described in detail above with reference to the accompanying drawings. The above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for designing and controlling the relative motion configuration of spacecraft formations in a lunar translational orbit, characterized in that... Includes the following steps: S1. For spacecraft formations in the Earth-Moon translation point orbit, the orbital amplitude difference and phase difference are defined as parameters describing the relative motion configuration of the formation. A two-level dynamic model is constructed, which includes a circular restricted three-body problem model and a high-precision ephemeris model. For different translation points and orbital families, a series of periodic orbits with different amplitudes are generated as the reference orbits for formation design. S2. Based on the aforementioned reference orbit, three types of natural relative motion configurations are defined according to different combinations of orbit amplitude differences and phase differences; S3. Rapidly generate multi-cycle periodic orbits and natural relative motion configurations under the circular restricted three-body problem model; S4. The results of the circular restricted three-body problem model are corrected for multiple orbits under a high-precision ephemeris model; S5. Establish a constrained optimization model of the controlled configuration based on the preset formation task requirements, and solve for the velocity increment required to maintain the configuration. Step S3 is as follows: S31. Under the dynamic model of the circular restricted three-body problem, the initial states of the reference orbit and the target orbit are obtained by loading the reference orbit data and interpolating. The reference orbit is the flight orbit of the master spacecraft, and the target orbit is the flight orbit of the slave spacecraft. The corresponding initial multi-cycle periodic orbits are generated by numerical integration. S32. The initial multi-cycle periodic orbit is corrected and spliced ​​into a continuous orbit through a two-stage multiple-targeting method. S33. Calculate the relative state of the master and slave spacecraft in the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system. The relative state is obtained by subtracting the states of the master spacecraft and the slave spacecraft in the same time period. Output the maximum / minimum relative distance and the maximum relative coordinate component, and then generate a three-dimensional configuration diagram to characterize the natural relative motion configuration and geometric features. The maximum / minimum relative distance is the norm of the relative state. The two-stage multiple-target firing method described in step S32 specifically includes: S321. Based on the preset orbital period and the number of calculated revolutions, discretize the entire orbit on the time axis and set a series of target points. Each firing point includes a position vector. Velocity vector and the corresponding flight time ; S322, First Stage Speed ​​Correction: Fix the position of all firing points and flight time For each firing point except the last firing point The initial velocity correction at this point is solved using a single-shot differential corrector. This makes the endpoint from that point and the next target point... Location The orbits overlap, thereby eliminating the spatial discontinuity in the trajectory at the firing point, and the initial velocity correction amount. The calculation formula is: ； in, This is the velocity-position sensitivity submatrix in the state transition matrix. This refers to the end position error; S323, Second Stage Spacetime Joint Correction: After completing the first stage velocity correction, by... Location of each firing point and time Find the partial derivative and construct a formula that includes the adjustment amount of all firing point positions. and time adjustment amount Global Jacobian matrix M The sum of the velocity discontinuities at each firing point of the trajectory With minimization as the objective, the target adjustment amount is calculated using the least squares method, i.e. = ， in, for: ； Then, the position and time of each firing point are updated using the solved correction values, and the process returns to step S322 until the position and velocity discontinuities at all firing points converge to the preset tolerance range.

2. The method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit as described in claim 1, characterized in that, In step S1, let Let represent the state vector of the spacecraft. Then, in the Earth-Moon rotating coordinate system, the normalized equations of motion for the circularly restricted three-body problem are expressed as: ； in, Represents the potential function Regarding respectively Partial derivatives, potential function Defined as: ； in, The Earth-Moon gravitational constant is... These represent the distances between the spacecraft and the two primary celestial bodies: , .

3. The method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit as described in claim 1, characterized in that, The process of establishing the high-precision ephemeris model in step S1 is as follows: The motion equations of the spacecraft are established in the lunar-centered inertial coordinate system. These equations include a circularly restricted three-body gravitational term and additional perturbation terms. Specifically, the additional perturbation terms include solar gravitational perturbation, solar radiation pressure perturbation, and Earth's non-spherical gravitational perturbation. The perturbation generated by the solar gravitational perturbation is expressed as: ； in, The gravitational constant of the Sun. This represents the position vector from the Moon to the Sun. It is the position vector of the spacecraft in the lunar center J2000 inertial coordinate system; The calculation process for the solar radiation pressure perturbation term is as follows: (1) Calculate the solar radiation pressure based on the relative positions of the spacecraft, the sun, and the moon: ； in The speed of light in a vacuum. Total solar irradiance: ； in Solar brightness power, This represents the distance from the spacecraft to the sun. (2) Determine whether the spacecraft is in the lunar shadow region. If the conditions for a solar eclipse are met... If the spacecraft is in the shadow region, the light pressure perturbation is zero; otherwise, the light pressure perturbation is calculated. ； in, The distance from the center of the sun to the spacecraft. The distance from the center of the sun to the point of tangency on the moon. The angle between the line connecting the center of the sun to the spacecraft and the line connecting the center of the sun to the center of the moon. The angle between the line connecting the center of the Sun and the center of the Moon and the line connecting the center of the Sun and the point of tangency of the Moon. The cross-sectional area of ​​the spacecraft. This is the radiation pressure coefficient. The value is between 0 and 2, where 0 represents no reflection, 1 represents a black body, and 2 represents complete reflection. m For spacecraft mass; For the non-spherical gravitational perturbation term of the Earth, the perturbation acceleration of the spacecraft by the Earth's non-spherical gravity is calculated using the spherical harmonic function expansion of the Earth's gravitational field potential function. The spherical harmonic function expansion is expressed as: in, It is the gravitational constant. The radial distance from the spacecraft's center of mass to the Earth's center. The equatorial radius of the Earth's reference ellipsoid. Let be the order of the spherical harmonic expansion. for Band harmonic coefficients, For Legendre polynomials of the first kind, The latitude of the spacecraft's center in the Earth-fixed coordinate system. Let be the number of spherical harmonic expansions. for Step Legendre polynomials of the first kind This represents the geocentric longitude of the spacecraft in the Earth-fixed coordinate system. for Step The second cosine term harmonic coefficient for Step The second sine term harmonic coefficient; make It is extracted into term J2 to characterize the effect of Earth's non-spherical gravitational perturbation; The perturbation acceleration of the spacecraft caused by Earth's non-spherical gravity is then calculated using the following steps: a. Calculate the Earth's gravitational potential function with respect to radial distance. Earth's latitude Earth's longitude partial derivatives , , ; b. Express the partial derivatives in the spherical coordinate system using the coordinate components of the rectangular coordinate system: ； ； ； in, ; c. Multiply the partial derivatives by the gradient components and then synthesize them to obtain the gravitational acceleration of the spacecraft in the body-fixed coordinate system: ； d. Then calculate the differential acceleration: ； in, Let be the rotation matrix from the lunar inertial coordinate system to the Earth-fixed coordinate system. Let be the position vector of the spacecraft in the lunar center J2000 inertial coordinate system. This is the position vector pointing from the Moon's center of mass to the Earth's center of mass.

4. The method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit as described in claim 1, characterized in that, The three types of natural relative motion configurations mentioned in step S2 are as follows: (1) Same amplitude but different phase condition: The master spacecraft and the slave spacecraft operate on the same translational periodic orbit, with a difference in amplitude. Furthermore, the two spacecraft have a non-zero initial phase difference. ; (2) Same phase, different amplitude condition: The master spacecraft and the slave spacecraft operate on two adjacent translational point periodic orbits in the family, and the two spacecraft have the same initial phase angle. Meanwhile, there is a non-zero orbital amplitude difference between the two spacecraft. ; (3) Different amplitude and different phase conditions: The master spacecraft and the slave spacecraft operate on two adjacent translational point periodic orbits in the family, and the two spacecraft have a non-zero initial phase difference. It also has a non-zero orbital amplitude difference. .

5. The method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit according to claim 1, characterized in that, Step S4 specifically involves: for the multi-cycle periodic orbits generated under the circular restricted three-body problem model in S3, further correcting the segmented states using a two-stage multi-targeting method under a high-precision dynamic model that includes perturbations, and transforming them to the Earth-Moon rotating coordinate system and the Sun-Earth rotating coordinate system, drawing the absolute configuration diagrams of the two spacecraft, and finally, calculating the three-dimensional relative position between the orbit pairs by interpolating the synchronous orbit time axis, extracting the maximum relative coordinate component as the relative motion configuration performance index, and drawing the relative configuration diagrams of the two spacecraft.

6. The method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit as described in claim 5, characterized in that: In the process of using a two-stage multiple target method under a high-precision dynamic model, when the time of the periodic orbit is segmented according to the orbital period and the number of calculated orbits to obtain the initial value of the segment, for the near-straight halo orbit NRHO, sensitive region segment points located near the lunar perihelion are removed. The sensitive region segment points are specifically time nodes that have strong nonlinear characteristics in the dynamic environment and cause the numerical iteration of the differential correction process to diverge. Subsequently, the remaining segment point states are converted to the lunar inertial coordinate system as the initial guess values ​​for the two-stage multiple target firing.

7. The method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit as described in claim 6, characterized in that, The formation task mentioned in step S5 includes: (1) Fixed position formation: The spacecraft maintains a fixed relative position with respect to the host spacecraft in a selected coordinate system, and the relative position deviation is less than a preset distance threshold; (2) Planar forced orbit formation: The spacecraft moves around the host spacecraft along a closed polygonal path or a quasi-closed path in a specific plane of the selected coordinate system. The vertices of the polygon are the control points. The spacecraft orbits the host spacecraft by activating maneuvering pulses at the control points. (3) Three-dimensional orbital formation: the spacecraft moves around the host spacecraft in three-dimensional space and satisfies the ratio constraint between the maximum and minimum relative distances.

8. The method for designing and controlling the relative motion configuration of spacecraft formations in the Earth-Moon translational orbit as described in claim 7, characterized in that, The solution process for the fixed-position formation controlled configuration is as follows: The reference orbit data is loaded and interpolated to obtain the initial state of the reference orbit. The initial state of the target orbit is constructed by adding a specific initial offset to the reference orbit. The corresponding initial multi-cycle periodic orbits are generated by numerical integration. Then, the velocity increment matching the configuration offset is obtained through a one-stage multiple target method, thereby obtaining the fixed position configuration diagram. The solution process for the planar forced orbiting formation controlled configuration and the space three-dimensional orbiting formation controlled configuration is as follows: The reference orbit data is loaded and interpolated to obtain the initial state of the reference orbit. An optimized offset is obtained through a two-stage optimization algorithm. The optimized offset is added to the reference orbit to obtain the standard orbit. Then, a one-stage multiple target method is used to obtain the velocity increment that matches the configuration, thereby obtaining the corresponding configuration diagram. The two-stage optimization algorithm for forced planar orbit formation and controlled configurations of three-dimensional spatial orbit formation specifically includes the following steps: (1) Determine design variables based on formation plane type: If the formation is restricted to a fixed coordinate plane, including the xy plane, xz plane, and yz plane, the design variables are: Corner position parameters , The angular position parameter controls the orientational distribution of the polygon vertices on the planar circumference, which is the number of configuration points. If the formation is in a spatial plane, the design variables are as follows: In addition to the angular position parameters, two additional planar orientation parameters are introduced. , used to determine the normal direction of the spatial formation plane; (2) Generate the target relative position of the spacecraft with respect to the host spacecraft based on design variables: For a fixed coordinate plane, the target's relative position is determined by the maximum relative distance of the formation. Determined by the angular position parameter, if it is in the xy plane, the target's relative position is specifically expressed as: ； If it is the xz / yz plane, then set the corresponding coordinate axis components to 0; For a spatial plane, the orientation parameter is used. Construct two orthogonal basis vectors in the plane Then, the relative position of the target is generated by combining the angular position parameters, that is: ； (3) Design the objective function with the goal of minimizing the cumulative velocity increment norm required for spacecraft to maintain formation: ； in To design the set of variables, For the first k The speed increment of the secondary maneuver; The constraints include: , A preset maneuverability threshold is set for the spacecraft; the relative position deviation is less than the maximum relative distance of the formation. ; (4) Solve through two-stage optimization: The first stage involves global optimization using a proxy. Within the range of design variable values, global sampling is performed. For each sample point, the following steps are executed: a. Generate the relative position of the target at each maneuver moment based on the design variables; b. Perform orbit propagation and position targeting under a high-precision ephemeris model to solve for the velocity increment required for each maneuver; c. Calculate the cumulative velocity increment and constraint violation, and iteratively update the proxy approximation model to guide the search to converge to a feasible region with a cost lower than the preset value, outputting candidate design variable solutions. The second stage involves local refinement of the sequential quadratic programming algorithm: using the candidate design variable solutions output in the first stage as initial values, the algorithm is used to optimize local constraints. A feasibility-first mechanism is enabled to suppress constraint violations during the iteration process, while further reducing the cumulative velocity increment. Finally, the diagonal position parameters and the plane orientation parameters in the spatial plane case are continuously corrected within the same design variable space until the convergence criterion is met, and the final design variable solution is output. Then, the target relative position at each maneuver moment is obtained through the final design variable solution as the optimization bias.

Images