Double-layer optimization method and system for direct transfer between earth-moon special orbits under ephemeris model

Abstract

The invention discloses a double-layer optimization method and system for direct transfer between earth-moon special orbits under an ephemeris model, relates to the technical field of detector transfer design between the earth-moon special orbits, and aims to solve the technical problem that a traditional method cannot be suitable for direct transfer between the earth-moon special orbits with higher precision requirements. According to the technical key points, the method comprises the following steps: constructing a high-precision ephemeris dynamic model, and correcting an obtained discrete orbit data initial value by using the high-precision ephemeris dynamic model; constructing a direct transfer optimization model based on the corrected discrete orbit data with the time label and a high-precision ephemeris dynamical model; splitting and reconstructing the direct transfer optimization model; and iteratively solving the split and reconstructed direct transfer double-layer optimization model by using a multivariate universe optimization algorithm and a sequential quadratic programming algorithm with adjustment capability to obtain an index evolution iteration curve and a transfer orbit curve. According to the method, the optimal design of the direct transfer orbit between the earth and moon special orbits under the ephemeris model with higher precision requirements is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of design of probe transfer between special earth-moon orbits, and in particular to a double-layer optimization method and system for direct transfer between special earth-moon orbits under an ephemeris model. Background Art

[0002] With the deepening of human space exploration activities, especially the growing demand for scientific research, resource development and future construction of lunar space stations in the Earth-Moon space, it is extremely important to achieve efficient, energy-saving and precise transfer in the Earth-Moon space. Among them, since the special orbits under the circular restricted three-body problem dynamics model of the Earth-Moon system such as HALO, NRHO and DRO are used as long-term orbits for lunar probes and lunar space stations in the lunar exploration plans of various countries, the direct transfer problem between these special Earth-Moon orbits has gradually become a research hotspot. However, the design of transfer orbits between special Earth-Moon orbits under the ephemeris model faces challenges such as complex multi-body gravitational environment, strict energy constraints and flight window restrictions. Traditional analytical methods and intelligent methods are often difficult to meet actual engineering needs, and there is also a technical problem that high-precision direct transfer cannot be achieved on special Earth-Moon orbits.

[0003] In order to solve the transfer problem between special orbits in the Earth-Moon system, the literature [1] proposed a two-stage "search-optimization" method suitable for the design of transfer orbits between NRHO and DRO. In the search stage, the traversal method is used to find the starting point that can approach the target orbit after applying the deorbit impulse. In the optimization stage, the pulse solution is performed by direct shooting method, and the initial and final orbit constraints are constructed by fitting. Although this method can effectively solve the transfer orbit between NRHO and DRO, it is not universal. The traversal method in the search stage greatly affects the algorithm calculation efficiency. In addition, the time window problem is not considered, which cannot meet the requirements of high-precision ephemeris model and actual engineering needs. Literature [2] conducted a dynamic analysis of NRHO and DRO and designed a transfer orbit based on the higher period orbit, nearby manifold and resonant orbit of NRHO orbit, and carried out conversion and verification under the high-precision ephemeris model. Although it realizes the transfer orbit design under the high-precision ephemeris model, its design method has high specificity and low universality due to the dynamic characteristics near the NRHO orbit. In addition, this method is mainly used for continuous thrust transfer orbit design and is not suitable for direct transfer problems between various special orbits. In order to solve the orbit transfer problem under the ephemeris model, reference [3] designed a Earth-Moon transfer orbit based on a free return orbit taking into account the ephemeris model for the lunar polar exploration problem. The approach from the two-body model to the ephemeris model was achieved through the improved multi-conic curve splicing method. Although this method can meet the dynamic accuracy requirements, due to the dynamic characteristics of special orbits, the idea of ​​correcting from the two-body to the ephemeris model is not applicable to the transfer problem between special orbits. At the same time, the transfer time window problem cannot be solved. Reference [4] also designed an Earth-Moon transfer orbit based on a free return orbit under the ephemeris model for the needs of manned lunar landing. The difference from reference [3] is that its free return orbit is corrected to meet the high-precision ephemeris model after solving the orbit from the circular restricted three-body problem dynamics model. The pulse transfer strategy is still corrected from the two-body model to meet the ephemeris model. A two-round successive correction strategy is proposed. The interior point method and SQP sequence quadratic programming algorithm are used to solve the pulse transfer strategy and transfer orbit. This method still fails to resolve the limitation of the two-body model as the initial value in the design and solution of the direct transfer orbit, and cannot be applied to the transfer problem between special Earth-Moon orbits with higher precision requirements.

[0004] In summary, although the above-mentioned existing optimization method for transfer between special orbits in the Earth-Moon system has a certain effect in realizing orbit transfer, the method used is still a traditional method. The traditional algorithm has low accuracy and cannot solve the technical problem of direct transfer between special orbits of the Earth and the Moon. Summary of the invention

[0005] In view of this, the present invention proposes a double-layer optimization method and system for direct transfer between special Earth-Moon orbits under an ephemeris model, so as to solve the technical problem that the traditional design method of transfer orbits between Earth-Moon orbits under an ephemeris model cannot be applied to direct transfer between special Earth-Moon orbits with higher precision requirements.

[0006] According to one aspect of the present invention, a two-layer optimization method for direct transfer between special Earth-Moon orbits under an ephemeris model is proposed, comprising the following steps:

[0007] Step 1: construct a high-precision ephemeris dynamics model, and use the high-precision ephemeris dynamics model to correct the acquired initial values ​​of discrete orbit data to obtain corrected discrete orbit data with time tags;

[0008] Step 2: construct a direct transfer optimization model based on the corrected discrete orbit data with time tags and the high-precision ephemeris dynamics model;

[0009] Step 3: split and reconstruct the direct transfer optimization model to obtain a direct transfer two-layer optimization model after split and reconstructing;

[0010] Step 4: Use a multiverse optimization algorithm with adjustment capability and a sequential quadratic programming algorithm to iteratively solve the split and reconstructed direct transfer two-layer optimization model until the expected index or the maximum number of iterations is reached, and obtain the index evolution iteration curve and the transfer trajectory curve.

[0011] Furthermore, the discrete track data in step 1 is track data between the starting track and the target track, including position data, speed data and time tags.

[0012] Furthermore, the high-precision ephemeris dynamics model described in step 1 is established as follows:

[0013]

[0014] In the formula, R EI =[X EI ,Y EI ,Z EI ] T represents the position vector of the spacecraft in the J2000 inertial coordinate system; μ E represents the gravitational constant of the Earth; + indicates the gravitational acceleration at the center of the earth; A N A is the perturbation acceleration caused by the gravity of other planets in the solar system outside the main celestial body; NSE A is the non-spherical perturbation acceleration of the main celestial body; R A is the solar pressure perturbation acceleration; D is the atmospheric drag perturbation acceleration.

[0015] Furthermore, in step 1, the high-precision ephemeris dynamics model is used to correct the initial value of the acquired discrete orbit data, including continuous position correction and continuous velocity correction; wherein,

[0016] The continuous position correction is: using the position differential correction equation to correct the initial value of the discrete orbit data obtained to obtain discrete orbit data with continuous position; the position differential correction equation is:

[0017] δv i0 =M r δr if

[0018] In the formula, represents the integral end point obtained after orbit recursion, i=1,2,...,n-1; represents the starting point of orbital integration; δr if Indicates the deviation between the i-th point integral end position and the i+1-th point integral start position; δv i0 Indicates the change in velocity at the starting point of the i-th integration point;

[0019] The speed continuous correction is: the trajectory after the position continuous correction has a speed jump at each segment connection point, and the jump is corrected by using the speed differential correction equation to obtain discrete track data with continuous speed; the speed differential correction equation is:

[0020]

[0021] In the formula, δr j0 represents the change in the starting point of the j-th point, j = 2, ..., n-1; δv kf represents the deviation between the velocity at the end of the integration at the kth point and the velocity at the start of the integration at the k+1th point, k=1,2,...,n-1; r j0 Indicates the starting point of the j-th point integration;

[0022] After the speed continuous correction is completed, the position of each end point is changed, and the position correction is continued; the iterative correction is repeated until the position error and the speed error reach the preset threshold value, and the discrete track data with continuous position and speed after correction is obtained.

[0023] Furthermore, step one also includes: after obtaining the discrete orbit data with continuous corrected position and velocity, adding a time tag to the discrete orbit data with continuous corrected position and velocity, wherein the time tag is adding a real ephemeris epoch as a time tag according to the integration step, thereby obtaining the discrete orbit data with continuous position and velocity with the time tag.

[0024] Furthermore, the direct transfer optimization model in step 2 is:

[0025] The objective function is to minimize the fuel J1: Energy J2 minimization:

[0026] Constraints:

[0027] Where N is the total number of pulses; ||Δv i ||2 represents the i-th pulse vector Δv i The amplitude of represents the state vector of the spacecraft in space, T represents the transpose; f1(X) and f2(X) represent the dynamic equations described by the circular restricted three-body problem dynamics model and the high-precision ephemeris dynamics model, respectively; Δv max Indicates the maximum limit of pulse amplitude; represents the starting state vector of the transfer trajectory; represents the terminal state vector of the transfer trajectory; Time labels indicating the integration tracks; and Respectively represent the spacecraft status before and after the i-th pulse; X c With X t They represent the discrete data sets of the starting orbit and the target orbit, that is, the discrete orbit data with continuous position and velocity with time tags; G(X i ) represents other state constraints.

[0028] Furthermore, the direct transfer two-layer optimization model after splitting and reconstruction in step 3 includes an upper layer optimization model and a lower layer optimization model; wherein the upper layer optimization model is expressed as follows:

[0029] The objective function is to minimize the fuel J1: Energy J2 minimization:

[0030] Constraints:

[0031] in, Δv means that the upper optimization model transmits the state vectors of the start and end points of the transfer trajectory to the lower optimization model, and the lower optimization model transmits the pulse and index back to the upper optimization model; J includes J1 and J2; Δv means the pulse vector;

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

[0033] The objective function is to minimize the fuel J1: Energy J2 minimization:

[0034] Constraints: ||Δv i ||2<=Δv max 、G(X i )<=0.

[0035] Furthermore, the specific steps of step four include: using a multiverse optimization algorithm with adjustment capability to solve the upper-level optimization model, and calculating the starting and end states and each pulse time; using the calculation results of solving the upper-level optimization model as known quantities, using a sequential quadratic programming algorithm to solve the pulse sequence of the lower-level optimization model, and calculating the minimum fuel consumption or minimum energy consumption; feeding back the minimum fuel consumption or minimum energy consumption to the upper-level optimization model, using a multiverse optimization algorithm with adjustment capability to reselect points, repeating iterations until the expected index or the maximum number of iterations is reached, and outputting the index iterative evolution curve and the transfer orbit curve.

[0036] Furthermore, the specific process of solving the upper optimization model using the multiverse optimization algorithm with adjustment capability is as follows:

[0037] Step 41: Initialize the multiverse matrix U; input the maximum number of iterations K and the expected index J d ; Initialize the minimum probability WEP of the wormhole existing in the universe min With maximum WEP max ; Initialize the upper and lower bounds ub and lb of each optimization variable; The multiverse matrix U is expressed as follows:

[0038]

[0039] Each row represents a universe, and the normalized expansion rate of each universe is set to NI (U a ), U a represents the a-th row of the matrix U; m is the number of defined multiverses;

[0040] Step 42: Input the data of each row of the multiverse matrix U into the lower optimization model to obtain the index function corresponding to each row: After sorting, the initial optimal index J is obtained best And the corresponding initial optimal universe U best ;

[0041] Step 43: Taking the expansion rate into account, move the objects in the universe as shown below:

[0042]

[0043] Where U a,b represents the ath row and bth column of matrix U; rand represents a random number between 0 and 1; U c,bWhen the random number is less than the expansion rate of the universe, a corresponding object in other universes is randomly selected;

[0044] Step 44. Without considering the expansion rate, each object in the universe realizes local changes and moves towards the current optimal universe, as shown in the following formula:

[0045]

[0046] Where U best represents the universe with the best current indicator function; rand1 and rand2 are random numbers in the range [0,1]; TDR represents the step size of the object moving towards the current optimal universe; WEP represents the probability of the existence of wormholes in the multiverse space;

[0047] Step 4 and 5: Update the wormhole existence probability WEP and the step length TDR moving toward the current optimal universe, as shown in the following formula:

[0048]

[0049] In the formula, k represents the current iteration number; λ is the adjustment factor used to adjust the moving step size; c1 and c2 are autonomous adjustment parameters, representing different degrees of exploration and exploitation, and the values ​​of c1 and c2 are adjusted according to the following strategy:

[0050] Strategy 1: c1: 2.55-(2.05 / K)k; c2: 1.25+(1 / K)k;

[0051] Strategy 2: c1: 2.55-(2.05 / K)k; c2: 0.5+2k 3 / K 3 ;

[0052] Strategy 3: c1: 2.5-2k 3 / K 3 ; c2: 1.25+(1 / K)k;

[0053] Strategy 4: c1: 2.5-2k 3 / K 3 ; c2: 0.5+2k 3 / K 3 ;

[0054] Step 46: Repeat step 42 to update the index corresponding to each row of the matrix universe; input the updated multiverse matrix U into step 43 for iterative optimization until the number of iterations is equal to the maximum number of iterations K or the optimization index is less than the expected index J d , output the optimal universe U at this time best and the optimal index J best .

[0055] According to another aspect of the present invention, a two-layer optimization system for direct transfer between special Earth-Moon orbits under an ephemeris model is proposed, comprising:

[0056] A data correction module is configured to construct a high-precision ephemeris dynamics model, and use the high-precision ephemeris dynamics model to correct the acquired initial values ​​of discrete orbit data to obtain corrected discrete orbit data with a time tag;

[0057] A model building module is configured to build a direct transfer optimization model based on the corrected discrete orbit data with time tags and the high-precision ephemeris dynamics model; split and reconstruct the direct transfer optimization model to obtain a split and reconstructed direct transfer dual-layer optimization model;

[0058] The optimization solution module is configured to iteratively solve the split and reconstructed direct transfer two-layer optimization model using a multiverse optimization algorithm with adjustment capability and a sequential quadratic programming algorithm until the expected index or the maximum number of iterations is reached, and the index evolution iteration curve and the transfer trajectory curve are obtained by solution.

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

[0060] In order to meet the transfer orbit design task requirements in lunar exploration and lunar space station construction, it is necessary to consider the transfer orbit design problem under higher precision. The transfer orbit design between special orbits of the Earth and the Moon under the ephemeris model will face challenges such as complex multi-body gravitational environment, strict energy constraints and flight window restrictions. Especially in actual aerospace engineering problems, the orbit transfer problem is not only an orbital dynamics problem, but also a spacecraft rendezvous problem that needs to consider time tags. The lunar probe needs to arrive at a certain position at a certain time to achieve rendezvous and docking with the lunar space station or another probe. For the direct transfer problem between any two special orbits of the Earth-Moon system under the high-precision ephemeris dynamics model, the transfer orbit optimization design is performed according to the method of the present invention.

[0061] The present invention proposes a double-layer optimization method and system for direct transfer between special earth-moon orbits under the ephemeris model, a modeling and calculation method for transfer orbit design directly under the ephemeris model, by constructing high-precision discrete orbit data with time tags, establishing a double-layer optimization model with minimum pulse as an indicator, and designing a double-layer optimization iterative solution algorithm. The purpose is to obtain a high-precision transfer orbit between special earth-moon orbits taking into account the rendezvous time window while consuming as little fuel as possible.

[0062] The high-precision ephemeris orbit data proposed in the present invention is used to construct a modeling method for direct transfer optimization problems. The discrete orbit data under the introduced high-precision ephemeris dynamics model is used as a constraint. This method not only ensures the accuracy of the initial and final states of the orbit, but also enhances the universality of the algorithm. It provides strong technical support for future long-term lunar exploration, manned lunar landing, and construction of a lunar space station. Through the optimized transfer orbit design, it will provide mission execution efficiency and feasibility.

[0063] The algorithm of the present invention is superior to the traditional algorithm, and proposes a combination of intelligent algorithm and numerical algorithm, combining the adaptive multiverse optimization algorithm AGMVO with adjustment capability with the sequential quadratic programming algorithm SQP, ensuring the convergence of the algorithm, improving the design accuracy and calculation efficiency of the transfer orbit, improving the efficiency of solving the double-layer optimization problem, and expanding the application scope of the algorithm. The present invention is based on the optimization model modeling method of discrete orbit data modeling and the orbit calculation method of the direct transfer double-layer optimization model based on the combination of the multiverse optimization algorithm MVO and the sequential quadratic programming algorithm SQP, which are used to solve the technical problem that the analytical method and intelligent method of the transfer orbit design between the Earth-Moon orbit under the traditional ephemeris model cannot be applied to the direct transfer between the Earth-Moon special orbits with higher precision requirements, and realizes the minimum pulse transfer orbit design that meets multiple constraints between various types of special orbits under the ephemeris model in the Earth-Moon system. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] 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.

[0065] Figure 1 It is a flow chart of a double-layer optimization method for direct transfer between special earth-moon orbits under an ephemeris model according to an embodiment of the present invention;

[0066] Figure 2 Schematic diagram of the rotation coordinate system of the center of mass of the earth-moon system in an embodiment of the present invention;

[0067] Figure 3 is a schematic diagram of the relationship between various coordinate systems in an embodiment of the present invention;

[0068] Figure 4 It is a schematic diagram of solving a double-layer optimization model based on a combination of an adjustable multiverse optimization algorithm AGMVO and a sequential quadratic programming algorithm SQP in an embodiment of the present invention;

[0069] Figure 5 Schematic diagram of iterative evolution curve of indicators in an embodiment of the present invention;

[0070] Figure 6 is an example diagram of a transfer orbit in a rotating coordinate system of the center of mass of the earth-moon system according to an embodiment of the present invention;

[0071] Figure 7 It is an example diagram of the transfer orbit in the J2000 coordinate system in an embodiment of the present invention. DETAILED DESCRIPTION

[0072] In order to enable those skilled in the art to better understand the scheme of the present invention, exemplary implementations or embodiments of the present invention will be described below in conjunction with the accompanying drawings. Obviously, the described implementations or embodiments are only implementations or embodiments of a part of the present invention, not all of them. Based on the implementations or embodiments of the present invention, all other implementations or embodiments obtained by ordinary technicians in the field without creative work should fall within the scope of protection of the present invention.

[0073] The embodiment of the present invention provides a double-layer optimization method for direct transfer between special earth-moon orbits under an ephemeris model, such as Figure 1 As shown, the method comprises the following steps:

[0074] Step 1: construct a high-precision ephemeris dynamics model, and use the high-precision ephemeris dynamics model to correct the acquired initial values ​​of discrete orbit data to obtain corrected discrete orbit data with time tags;

[0075] Step 2: construct a direct transfer optimization model based on the corrected discrete orbit data with time tags and the high-precision ephemeris dynamics model;

[0076] Step 3: split and reconstruct the direct transfer optimization model to obtain a direct transfer two-layer optimization model after split and reconstructing;

[0077] Step 4: Use a multiverse optimization algorithm with adjustment capability and a sequential quadratic programming algorithm to iteratively solve the split and reconstructed direct transfer two-layer optimization model until the expected index or the maximum number of iterations is reached, and obtain the index evolution iteration curve and the transfer trajectory curve.

[0078] The method starts with step one, in which a high-precision ephemeris dynamics model is constructed, and the acquired initial values ​​of discrete orbit data are corrected using the high-precision ephemeris dynamics model to obtain corrected discrete orbit data with time tags.

[0079] According to an embodiment of the present invention, firstly, a circular restricted three-body problem dynamics model is established in the Earth-Moon mass-center rotating coordinate system, and discrete data of the starting orbit and the target orbit under the circular restricted three-body problem dynamics model are calculated as orbit data. The above process is the prior art; then, a high-precision ephemeris dynamics model is established in the J2000 inertial coordinate system, and the orbit data under the circular restricted three-body problem dynamics model is corrected under the high-precision ephemeris dynamics model, and discrete orbit data with time tags that meet the requirements of ephemeris dynamics are calculated.

[0080] In order to facilitate the description of the transfer orbit between the special orbits of the Earth and the Moon, the relevant concepts are defined and the Earth-Moon mass-center rotating coordinate system, the Earth-centered J2000 inertial coordinate system, the Earth-centered fixed coordinate system and the Moon-centered inertial system are established as follows:

[0081] like Figure 2 As shown, the Earth-Moon mass center rotating coordinate system / rendezvous coordinate system: Under the three-body system, the rendezvous coordinate system B-xyz can be defined, whose coordinate origin is at the mass center of the three-body system, that is, the common mass center B of the Earth and the Moon. The x-axis points from the system mass center B to the mass center of the second main celestial body, that is, the moon. The xy coordinate plane is the plane of relative motion of the two main celestial bodies, and the x-axis, y-axis and z-axis form a right-handed coordinate system.

[0082] Geocentric J2000 inertial coordinate system: In a multi-body system, the J2000 inertial coordinate system OX can be defined E Y E Z E , the origin of the coordinate system O is located at the center of mass of the earth, Z E The axis points to the North Pole of the Earth (the direction of the celestial north pole), X E The axis points to the vernal equinox (one of the intersections of the ecliptic and the celestial equator, at the time of J2000.0), Y E The axes are determined according to the right-hand rule, forming a right-hand coordinate system.

[0083] Geocentric fixed coordinate system: In a multi-body system, the geocentric fixed coordinate system Ox can be defined E y E z E , the origin O of the coordinate system is located at the center of mass of the earth, z E The axis points to the North Pole of the Earth (the direction of the celestial north pole), E The axis points to the intersection of the prime meridian and the equator, y E The axes are determined according to the right-hand rule, forming a right-hand coordinate system.

[0084] Moon-centered inertial coordinate system: In a multi-body system, the moon-centered inertial coordinate system OX can be defined M Y M Z MThe origin O of the coordinate system is located at the mass center of the moon, and the coordinate axes in the three directions are parallel to the corresponding coordinate axes of the geocentric J2000 inertial coordinate system.

[0085] The relationship between the coordinate systems is as follows Figure 3 shown.

[0086] In order to conveniently characterize the special orbits in the Earth-Moon system, a circular restricted three-body problem dynamic model is established in the Earth-Moon mass center rotating coordinate system as follows:

[0087]

[0088] Where: x, y, z represent the position of the spacecraft in the three-body combined coordinate system; μ represents the mass ratio parameter of the two main celestial bodies; r1 and r2 represent the distances from the spacecraft to the two main celestial bodies, respectively.

[0089] The circular restricted three-body problem dynamics model has five equilibrium points, also called libration points. There are also many special periodic solutions around the libration points and the main celestial body, which constitute the special periodic orbits and quasi-periodic orbits in the circular restricted three-body problem dynamics model. The Halo orbit family and its special type of nearly linear halo orbit will have outstanding application value in the future, and are of great significance for lunar exploration and landing at the lunar north and south poles. The large-amplitude retrograde periodic orbit DRO has high stability, so it is suitable for the long-term on-orbit mission requirements of the lunar space station or probe.

[0090] When a special orbit type and condition requirements are given, a specific periodic orbit can be generated according to existing special orbit design methods, such as LP perturbation method, elliptic extension method, etc., and the orbit data under the circular restricted three-body problem dynamics model, i.e., discrete orbit data, is output; the discrete orbit data is the orbit data between the starting orbit and the target orbit, including position data, velocity data and time label; it can be represented by the following matrix:

[0091]

[0092] Among them, x i ,y i 、z i (i=1,2,3,...,n) represents the position data in three directions under the circular restricted three-body problem dynamics model; Represents the velocity data in three directions under the circular restricted three-body problem dynamics model.

[0093] The circular restricted three-body problem dynamics model is only suitable for the conceptual design stage because it assumes ideal conditions and ignores the eccentricity of the moon's orbit around the earth and other orbital perturbations. In the final stage of transfer orbit design, a more accurate ephemeris orbit dynamics model is required. Under the ephemeris model, the positions of different celestial bodies at each moment are obtained by looking up the ephemeris table, and the dynamics equations of the navigation vehicle need to be deduced in the inertial system.

[0094] The equations under the high-precision ephemeris dynamics model of the J2000 inertial coordinate system are as follows:

[0095]

[0096] In the formula, R EI =[X EI ,Y EI ,Z EI ] T represents the position vector of the spacecraft in the J2000 inertial coordinate system; μ E represents the gravitational constant of the Earth; + indicates the gravitational acceleration at the center of the earth; A N A is the perturbation acceleration caused by the gravity of other planets in the solar system outside the main celestial body; NSE A is the non-spherical perturbation acceleration of the main celestial body; R A is the solar pressure perturbation acceleration; D is the atmospheric drag perturbation acceleration.

[0097] The models of each perturbation term are as follows:

[0098] (1) Non-spherical perturbation of the main celestial body: The main celestial body is usually a non-uniform sphere with uneven density. The gravitational potential function is generally expressed by spherical harmonics. In the Earth-Moon system, the main celestial bodies include the Earth and the Moon. Both are usually non-uniform spheres with uneven density. Their gravitational potential functions are generally established using spherical harmonics in the Earth / Moon fixed system. Currently, multiple Earth-Moon gravitational field models have been established. Commonly used Earth gravitational field models include SAO-III, JGM-3, WGS-84, etc., and Moon gravitational field models include LP165 and GLGM-2, etc.;

[0099] (2) Perturbation acceleration caused by the gravity of other planets in the solar system outside the main celestial body: The gravity of other celestial bodies such as Mars and Mercury; in the Earth-Moon space, it mainly includes the gravity of the moon (in the moon-centered inertial system, the gravity of the earth is the perturbation gravity), the gravity of the sun, and the perturbation gravity of other planets in the solar system such as Venus and Jupiter, which are calculated as follows:

[0100]

[0101] In the formula, N represents the number of perturbed celestial bodies in the calculation, R i is the position vector of the spacecraft relative to the i-th perturbing celestial body, ρi is the position vector of the ith perturbed celestial body relative to the center of the earth. When designing the orbit of a spacecraft in the cis-lunar space, generally only the gravitational perturbations of the sun and the moon are considered. The position vectors of the sun and the moon relative to the center of the earth are obtained from the DE430 ephemeris published by JPL.

[0102] (3) Solar light pressure perturbation: Solar light pressure perturbation has an important influence on the long-term maintenance of the Earth-Moon space orbit and long-term orbit transfer. The general method for calculating solar light pressure is as follows:

[0103]

[0104] In the formula, C r is the reflection coefficient, S is the area of ​​the aircraft, m is the mass of the aircraft, K is the shadow coefficient (all obscured is 0, no obscuration is 1), is the light pressure flux at 1AU, c is the speed of light, r0 is 1AU, r s is the position vector from the sun to the current orbit point.

[0105] (4) Atmospheric drag perturbation acceleration:

[0106]

[0107] In the formula, C D is the drag coefficient, the nominal value can be 2.2; S / m is the ratio of the cross-section to the mass in the aerodynamic direction; is the atmospheric density; V D is the speed of the spacecraft relative to the atmosphere. If it is assumed that the atmosphere rotates synchronously with the Earth, then V D The calculation formula is:

[0108]

[0109] In the formula, R E and are the position and velocity vector of the spacecraft relative to the center of the earth; ω E is the angular velocity of the Earth's rotation.

[0110] Define R = [X, Y, Z] T and are the position vector and state vector of the detector in the inertial system, r = [x, y, z] T and are the position vector and state vector of the detector in the rotating system. Then the state vector of the rotating system is converted to the inertial system to satisfy the following relationship:

[0111]

[0112] in is the state transfer matrix from the rotating system to the inertial system, and its specific expression is:

[0113]

[0114] Then the transformation relationship from the rotating coordinate system state vector to the inertial coordinate system is:

[0115]

[0116] The orbital data under the circular restricted three-body model are transformed into the initial value of discrete orbital data under the high-precision ephemeris model through coordinate transformation; the coordinate transformation is the transformation of the state vector from the rotating system to the J2000 inertial coordinate system; therefore, the initial value of discrete orbital data under the high-precision ephemeris dynamics model can be expressed by the following matrix:

[0117]

[0118] Furthermore, the high-precision ephemeris dynamics model is used to correct the orbital data of the transfer orbit between the Earth and the Moon, and discrete orbital data with time tags are calculated. The specific process is as follows:

[0119] (1) Position continuous correction: The position differential correction equation is used to correct the initial value of the discrete orbit data obtained to obtain the position continuous discrete orbit data; the initial value of the discrete orbit data or the state quantity of each point of the previous iteration result is used as the starting point of the integration, and it is redefined as i=1,2,...,n-1, the end point of the integral is obtained after orbit recursion: i=1,2,…,n-1.

[0120] The position deviation of the i-th segment can be defined as δr if =[X if ,Y if ,Z if ] T -[X (i+1)f ,Y (i+1)f ,Z (i+1)f ] T , use the following differential correction equation to obtain the change in speed at point i The position differential correction equation is:

[0121] δv i0 =M r δr if

[0122] in, represents the integral end point obtained after orbit recursion, i=1,2,...,n-1; represents the starting point of orbital integration; δr if Indicates the deviation between the i-th point integral end position and the i+1-th point integral start position; δv i0Represents the change in velocity at the starting point of the i-th point integration; the partial derivative result is obtained from the state transfer matrix of the high-precision ephemeris dynamics model. Repeat this position correction operation for n-1 segments of the trajectory until each segment is positionally continuous. At this time, the discrete orbit data with positional continuity can be obtained as follows:

[0123]

[0124] in, i=1,2,...,n-1.

[0125] (2) Speed ​​Continuity Correction: After the first level differential correction of position continuous correction, the trajectory is continuous in position, but there will be a speed jump at each segment connection point. The trajectory after the position continuous correction has a speed jump at each segment connection point. The speed differential correction equation is used to correct the jump to obtain discrete orbit data with continuous speed; the speed deviation of each point can be defined as Use the following differential correction equation to obtain the change in position of point i: i0 =[δX i0 ,δY i0 ,δZ i0 ] T ; The speed differential correction equation is:

[0126]

[0127] in, δr j0 represents the change in the starting point of the j-th point, j = 2, ..., n-1; δv kf represents the deviation between the velocity at the end of the integration at the kth point and the velocity at the start of the integration at the k+1th point, k=1,2,...,n-1; r j0 represents the starting point of the j-th point integration; the partial derivative result is obtained from the state transfer matrix of the high-precision ephemeris dynamics model. At this time, the discrete orbit data with continuous position can be obtained as follows:

[0128]

[0129] in, i=2,...,n-1.

[0130] (3) After completing the speed correction, since the positions of the end points of each segment are also changed, the position correction must be continued. This is repeated iteratively until the position error and the speed error reach the preset threshold and are both within the allowable range. It can be considered that discrete orbit data with continuous position and speed are obtained, and the correction is completed.

[0131] In addition, it is necessary to add a time tag to each discrete point data of the starting orbit and the target orbit; after obtaining the discrete orbit data with continuous corrected position and velocity, add a time tag to the discrete orbit data with continuous corrected position and velocity, and the time tag is to add the real ephemeris epoch as the time tag according to the integration step, so as to obtain the discrete orbit data with continuous position and velocity with the time tag; the discrete orbit data is:

[0132]

[0133] In the formula, 0 X i , 0 Y i , 0 Z i , f X i , f Y i , f Z i (i=1,2,3,...,n) respectively represent the position data of the starting orbit and the target orbit in three directions of the high-precision ephemeris dynamics model; Respectively represent the velocity data of the starting orbit and the target orbit in three directions of the high-precision ephemeris dynamics model; 0 t i , f t i (i=1,2,3,...,n) represent the time labels of the starting orbit and the target orbit respectively, generally in the form of Julian day.

[0134] Then, step 2 is performed, in which a direct transfer optimization model is constructed based on the corrected discrete orbit data and the high-precision ephemeris dynamics model.

[0135] According to an embodiment of the present invention, the direct transfer optimization model includes an objective function and a constraint model; the constraint model includes: establishing a dynamic model constraint based on the dynamics considering the ephemeris model; establishing an initial and final state constraint model based on discrete orbit data under a high-precision ephemeris dynamics model; establishing other constraint models such as the continuity of the position before and after the velocity pulse, the upper limit of the velocity pulse, and the pulse point limit; and establishing a minimum pulse index function model for energy optimization or fuel optimization.

[0136] Establish the initial optimization problem of the direct transfer optimization model between the special Earth-Moon orbit, that is, establish the constraints and objective function.

[0137] (1) Dynamic model constraints

[0138] The circular restricted three-body problem dynamics model and the high-precision ephemeris dynamics model are transformed and can be written as follows:

[0139]

[0140] Where: The state vector describing the spacecraft in space, f1(X) and f2(X) are expressed as a circular restricted three-body problem dynamics model and a high-precision ephemeris dynamics model, respectively.

[0141] When there is no pulse effect, the spacecraft should always satisfy the dynamic equation, which can be expressed in integral form as follows:

[0142]

[0143] (2) Initial and final state constraints

[0144] The initial and final state constraints require that the starting point of the transfer track is located on the starting track and the end point is located on the target track. The data table consisting of the starting track and the target track based on the discrete track data is given, and the feasible number set of the starting point and the feasible number set of the end point can be defined as follows:

[0145]

[0146] Where p and q are the number of discrete points of orbital data.

[0147] In order to ensure the accuracy of orbital discrete points and the efficiency of integral calculation, orbital integral or recursive algorithm generally adopts the variable order method for solving non-rigid differential equations, that is, the integral function ode113 in MATLAB, which can only design the generation of discrete points by changing the time of the required discrete points. The final generated discrete points may be too dense or too sparse in some parts of the given orbital parameter table due to unreasonable time series design, which will cause the search burden on the transfer orbit generation algorithm or miss the better solution. Therefore, two related methods for screening discrete points are given below:

[0148] (1) Angle averaging sequence: taking the starting point as the reference, checking the angle of the position vector between each discrete point and the reference point, and interpolating, and finally selecting points with an appropriate interval angle step (generally 1° interval), and finally forming a set of several discrete data;

[0149] (2) Accuracy averaging sequence: First, the starting point is used as the reference point. The state deviation between the next point and the reference point is checked by comprehensively considering the position and velocity. If the deviation between the two points is within the accuracy requirement, it is considered that the state difference between the two points is not large. This point is eliminated until a point is found whose state deviation with the reference point is greater than the accuracy requirement. In other words, it is considered that this point is much different from the reference point and will affect the transfer trajectory. This point is retained and considered as the reference point. The above operation is repeated until all points are screened. It can be expressed mathematically as follows:

[0150]

[0151] Where i is the number of state rows / columns in the original data table, n is the number of discrete points in the original data; j is the number of state rows / columns in the filtered data table; σ is the given state deviation accuracy.

[0152] Define X cj With X tl They correspond to discrete orbital data, and each point should contain time, position and speed, a total of 7 dimensions, namely:

[0153]

[0154] The transfer orbit starting point and orbit end point are X0 and X n , that is, the spacecraft state corresponding to the first pulse time and the last pulse time. The initial and final state constraints require that the state of the spacecraft before the pulse is applied at the initial time is a point in the initial orbit data, and the state after the pulse is applied at the terminal time is a point in the target orbit data. The constraints can be expressed as follows:

[0155]

[0156] (3) Other constraints

[0157] The position of each pulse moment is required to be continuous, and the pulse before and after can be expressed as follows:

[0158]

[0159] In addition, there may be other restrictions on the pulse point, such as the position of the pulse point, etc. It can be written in the following form:

[0160] G(X i )<=0

[0161] Considering the pulse upper limit, it can be expressed as the following constraint:

[0162] ||Δv i ||2≤Δv max

[0163] The indicators are expressed in fuel-optimal or energy-optimal form as follows:

[0164] or

[0165] In summary, the minimum pulse transfer problem (i.e., direct transfer optimization model) can be established as follows:

[0166] The objective function is to minimize the fuel J1: Energy J2 minimization:

[0167]

[0168] Where J1 represents the optimal fuel; J2 represents the optimal energy; n is the number of discrete points of the original data; N is the total number of pulses; ||Δv i ||2 represents the i-th pulse vector Δv i The amplitude of; i represents the number of pulses; Describes the state vector of the spacecraft in space, where T represents the transpose of the state vector of the spacecraft in space; f1(X) and f2(X) represent the dynamic equations described by the circular restricted three-body problem dynamics model and the high-precision ephemeris dynamics model, respectively; Δv i represents the i-th pulse vector, Δv max Indicates the maximum limit of pulse amplitude; represents the starting state vector of the transfer trajectory; represents the terminal state vector of the transfer trajectory; Indicates the time labels of the integration orbit; and Respectively represent the spacecraft status before and after the i-th pulse; X c With X t They represent the discrete data sets of the starting orbit and the target orbit, that is, the discrete orbit data with continuous position and velocity with time tags; G(X i ) represents other state constraints.

[0169] Then, step three is executed. In step three, the direct transfer optimization model is split and reconstructed to obtain a direct transfer two-layer optimization model after split and reconstructed.

[0170] According to an embodiment of the present invention, a specific method of constructing a dual-layer optimization model for direct transfer after splitting between special orbits of the Earth-Moon system is as follows:

[0171] All optimization variables are divided into two parts: the states of the starting point and the end point are upper-level variables, and the pulse sequence is the lower-level variable. Correspondingly, the constraints of the original optimization problem can also be divided into two parts, thus reconstructing the direct transfer two-level optimization model as shown below. The upper-level optimization model is expressed as:

[0172] The objective function is to minimize the fuel J1: Energy J2 minimization:

[0173]

[0174] in, Δv means that the upper optimization model transmits the state vectors of the start and end points of the transfer trajectory to the lower optimization model, and the lower optimization model transmits the pulse and index back to the upper optimization model; J includes J1 and J2; Δv means the pulse vector;

[0175] The lower optimization model is expressed as:

[0176] The objective function is to minimize the fuel J1: Energy J2 minimization:

[0177] ||Δv i ||2<=Δv max 、G(X i )<=0

[0178] The dual-level optimization essentially splits the minimum pulse direct transfer problem into two steps:

[0179] The first step is to select points on the starting orbit and the target orbit, and calculate the transfer time between the two points based on the time tags of these two points, that is, the optimization variables of the upper optimization model are the starting parking orbit parameters, the epoch parameters of the target special orbit, and the time of each pulse; the second step is to establish an optimal control problem model, and only optimize the control quantity, that is, optimize the pulse sequence, when the starting point, the end point and the transfer time are determined, that is, the optimization variables of the lower optimization model are each pulse sequence, and the optimization index is fuel optimization or energy optimization.

[0180] Then execute step four, in which the split and reconstructed direct transfer two-layer optimization model is iteratively solved using a multiverse optimization algorithm with adjustment capabilities and a sequential quadratic programming algorithm until the expected index or the maximum number of iterations is reached, and the index evolution iteration curve and the transfer trajectory curve are obtained. The specific steps include: using a multiverse optimization algorithm with adjustment capabilities to solve the upper-level optimization model, and calculating the starting and end states and each pulse time; using the calculation results of solving the upper-level optimization model as known quantities, using a sequential quadratic programming algorithm to solve the pulse sequence of the lower-level optimization model, and calculating the minimum fuel consumption or minimum energy consumption; feeding the minimum fuel consumption or minimum energy consumption back to the upper-level optimization model, using a multiverse optimization algorithm with adjustment capabilities to reselect points, repeating the iteration until the expected index or the maximum number of iterations is reached, and outputting the index iteration evolution curve and the transfer trajectory curve.

[0181] According to an embodiment of the present invention, an iterative solution of a direct transfer two-layer optimization algorithm based on a combination of a multiverse optimization algorithm AGMVO with an adjustable capability and a sequential quadratic programming algorithm SQP is performed until the expected index or the maximum number of iterations is reached, and an index iteration curve, a transfer orbit under a circular restricted three-body problem dynamics model, and a transfer orbit under a high-precision ephemeris dynamics model are obtained; specifically, the method includes: solving the point selection problem of an upper optimization model in the direct transfer two-layer optimization model based on the multiverse optimization algorithm AGMVO with an adjustable capability to obtain a selected starting point and a target point; solving the pulse solution problem of a lower optimization model in the direct transfer two-layer optimization model based on the sequential quadratic programming algorithm SQP, taking the result of the upper optimization model as a known quantity, solving the corresponding pulse sequence, and calculating the corresponding index; feeding the index back to the upper optimization model, reselecting the point using the optimization mechanism of the multiverse optimization algorithm AGMVO with an adjustable capability, repeating the iteration until the expected index or the maximum number of iterations is reached, and outputting the index evolution iteration curve, the transfer orbit under the three-body model, and the transfer orbit under the ephemeris model.

[0182] The upper-level optimization model solution method uses intelligent algorithms such as the multiverse optimization algorithm MVO with autonomous capabilities, and the lower-level optimization model solution method uses numerical algorithms such as the sequential quadratic programming algorithm SQP.

[0183] Since there is a fixed mapping relationship between the number of rows / columns and the data in this row / column, the mapping relationship between different rows / columns is different. Each row / column number and its corresponding data can be described as a specific unknown function:

[0184] X cj =F cj (j), X tl =F tl (l)

[0185] For analytical iterative methods, functions that cannot be explicitly expressed are often not conducive to solving. However, thanks to the universality of intelligent algorithms, data can be obtained by looking up tables, and dimensionality reduction can be achieved by finding the number of rows or columns: the original 12-dimensional variables of the starting and end states to be optimized are reduced to 2 dimensions through the natural constraints of the table data, that is, a row or column in the starting orbit data table and a row or column in the target orbit data table (for a family of orbital data, constraints can be imposed through a three-dimensional table, and the optimization variables are reduced to 4 dimensions). Therefore, the upper-level optimization model can be reconstructed as follows:

[0186] minimize or

[0187]

[0188] Wherein, j represents the position of the starting point in the starting track data, l represents the position of the target point in the target track data, p represents the total number of rows of the starting track data, and q represents the total number of rows of the target track data.

[0189] Since the relationship between the optimization variables and optimization indicators of the upper optimization model is unclear and the initial value of the optimization variable is difficult to design, an intelligent algorithm is selected for solution. Taking into account the computational efficiency and convergence, the multiverse optimization algorithm MVO with adjustment capability is selected for optimization.

[0190] In the traditional multiverse optimization algorithm MVO, all universes are the same in terms of the effects of local search and global search. Therefore, all universes can be considered as a group with the same strategy, which greatly reduces the optimal search ability and makes the algorithm more likely to fall into the local optimal solution. The position and velocity updates of each particle can be expressed by the following formula; for the upper-level optimization model, the multiverse optimization algorithm (Multi-Verse Optimizer, MVO) is selected for improvement, and an adaptive multiverse optimization algorithm (Adaptive Genetic Multi-Verse Optimizer, AGMVO) is proposed, which improves the algorithm's optimization ability and enhances the algorithm's ability to jump out of the local optimum. The present invention uses AGMVO to solve the upper-level optimization model, and the optimization variables are the starting orbit data number j, the target orbit data number l, and the transfer time t; initialize the multiverse matrix U; input the maximum number of iterations K and the expected index J d ; Initialize the minimum probability WEP of the wormhole existing in the universe min With maximum WEP max ; Initialize the upper and lower bounds ub and lb of each optimization variable; The specific process is as follows:

[0191] (1) Input the maximum number of iterations K and the expected index J d ; Initialize the minimum probability WEP of the wormhole existing in the universe min With maximum WEP max ; Initialize the upper and lower bounds ub and lb of each optimization variable; define the number of multiverses m, and initialize the multiverse matrix, that is, the multiverse matrix U is expressed as follows:

[0192]

[0193] Each row represents a universe, and the normalized expansion rate of each universe is set to NI (U a ), U a represents the a-th row of the matrix U; m is the number of defined multiverses.

[0194] (2) Input the data of each row of the multiverse matrix U into the lower optimization model to obtain the indicator function corresponding to each row:

[0195]

[0196] After sorting, the optimal index J is initialized best and its corresponding initialized optimal universe U best .

[0197] (3) Taking the expansion rate into account, the movement of objects in the universe is as follows:

[0198]

[0199] Among them, U a,b represents the ath row and bth column of matrix U; rand represents a random number between 0 and 1; U c,b When the random number is less than the expansion rate of the universe, a corresponding object from other universes is randomly selected. Through this movement mechanism, objects are constantly exchanged between multiverses without disturbance.

[0200] (4) Without considering the expansion rate, each object in the universe undergoes local changes and moves toward the current optimal universe, as shown in the following formula:

[0201]

[0202] Among them, U best represents the universe with the best current indicator function; lb and ub represent the lower and upper limits of the optimization variables respectively; rand1 and rand2 are random numbers in the range of [0.1]; WEP represents the probability of the existence of wormholes in the multiverse space, WEP min Indicates the minimum probability of setting, WEP max represents the maximum probability of the setting; TDR represents the step size of the object moving towards the current optimal universe.

[0203] (5) Update the wormhole existence probability WEP and the step length TDR of moving toward the current optimal universe, as shown in the following formula:

[0204]

[0205] In the formula, λ is the adjustment factor used to adjust the moving step size; K represents the maximum number of iterations, k represents the current number of iterations; c1 and c2 are autonomous adjustment parameters, representing different degrees of exploration and exploitation. In particular, the present invention uses different autonomous groups with common goals to enable the algorithm to have more random or directional searches. The core idea of ​​the multiverse optimization algorithm AGMVO with adjustment capability is to use different strategies to update c1 and c2, that is, to adjust the values ​​of c1 and c2. Four groups can be designed, that is, four c1 and c2 change strategies can be designed. The update function is shown in Table 1.

[0206] Table 1 Four update strategies for c1 and c2

[0207]

[0208] (6) Repeat the operation in step (2) to update the index corresponding to each row of the matrix universe; input the updated matrix U into step (3) for iterative optimization until the number of iterations is equal to the maximum number of iterations K or the optimization index is less than the expected index J d , output the optimal universe U at this time best and the optimal index J best , the optimization is completed.

[0209] After the splitting process, the free variables of the lower optimization model are only the 3(n+1)-dimensional pulse sequence Δv. According to the starting and ending states and each pulse moment obtained by the upper optimization model, the optimal control problem can be established as follows:

[0210] minimize or

[0211] ||Δv i ||2<=Δv max

[0212] After the upper optimization model obtains the initial and final states and the transfer time, it is still a nonlinear optimization problem with constraints, but compared with the original problem, the number of variables is greatly reduced and the difficulty of solving is greatly reduced. When the transfer process is only two pulses, the lower optimization model is reconstructed as a Lambert problem under the ephemeris model, which can be solved by algorithms such as the Newton-homotopy method, quasi-linearization-local variation iteration method, etc. The present invention solves the minimum pulse through the Sequential Quadratic Programming (SQP) algorithm.

[0213] According to the pulse calculation index solved by the sequential quadratic programming algorithm SQP, the index is fed back to the upper optimization model, and the point is re-selected using the AGMVO optimization mechanism with adjustable multiverse optimization algorithm. The direct transfer two-layer optimization model is subjected to an iterative solution process based on the direct transfer two-layer optimization algorithm combined with the AGMVO with adjustable multiverse optimization algorithm and the sequential quadratic programming algorithm SQP, such as Figure 4 As shown, the iteration is repeated until the expected index or the maximum number of iterations is reached, and the index iteration evolution curve and the transfer trajectory curve are output.

[0214] Another embodiment of the present invention provides a two-layer optimization system for direct transfer between special Earth-Moon orbits under an ephemeris model, the system comprising:

[0215] A data correction module is configured to construct a high-precision ephemeris dynamics model, and use the high-precision ephemeris dynamics model to correct the acquired initial values ​​of discrete orbit data to obtain corrected discrete orbit data with a time tag;

[0216] A model building module is configured to build a direct transfer optimization model based on the corrected discrete orbit data with time tags and the high-precision ephemeris dynamics model; split and reconstruct the direct transfer optimization model to obtain a split and reconstructed direct transfer dual-layer optimization model;

[0217] The optimization solution module is configured to iteratively solve the split and reconstructed direct transfer two-layer optimization model using a multiverse optimization algorithm with adjustment capability and a sequential quadratic programming algorithm until the expected index or the maximum number of iterations is reached, and the index evolution iteration curve and the transfer trajectory curve are obtained by solution.

[0218] In this embodiment, preferably, the specific process of solving the upper optimization model by the multiverse optimization algorithm with adjustment capability is as follows:

[0219] Step 41: Initialize the multiverse matrix U; input the maximum number of iterations K and the expected index J d ; Initialize the minimum probability WEP of the wormhole existing in the universe min With maximum WEP max ; Initialize the upper and lower bounds ub and lb of each optimization variable; The multiverse matrix U is expressed as follows:

[0220]

[0221] Each row represents a universe, and the normalized expansion rate of each universe is set to NI (U a ), U a represents the a-th row of the matrix U; m is the number of defined multiverses;

[0222] Step 42: Input the data of each row of the multiverse matrix U into the lower optimization model to obtain the index function corresponding to each row: After sorting, the initial optimal index J is obtained best And the corresponding initial optimal universe U best ;

[0223] Step 43: Taking the expansion rate into account, move the objects in the universe as shown below:

[0224]

[0225] Where U a,b represents the ath row and bth column of matrix U; rand represents a random number between 0 and 1; U c,bWhen the random number is less than the expansion rate of the universe, a corresponding object in other universes is randomly selected;

[0226] Step 44. Without considering the expansion rate, each object in the universe realizes local changes and moves towards the current optimal universe, as shown in the following formula:

[0227]

[0228] Where U best represents the universe with the best current indicator function; rand1 and rand2 are random numbers in the range [0,1]; TDR represents the step size of the object moving towards the current optimal universe; WEP represents the probability of the existence of wormholes in the multiverse space;

[0229] Step 4 and 5: Update the wormhole existence probability WEP and the step length TDR moving toward the current optimal universe, as shown in the following formula:

[0230]

[0231] In the formula, k represents the current iteration number; λ is the adjustment factor used to adjust the moving step size; c1 and c2 are autonomous adjustment parameters, representing different degrees of exploration and exploitation, and the values ​​of c1 and c2 are adjusted according to the following strategy:

[0232] Strategy 1: c1: 2.55-(2.05 / K)k; c2: 1.25+(1 / K)k;

[0233] Strategy 2: c1: 2.55-(2.05 / K)k; c2: 0.5+2k 3 / K 3 ;

[0234] Strategy 3: c1: 2.5-2k 3 / K 3 ; c2: 1.25+(1 / K)k;

[0235] Strategy 4: c1: 2.5-2k 3 / K 3 ; c2: 0.5+2k 3 / K 3 ;

[0236] Step 46: Repeat step 42 to update the index corresponding to each row of the matrix universe; input the updated multiverse matrix U into step 43 for iterative optimization until the number of iterations is equal to the maximum number of iterations K or the optimization index is less than the expected index J d , output the optimal universe U at this time best and the optimal index J best .

[0237] The functions of a double-layer optimization system for direct transfer between special Earth-Moon orbits under an ephemeris model in an embodiment of the present invention can be described by the double-layer optimization method for direct transfer between special Earth-Moon orbits under an ephemeris model described above. 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.

[0238] The technical effect of the present invention is further verified through experiments.

[0239] The transfer between the NRHO orbiting the L1 point and the Halo orbiting the L2 point in the Earth-Moon system under the ephemeris model is used as an experimental example to simulate and verify the technical effect of the present invention. The orbit transfer mission is designed as follows: the lunar space station is deployed in the north-bound NRHO orbit around the L1 point. After the cargo cabin is replenished with the space station, it is expected to transfer to the north-bound Halo orbit around the L2 point and rendezvous and dock with the orbiter on the Halo orbit. The present invention is suitable for multi-pulse transfer between orbits. This example takes the direct transfer of two pulses as an example, and the relevant parameters are selected as shown in Table 2 below:

[0240] Table 2 Simulation parameters

[0241]

[0242] The maximum number of evolutionary iterations is selected as 50, and modeling, transformation and solution are performed according to the method described in the embodiment of the present invention, and finally the index evolution iteration curve under the maximum number of iterations of 50, the transfer orbit in the center of mass rotation coordinate system and the transfer orbit in the J2000 coordinate system are obtained, as shown in FIG. Figure 5-Figure 7 shown.

[0243] It should be noted that the transfer trajectory curve obtained by the present invention is a pre-designed trajectory of the spacecraft before the actual flight mission, which can be used as the nominal trajectory curve in actual control. When used specifically, it is necessary to make mid-course corrections according to the actual flight situation. By controlling based on the results of the embodiment of the present invention as the standard, a theoretically optimal fuel flight trajectory can be obtained.

[0244] Although the present invention has been described according to a limited number of embodiments, it will be apparent to those skilled in the art, with the benefit of the above description, that other embodiments are contemplated within the scope of the invention thus described. The disclosure of the present invention is intended to be illustrative rather than restrictive of the scope of the invention, which is defined by the appended claims.

[0245] The present invention cites the following documents:

[0246] [1]WANG Y, ZHANG R, ZHANG C, et al.Transfers between NRHOs and DROs in the Earth-Moon system[J]. Acta Astronautica, 2021, 186: 60-73.

[0247] [2]ZIMOVAN-SPREEN EM, HOWELL KC, DAVISD C.Dynamical Structures NearbyNRHOs with Applications to Transfer Design in Cislunar Space[J].The Journal of the Astronautical Sciences, 2022,69(3):718-744.

[0248] [3] Lu Lin, Li Haiyang, Zhou Wanmeng, et al. Design of Earth-Moon transfer orbit for manned lunar polar exploration[J]. Acta Astronautics, 2021, 42(4): 425-436.

[0249] [4] Peng Qibo, He Xingji, Chen Tianji, et al. Orbit design for the full flight process of Earth-Moon free return based on the ephemeris model[J]. Acta Astronautics, 2023, 44(1): 43-51.

Claims

1. A two-layer optimization method for direct transfer between special Earth-Moon orbits under the ephemeris model, characterized in that: The following steps are involved: Step 1: construct a high-precision ephemeris dynamics model, and use the high-precision ephemeris dynamics model to correct the acquired initial values ​​of discrete orbit data to obtain corrected discrete orbit data with time tags; Step 2: construct a direct transfer optimization model based on the corrected discrete orbit data with time tags and the high-precision ephemeris dynamics model; Step 3: split and reconstruct the direct transfer optimization model to obtain a direct transfer two-layer optimization model after split and reconstructing; Step 4: Use a multiverse optimization algorithm with adjustment capability and a sequential quadratic programming algorithm to iteratively solve the split and reconstructed direct transfer two-layer optimization model until the expected index or the maximum number of iterations is reached, and obtain the index evolution iteration curve and the transfer trajectory curve.

2. The double-layer optimization method for direct transfer between special Earth-Moon orbits under an ephemeris model according to claim 1 is characterized in that: The discrete track data in step 1 is the track data between the starting track and the target track, including position data, speed data and time tags.

3. The double-layer optimization method for direct transfer between special Earth-Moon orbits under an ephemeris model according to claim 2 is characterized in that: The high-precision ephemeris dynamics model described in step 1 is established as follows: In the formula, R EI =[X EI ,Y EI ,Z EI ] T represents the position vector of the spacecraft in the J2000 inertial coordinate system; μ E represents the gravitational constant of the Earth; + indicates the gravitational acceleration at the center of the earth; A N A is the perturbation acceleration caused by the gravity of other planets in the solar system outside the main celestial body; NSE A is the non-spherical perturbation acceleration of the main celestial body; R A is the solar pressure perturbation acceleration; D is the atmospheric drag perturbation acceleration.

4. The double-layer optimization method for direct transfer between special Earth-Moon orbits under an ephemeris model according to claim 3 is characterized in that: In step 1, the high-precision ephemeris dynamics model is used to correct the initial value of the acquired discrete orbit data, including continuous position correction and continuous velocity correction; wherein, The continuous position correction is: using the position differential correction equation to correct the initial value of the discrete orbit data obtained to obtain discrete orbit data with continuous position; the position differential correction equation is: dv i0 =M r δr if In the formula, represents the integral end point obtained after orbit recursion, i=1,2,...,n-1; represents the starting point of orbital integration; δr if Indicates the deviation between the i-th point integral end position and the i+1-th point integral start position; δv i0 Indicates the change in velocity at the starting point of the i-th integration point; The speed continuous correction is: the trajectory after the position continuous correction has a speed jump at each segment connection point, and the jump is corrected by using the speed differential correction equation to obtain discrete track data with continuous speed; the speed differential correction equation is: In the formula, δr j0 represents the change in the starting point of the j-th point, j = 2, ..., n-1; δv kf represents the deviation between the velocity at the end of the integration at the kth point and the velocity at the start of the integration at the k+1th point, k=1,2,...,n-1; r j0 Indicates the starting point of the j-th point integration; After the speed continuous correction is completed, the position of each end point is changed, and the position correction is continued; the iterative correction is repeated until the position error and the speed error reach the preset threshold value, and the discrete track data with continuous position and speed after correction is obtained.

5. The double-layer optimization method for direct transfer between special Earth-Moon orbits under the ephemeris model according to claim 4 is characterized in that: Step one also includes: after obtaining the discrete orbit data with continuous corrected position and velocity, adding a time tag to the discrete orbit data with continuous corrected position and velocity, wherein the time tag is a time tag that adds a real ephemeris epoch according to the integration step, thereby obtaining the discrete orbit data with continuous position and velocity with the time tag.

6. The double-layer optimization method for direct transfer between special Earth-Moon orbits under the ephemeris model according to claim 5 is characterized in that: The direct transfer optimization model described in step 2 is: The objective function is to minimize the fuel J1: Energy J2 minimization: Constraints: Where N is the total number of pulses; ||Δv i ||2 represents the i-th pulse vector Δv i The amplitude of represents the state vector of the spacecraft in space, T represents the transpose; f1(X) and f2(X) represent the dynamic equations described by the circular restricted three-body problem dynamics model and the high-precision ephemeris dynamics model, respectively; Δv max Indicates the maximum limit of pulse amplitude; represents the starting state vector of the transfer trajectory; represents the terminal state vector of the transfer trajectory; Time labels representing the integration tracks; and Respectively represent the spacecraft status before and after the i-th pulse; X c With X t They represent the discrete data sets of the starting orbit and the target orbit, that is, the discrete orbit data with continuous position and velocity with time tags; G(X i ) represents other state constraints.

7. The double-layer optimization method for direct transfer between special Earth-Moon orbits under the ephemeris model according to claim 6 is characterized in that: The direct transfer two-layer optimization model after splitting and reconstruction in step 3 includes an upper layer optimization model and a lower layer optimization model; wherein the upper layer optimization model is expressed as follows: The objective function is to minimize the fuel J1: Energy J2 minimization: Constraints: in, It means that the upper optimization model transmits the state vectors of the starting point and the end point of the transfer trajectory to the lower optimization model, and the lower optimization model transmits the pulse and index back to the upper optimization model; J includes J1 and J2; Δv represents the pulse vector; The lower layer optimization model is expressed as follows: The objective function is to minimize the fuel J1: Energy J2 minimization: Constraints:

8. The double-layer optimization method for direct transfer between special Earth-Moon orbits under the ephemeris model according to claim 7 is characterized in that: The specific steps of step four include: using a multiverse optimization algorithm with adjustment capability to solve the upper optimization model, and calculating the starting and end states and each pulse time; using the calculation result of solving the upper optimization model as a known quantity, using a sequential quadratic programming algorithm to solve the pulse sequence of the lower optimization model, and calculating the minimum fuel consumption or minimum energy consumption; feeding back the minimum fuel consumption or minimum energy consumption to the upper optimization model, using a multiverse optimization algorithm with adjustment capability to reselect points, repeating iterations until the expected index or the maximum number of iterations is reached, and outputting the index iterative evolution curve and the transfer orbit curve.

9. The double-layer optimization method for direct transfer between special Earth-Moon orbits under the ephemeris model according to claim 8 is characterized in that: The specific process of solving the upper optimization model using the multiverse optimization algorithm with adjustment capability is as follows: Step 41: Initialize the multiverse matrix U; input the maximum number of iterations K and the expected index J d ; Initialize the minimum probability WEP of the wormhole existing in the universe min With maximum WEP max ; Initialize the upper and lower bounds ub and lb of each optimization variable; The multiverse matrix U is expressed as follows: Each row represents a universe, and the normalized expansion rate of each universe is set to NI (U a ), U a represents the a-th row of the matrix U; m is the number of defined multiverses; Step 42: Input the data of each row of the multiverse matrix U into the lower optimization model to obtain the index function corresponding to each row: After sorting, the initial optimal index J is obtained best And the corresponding initial optimal universe U best ; Step 43: Taking the expansion rate into account, move the objects in the universe as shown below: Where U a,b represents the ath row and bth column of matrix U; rand represents a random number between 0 and 1; U c,b When the random number is less than the expansion rate of the universe, a corresponding object in other universes is randomly selected; Step 44. Without considering the expansion rate, each object in the universe realizes local changes and moves towards the current optimal universe, as shown in the following formula: Where U best represents the universe with the best current indicator function; rand1 and rand2 are random numbers in the range [0,1]; TDR represents the step size of the object moving towards the current optimal universe; WEP represents the probability of the existence of wormholes in the multiverse space; Step 4 and 5: Update the wormhole existence probability WEP and the step length TDR moving toward the current optimal universe, as shown in the following formula: In the formula, k represents the current iteration number; λ is the adjustment factor used to adjust the moving step size; c1 and c2 are autonomous adjustment parameters, representing different degrees of exploration and exploitation, and the values ​​of c1 and c2 are adjusted according to the following strategy: Strategy 1: c1: 2.55-(2.05 / K)k; c2: 1.25+(1 / K)k; Strategy 2: c1: 2.55-(2.05 / K)k; c2: 0.5+2k 3 / K 3 ; Strategy 3: c1: 2.5-2k 3 / K 3 ; c2: 1.25+(1 / K)k; Strategy 4: c1: 2.5-2k 3 / K 3 ; c2: 0.5+2k 3 / K 3 ; Step 46: Repeat step 42 to update the index corresponding to each row of the matrix universe; input the updated multiverse matrix U into step 43 for iterative optimization until the number of iterations is equal to the maximum number of iterations K or the optimization index is less than the expected index J d , output the optimal universe U at this time best and the optimal index J best .

10. A two-layer optimization system for direct transfer between special Earth-Moon orbits under the ephemeris model, characterized in that: include: A data correction module is configured to construct a high-precision ephemeris dynamics model, and use the high-precision ephemeris dynamics model to correct the acquired initial values ​​of discrete orbit data to obtain corrected discrete orbit data with a time tag; A model building module is configured to build a direct transfer optimization model based on the corrected discrete orbit data with time tags and the high-precision ephemeris dynamics model; split and reconstruct the direct transfer optimization model to obtain a split and reconstructed direct transfer dual-layer optimization model; The optimization solution module is configured to iteratively solve the split and reconstructed direct transfer two-layer optimization model using a multiverse optimization algorithm with adjustment capability and a sequential quadratic programming algorithm until the expected index or the maximum number of iterations is reached, and the index evolution iteration curve and the transfer trajectory curve are obtained by solution.

Images