Optimization method for multi-satellite optimal orbit in a planetary system

Abstract

The application discloses a kind of multi-satellite optimal orbit optimization method in planetary system inside borrowing force fuel consumption, belong to aerospace technical field.The application implementation method is: selecting the sequence of borrowing force satellite satisfying detection demand, with each transfer orbit period as variable, only adding maneuver when leaving initial orbit and being captured or entering target orbit by target satellite, the optimization target is minimum total velocity increment, and the orbit optimization problem is constructed under planetary mass center inertia system.According to the state of spacecraft in the parameter solving of transfer orbit between borrowing force satellites, the state of next satellite borrowing force is solved, the parameter of transfer orbit after borrowing force can be matched from the state to realize continuous solving, until the state when reaching target satellite is obtained.The optimization problem is solved using global optimization algorithm, the relative velocity of reaching each borrowing force satellite and the transfer orbit parameters before and after borrowing force are obtained, and the transfer of each segment transfer orbit under reference ephemeris is realized according to the transfer orbit parameters.

Description

Technical Field

[0001] This invention relates to a spacecraft orbit optimization method, and more particularly to a multi-satellite orbit optimization method within a planetary system, belonging to the field of aerospace technology. Background Technology

[0002] Multi-satellite levered orbits are a crucial approach for achieving low-energy transfers within target planetary systems during deep space exploration missions. Utilizing the levered cooperation between different satellites can effectively reduce spacecraft burnup during the transfer process. However, due to the large number of satellites and their short orbital periods within planetary systems, the levered orbit sequences are diverse, resulting in a significantly increased number of decision variables compared to interplanetary transfer levered orbits. This makes it difficult to solve using conventional orbit optimization methods. Therefore, simplifying the calculation model for levered orbits within planetary systems and proposing a fast optimization method for finding the burnup-optimal levered orbit is a pressing issue.

[0003] Research using Tr p Graph-based multi-satellite orbit optimization can improve orbit design efficiency and reduce the difficulty of solving complex multi-satellite orbit sequence problems. p The graph essentially depicts the possible states of a spacecraft over a period T and a pericenter radius r. p In the relationship diagram, without considering the inclination angle, the shape of the orbit can be determined by just two parameters. This is relevant when considering the relative glide satellite V. ∞ Under certain conditions, all states achievable through resonance-based leverage can be represented by a single curve in the diagram. Different V values ​​for the same leveraged satellite... ∞ The corresponding curves do not intersect. However, intersecting curves may exist between two different satellites using the same glide path. The state corresponding to the intersection point allows for the transfer between satellites without considering phase. Therefore, how to utilize Tr... p Accurately and quickly solving for the leverage sequence is a key problem that needs to be solved in multi-satellite leverage orbit optimization.

[0004] In the existing research on multi-satellite orbit optimization within planetary systems, the first technique [1] (Strange N J, Longuski J M. Graphical method for gravity-assist trajectory design[J]. Journal of Spacecraft and Rockets, 2002, 39(1):9-16.) proposed using Tr p The diagram helps assess potential celestial-assisted orbits and indicates that, from an energy perspective, such orbits are feasible. If in Tr pThe orbit that does not exist in the figure is definitely not feasible. The method can help reduce the feasible solution search range and improve the orbit design efficiency, but is only applicable to interplanetary gravity-assisted transfer orbit design and has high artificial participation.

[0005] In the prior art[2](Jingquan, Li Mingtao, Wang Youliang. State sequence search and optimization of transfer trajectory for Jupiter IV exploration mission [J]. Journal of Beijing University of Aeronautics and Astronautics, 1-14.), in view of the high artificial participation in designing the Jupiter IV orbit around the transfer orbit using the graphical analysis method, a state sequence search algorithm for the transfer orbit based on graphical analysis is proposed. Through analysis of the graphical properties and resonant gravity assistance, the change mode of the orbit state is summarized, and a search algorithm for iteratively solving the state transition sequence is proposed. The method solves the problem of high artificial participation when using the graphical analysis method, but in the implementation process, only the search is performed on the feature point set corresponding to a specific period in a specific curve in the figure, and only the feasible solution but not the optimal solution can be obtained. SUMMARY

[0006] To solve the problem that the existing method cannot solve the multi-satellite gravity-assisted fuel-optimal orbit in a planetary system, the purpose of the present application is to provide a multi-satellite gravity-assisted fuel-optimal orbit optimization method in a planetary system. p The method deduces the relationship between the gravity-assisted parameters and the transfer orbit state, establishes a multi-satellite gravity-assisted orbit parameter recursive model, takes the transfer orbit parameters between different satellites as variables, takes the minimum total velocity increment as the objective function to optimize the orbit, and obtains the multi-satellite gravity-assisted fuel-optimal orbit in a planetary system.

[0007] The purpose of the present application is realized through the following technical solutions:

[0008] The multi-satellite gravity-assisted fuel-optimal orbit optimization method in a planetary system disclosed by the present application selects a gravity-assisted satellite sequence that meets the exploration requirements, takes the period of each transfer orbit as a variable, adds a maneuver only when leaving the initial orbit and being captured by the target satellite or entering the target orbit, optimizes the objective function to be the minimum total velocity increment, and constructs an orbit optimization problem in the planetary center inertia system. When solving the optimization problem, the state of the next satellite gravity assistance can be solved according to the parameters of the spacecraft in the transfer orbit between the gravity-assisted satellites, the parameters of the transfer orbit after gravity assistance can be matched from the state, and continuous solving is realized until the state of reaching the target satellite is obtained. Finally, a genetic algorithm is used to solve the optimization problem, and the relative velocity of reaching each gravity-assisted satellite and the transfer orbit parameters before and after gravity assistance are obtained, that is, the multi-satellite gravity-assisted fuel-optimal orbit in a planetary system is obtained, and the transfer of each segment of the transfer orbit under the reference ephemeris is realized according to the transfer orbit parameters.

[0009] The multi-satellite gravity-assisted fuel-optimal orbit optimization method in a planetary system disclosed by the present application comprises the following steps:

[0010] Step one, the initial orbit is a large elliptical orbit around a planet, determine the target satellite and its surrounding circular orbit height, select the sequence of the flyby satellite to meet the detection requirements, take the transfer orbit period between the first flyby satellite and the adjacent flyby satellite in the sequence as the variable, only add the impulse when leaving the initial orbit and being captured by the target satellite, the optimization target is to minimize the total velocity increment, and construct the orbit optimization problem in the planet center inertia system;

[0011] The initial orbit period is T0, and the pericenter radius is r p，0 ; The sequence of flyby satellites is M1, M2, …, M n , and the sequence of flyby satellites is allowed to repeat, but the adjacent flyby satellites are different; In order to simplify the expression, the target satellite is also included in the flyby satellite sequence but does not participate in flyby, which is denoted as M n ; The transfer orbit periods as variables are T1, T2, …, T n , where T1 is the transfer orbit period between the initial orbit and the flyby satellite M1, and T2, …, T n are the transfer orbit periods between adjacent flyby satellites in the sequence;

[0012] Only two impulse maneuvers are added in the entire transfer orbit process, the first is to add a maneuver at the apocenter of the initial orbit to adjust the pericenter height, and then enter the transfer orbit before the flyby satellite M1, this maneuver is denoted as PEM, the second is the insertion maneuver when entering the target satellite's orbit, denoted as OIM; In summary, the total velocity increment of the transfer orbit is estimated as:

[0013] Δv Seq = Δv PEM + Δv OIM (1)

[0014] Where Δv PEM is the apocenter maneuver size, and Δv OIM is the insertion maneuver size.

[0015] There are two constraints in the optimization process, first, the period T i of each transfer orbit in the flyby sequence needs to be within the range [T min,i , T max,i ] that can be reached after the satellite flyby, and the other is that the pericenter radius r p,j of each transfer orbit needs to be higher than the set safety radius R lim ;

[0016] The orbit optimization problem is modeled as follows in the planet center inertia system:

[0017] x = (T1, T2,..., Tn) n )

[0018] minf(x) = Δv Seq

[0019] s.t. T min,i ≤ T i ≤ T max,i (i = 2, 3,..., n)

[0020] r p,j ≥ R lim (j = 1, 2,..., n)

[0021] Step two, add a maneuver at the apocenter of the initial orbit, the spacecraft enters the transfer orbit and finally intersects with the M1, according to the period of the transfer orbit, the semi-major axis of the transfer orbit is obtained, and then the apocenter velocity of the initial orbit and the transfer orbit is obtained respectively, and the size of the apocenter maneuver is solved by the velocity difference;

[0022] The relationship between the semi-major axis and the period is shown in formula (2),

[0023]

[0024] In the formula, μ is the planetary gravitational constant; T is the orbital period; a is the orbital semi-major axis.

[0025] According to formula (2), the initial orbit semi-major axis a0 and the transfer orbit semi-major axis a1 are obtained from T0 and T1 respectively.

[0026] The apocenter radius r p，0 is obtained from the initial orbit semi-major axis a0 and the pericenter radius r a,0 :

[0027] r a，0 = 2a0-r p，0 (3)

[0028] The transfer orbit apocenter is the same as the initial orbit apocenter, that is, r a，1 = r a,0 .

[0029] The apocenter velocity formula is shown in formula (4),

[0030]

[0031] In the formula, r a is the orbital pericenter radius, and V ra is the orbital apocenter velocity.

[0032] a0 and r a,0 and a1 and r a,1 are brought into formula (4) respectively, and the initial orbit apocenter velocity Vra,0 and the transfer orbit apocenter velocity V ra,1 ;

[0033] The velocity increment direction applied at the apocenter is the velocity tangent direction, so Δv PEM :

[0034] Δv PEM = |V ra,1 -V ra,0 | (5)

[0035] Step three, the velocity of the spacecraft relative to the M1 is calculated according to the transfer orbit parameters before reaching the M1 and the orbit parameters of the M1;

[0036] The eccentricity e1 is calculated from the semi-major axis a1 of the transfer orbit before reaching the M1 and the apocenter radius r a，1 :

[0037]

[0038] Based on the Tisserand criterion, formula (7) is derived to calculate the eccentricity of the transfer orbit from the velocity state relative to the M1;

[0039]

[0040] In formula (7), e is the eccentricity of the transfer orbit, a is the semi-major axis of the transfer orbit, R is the semi-major axis of the orbit of the M1, V is the velocity of the spacecraft, V ∞ is the velocity of the spacecraft relative to the M1, and α is the included angle between the V ∞ vector and the velocity V P vector of the M1.

[0041] Formula (8) is derived from formula (7):

[0042]

[0043] According to the relative relationship between the velocity of the spacecraft and the velocity of the M1, the following is obtained:

[0044]

[0045] In formula (9), V P is the velocity of the M1 in the revolution around the central celestial body.

[0046] Formula (10) is derived from formula (9):

[0047]

[0048] By combining formula (8) and formula (10) and eliminating α, the following is obtained:

[0049]

[0050] Calculate the velocity V of the spacecraft when it reaches the satellite using a1. 1,in :

[0051]

[0052] R1 is the semi-major axis of the orbit of satellite M1;

[0053] From equation (11), we obtain the expression for satellite M1 before leveraging the force:

[0054]

[0055] Among them, V ∞,1 V is the relative velocity of the spacecraft when it arrives at the lever satellite M1. P,1 To take advantage of the orbital speed of satellite M1 around the central celestial body;

[0056] In equation (13), except for V ∞,1 All quantities other than those mentioned above are known quantities. According to equation (13), the V that meets the conditions is obtained. ∞,1 ;

[0057] Step 4: Based on the spacecraft's arrival at satellite M i relative velocity V ∞,i (i = 1, 2, ..., n-1), the velocity state after leveraging the force is calculated for different leveraging angles using the leveraging model, and the mapping relationship between the orbital period and the radius of the pericenter after leveraging the force is obtained, i.e., Tr. p The image shows satellite M. i In V ∞,i The corresponding Tr p The curve; by interpolating the period of the transfer orbit after leveraging the glide path, the pericentric radius of the transfer orbit is obtained, and then the semi-major axis and eccentricity of the transfer orbit are calculated; based on the parameters of the transfer orbit, the arrival time of the spacecraft at the next lever-guided satellite M can be determined. i+1 relative velocity V ∞,i+1 .

[0058] According to equation (9), the spacecraft passes through satellite M. i Speed ​​V after leveraging i,out :

[0059]

[0060] Because V in equation (14) P,i and V ∞,i All are known quantities, α out For the spacecraft to leave the lever satellite M i The direction of the velocity vector and the lever satellite Mi The included angle between the revolution velocity vector direction, thus V i,out is expressed as a function only related to α out :

[0061] V 1,out = g(α out ) (15)

[0062] α out is any value in the theoretical range [0, 2π]; thus different α out corresponds to the satellite M i after the transfer orbit period T i+1 and the perihelion radius r p,i+1 , and then T i+1 -r p,i+1 relationship curve; according to the relationship curve, the range of T i+1 value is determined [T min,i+1 , T max,i+1 ], and the variable T i+1 is calculated by interpolation to obtain r p,i+1 .

[0063] The semi-major axis a i+1 of the transfer orbit is obtained from T i+1 :

[0064] The eccentricity e p,i+1 of the transfer orbit is obtained from r i+1 and a i+1 :

[0065]

[0066]

[0067] All quantities in formula (19) except V ∞,i+1 are known quantities, and V ∞,i+1 that meets the conditions is obtained according to formula (19).

[0068] Step five, repeat step four until the relative velocity V ∞,n to the target satellite is solved;

[0069] Step six, according to the V ∞,n to the target satellite and the circular orbit height r h , the size of the insertion maneuver into the circular orbit is solved;

[0070] The radius R circle of the circular orbit is expressed as:

[0071] R circle = radn +r h (20)

[0072] where rad n denotes the physical radius of the target satellite;

[0073] The linear velocity v circle of the circular orbit is solved according to formula (21):

[0074]

[0075] where μ n is the gravitational constant of the target satellite;

[0076] The velocity v in of the spacecraft when reaching the height of the circular orbit is:

[0077]

[0078] The insertion maneuver Δv OIM of the spacecraft into the circular orbit is solved according to formula (23):

[0079] Δv OIM = |v circle -v in | (23)

[0080] Step seven, a global optimization algorithm is used to solve the orbit optimization problem, and the relative velocity V ∞,i reaching each borrowing satellite and the corresponding transfer orbit period T i are obtained; without considering ephemeris, the V ∞,i and T i complete each segment of the transfer orbit optimization.

[0081] Step eight, if the orbit conversion before and after borrowing cannot be realized by one-time borrowing, the spacecraft can realize the orbit conversion before and after borrowing by multiple resonance borrowings of the same borrowing satellite.

[0082] Beneficial effects:

[0083] 1. The multi-satellite borrowing fuel consumption optimal orbit optimization method in the planetary system disclosed by the application is based on T-r p figure to deduce the relationship formula between the borrowing state and the transfer orbit, establish a multi-satellite borrowing orbit state recursion model, use the transfer orbit period between adjacent satellites in the borrowing sequence as a variable, realize the fast solution of the multi-satellite borrowing orbit of the known sequence, and improve the optimization efficiency of the multi-satellite borrowing fuel consumption optimal orbit in the planetary system.

[0084] 2. The present invention discloses a method for optimizing the optimal orbit for multiple satellites leveraging transfer orbits within a planetary system by establishing an optimization model for multiple satellites leveraging transfer orbits within a planetary system and solving for the optimal transfer orbit with the transfer orbit period of each segment as a variable. This method is applicable to the optimization of multiple satellites leveraging transfer orbits within any planetary system and has universality.

[0085] 3. The present invention discloses a method for optimizing the orbit of multiple satellites in a planetary system by leveraging fuel consumption. It selects a sequence of satellites that meet the exploration requirements, uses the period of each transfer orbit as a variable, and adds maneuvers only when leaving the initial orbit and when being captured by the target satellite or entering the target orbit. The optimization objective is to minimize the total velocity increment. The orbit optimization problem is constructed in the planetary center-of-mass inertial frame, which improves the convergence of the optimization results and makes the optimization results close to the theoretical optimal solution. Attached Figure Description

[0086] Figure 1 This is a flowchart of a method for optimizing the orbit of multiple satellites in a planetary system by leveraging fuel consumption, according to the present invention.

[0087] Figure 2 This is the Tr corresponding to the optimization result in step seven of the specific implementation method of the present invention. p picture;

[0088] Figure 3 This is a schematic diagram of the leverage track corresponding to the optimization result in step seven of the specific implementation of the present invention. Detailed Implementation

[0089] To better illustrate the purpose and advantages of the present invention, the invention will be further described below in conjunction with the accompanying drawings and examples.

[0090] like Figure 1 As shown in this embodiment, a method for optimizing the orbit of multiple satellites in a planetary system by leveraging fuel consumption is disclosed. The specific implementation steps are as follows:

[0091] Taking the design of a transfer orbit within the Jupiter system from the initial capture orbit of Jupiter to the orbit around Ganymede as an example, the initial orbit is a highly elliptical orbit with a period T0 of 200 days and a pericentric radius r. p，0 With a radius 10 times that of Jupiter, Callisto orbits at an altitude of 500 km. Based on experience, the sequence of satellites using the booster orbit is set as Ganymede-Callisto-Gallisto-Callisto (GCGC). The optimal burnup design method proposed in this patent is used to optimize the booster orbit for multiple satellites.

[0092] Step 1: Construct an orbit optimization problem in the Jupiter center-of-mass inertial frame. The initial orbit is a highly elliptical orbit around Jupiter. The transfer orbit period before entering the first borrowed satellite (G) and the transfer orbit period between adjacent borrowed satellites in the sequence are used as variables. Maneuvers are added only when leaving the initial orbit and when captured by Callisto. The optimization objective is to minimize the total velocity increment.

[0093] The initial orbital state S0 has a period T0 = 200 days and a pericentric radius r. p，0 =10 rad Jupiter , of which rad Jupiter denoted as Jupiter's physical radius. The sequence of the lever satellites is GCGC. The periods of each transfer orbit, which are variables, are T1, T2, T3, and T4, where T1 is the period of the transfer orbit from the initial orbit to the first lever satellite (G), and T2, T3, and T4 are the periods of the transfer orbits between adjacent lever satellites in the sequence.

[0094] Only two pulse maneuvers were added during the entire transfer orbit process. The first was a maneuver added at the apocentric point of the initial orbit to adjust the pericentric altitude, before entering the transfer orbit prior to the first steerable satellite (G). This maneuver is denoted as PEM. The second was an insertion maneuver, denoted as OIM, when entering the target satellite's orbit. In summary, the total velocity increment of the transfer orbit is estimated as follows:

[0095] Δv Seq =Δv PEM +Δv CIM (twenty four)

[0096] The optimization process requires two constraints. The first is the period T of each transfer trajectory segment in the leverage sequence. i All need to be within the achievable range of the cycle after satellite leverage [T] min,i ,T max,i Inside, another reason is for safety; the radius r of the pericenter of each transfer orbit segment... p,j All must be higher than the set safety radius R lim =8rad Jupiter .

[0097] The optimization problem is modeled as follows:

[0098] x=(T1,T2,…,T n )

[0099] minf(x)=Δv Seq

[0100] st T min,i ≤T i ≤T max,i (i = 2, 3, ..., n)

[0101] r p,j ≥R lim (j=1,2,…,n)

[0102] Step two, add a maneuver at the apocenter of the initial orbit to enter the transfer orbit, which will eventually intersect with the first gravity assist satellite (G). According to the period of the transfer orbit, the semi-major axis of the transfer orbit can be obtained, and then the apocenter velocities of the initial orbit and the transfer orbit can be obtained respectively. The size of the apocenter maneuver can be solved by the velocity difference.

[0103] According to the relationship between the semi-major axis and the period, as shown in equation (25):

[0104]

[0105] In equation (25), μ J is the gravitational constant of Jupiter.

[0106] According to equation (25), the semi-major axis a0 of the initial orbit and the semi-major axis a1 of the transfer orbit can be obtained respectively from T0 and T1.

[0107] From the semi-major axis a0 of the initial orbit and the pericenter radius r p，0 , the apocenter radius r a,0 of the initial orbit can be obtained.

[0108] r a，0 =2a0-r p，0 (26)

[0109] The apocenter of the transfer orbit is the same as that of the initial orbit, i.e. r a，1 =r a,0 . According to the apocenter velocity formula as shown in equation (27):

[0110]

[0111] From the semi-major axis and the apocenter radius of the initial orbit and the transfer orbit, the apocenter velocity V ra,0 of the initial orbit and the apocenter velocity V ra,1 of the transfer orbit can be obtained.

[0112] When adjusting the pericenter radius of the transfer orbit, the velocity increment direction is the tangential direction of the velocity, so Δv PEM can be obtained:

[0113] Δv PEM =|V ra,1 -V ra,0 | (28)

[0114] Step three, according to the parameters of the transfer orbit before reaching the first gravity assist satellite (G) and the parameters of the gravity assist satellite orbit, the velocity of the spacecraft relative to the first gravity assist satellite (G) can be obtained.

[0115] The semi-major axis a1 of the transfer orbit before reaching the first satellite (G) and the apocenter radius r a，1 The eccentricity e1 can be obtained:

[0116]

[0117] Based on the Tisserand criterion, formula (30) can be derived, which can be solved by the velocity state of the spacecraft relative to the satellite to obtain the eccentricity of the transfer orbit.

[0118]

[0119] In formula (30), e is the eccentricity of the transfer orbit, a is the semi-major axis of the transfer orbit, R is the semi-major axis of the orbit of the satellite, V is the velocity of the spacecraft, V ∞ is the velocity of the spacecraft relative to the satellite, and a is the included angle between the V ∞ vector and the satellite velocity V P vector.

[0120] Formula (31) can be derived from formula (30):

[0121]

[0122] According to the relative relationship between the spacecraft velocity and the satellite velocity, we have:

[0123]

[0124] Formula (33) can be derived from formula (32):

[0125]

[0126] By combining formula (31) and formula (33), we can eliminate a to obtain:

[0127]

[0128] The velocity V 1,in of the spacecraft when it reaches the satellite can be obtained from the semi-major axis a1 of the transfer orbit before the first satellite (G):

[0129]

[0130] where R1 is the semi-major axis of the orbit of the first satellite (G).

[0131] From formula (34), we can obtain the expression before the first satellite:

[0132]

[0133] where V ∞,1 is the relative velocity of the spacecraft to the first gravitational assist satellite (G) at the time of arrival, V P,1 is the orbital velocity of the first gravitational assist satellite (G) around the central satellite.

[0134] In equation (36), all quantities except V ∞,1 are known, so V ∞,1 that satisfies the condition can be solved.

[0135] Step four, based on the relative velocity V i of the spacecraft to the gravitational assist satellite M ∞,i (i = 1, 2,..., n-1), the speed state after gravitational assistance corresponding to different gravitational assistance angles is calculated through the gravitational assistance model, and the mapping relationship between the orbital period and the pericenter radius after gravitational assistance, i.e. T-r p is obtained. The T-r i curve corresponding to V ∞,i of the gravitational assist satellite M p is shown in the figure; the pericenter radius of the transfer orbit is obtained by interpolation calculation using the period of the transfer orbit as a variable, and then the semi-major axis and eccentricity of the transfer orbit are obtained; the relative velocity V i+1 of the spacecraft to the next gravitational assist satellite M ∞,i+1 is obtained according to the parameters of the transfer orbit.

[0136] According to equation (32), the speed V i of the spacecraft after gravitational assistance of the gravitational assist satellite M i,out is obtained:

[0137]

[0138] Since V P,i and V ∞,i in equation (37) are known quantities, and α out is the angle between the direction of the velocity vector of the spacecraft when leaving the gravitational assist satellite M i and the direction of the orbital velocity vector of the gravitational assist satellite M i , V i,out is expressed as a function only related to α out :

[0139] V 1,out = g(α out ) (38)

[0140] α out is arbitrarily selected within the theoretical range [0, 2π]; therefore, the period T i+1 and the pericenter radius r of the transfer orbit after gravitational assistance of the gravitational assist satellite M i corresponding to different α out can be solved.p,i+1 , and then T i+1 is obtained p,i+1 , and then the relationship curve is determined i+1 , and then the range of values of T min,i+1 , T max,i+1 ] is determined i+1 , and then the variable T p,i+1 is interpolated to obtain r .

[0141] The semi-major axis a i+1 of the transfer orbit is obtained from T i+1 :

[0142]

[0143] The eccentricity e p,i+1 of the transfer orbit is obtained from r i+1 and a i+1 :

[0144]

[0145] The velocity V i+1 of the spacecraft when it reaches the M i+1 satellite is obtained from a i+1,in :

[0146]

[0147] In the formula, R i+1 is the semi-major axis of the orbit of the M i+1 satellite.

[0148] According to formula (36), formula (42) is obtained:

[0149]

[0150] In the formula, V P,i+1 is the orbital velocity of the M i+1 satellite around the central celestial body.

[0151] All quantities in formula (42) except V ∞,i+1 are known quantities, so V ∞,i+1 that meets the conditions is obtained.

[0152] Step five, repeat step four until the relative velocity V ∞,n to the target satellite (C) is solved.

[0153] Step six, according to V ∞,n to the target satellite and the circular orbit height r h , the size of the insertion maneuver into the circular orbit is solved.

[0154] R circle may be expressed as:

[0155] R circle = rad C + r h = 2400 km + 500 km = 2900 km (43)

[0156] where rad C represents the physical radius of Jupiter IV.

[0157] The linear velocity v circle of the circular orbit can thus be solved:

[0158]

[0159] where μ C is the gravitational constant of the target satellite.

[0160] The velocity v in of the spacecraft when it reaches the height of the circular orbit is:

[0161]

[0162] Therefore, the insertion maneuver Δv OIM of the spacecraft into the circular orbit is solved:

[0163] Δv OIM = |v circle -v in | (46)

[0164] Step seven, the optimization problem constructed by the first six steps is solved by using a global optimization algorithm to obtain the relative velocity V ∞,i to each borrowed satellite and the corresponding transfer orbit period T i ; without considering ephemeris, the transfer orbit design of each segment is completed through V ∞,i and T i .

[0165] The initial parameters of the borrowed sequence obtained by solving the optimization problem are shown in the following table, and the sizes of the two maneuvers are Δv PEM = 0.11 km / s and Δv OIM = 1.03 km / s, so the total speed increment is 1.14 km / s.

[0166] Table 1 Summary of preliminary design parameters of borrowed sequence

[0167]

[0168]

[0169] Step eight, if the transfer of orbit before and after the power cannot be achieved by one time, the spacecraft through the same satellite resonance multiple times to achieve the transfer of orbit before and after the power.

[0170] The T-r of the orbit of the power p Figure as Figure 2 The initial sequence is designed. The transfer of orbit period from 204.01 days to 9.27 days cannot be achieved by one time, so the transfer is achieved by five times of resonance power of Jupiter 3. Similarly, the transfer of orbit period from 9.27 days to 13.02 days cannot be achieved by one time, so the transfer is achieved by two times of resonance power of Jupiter 4.

[0171] The power sequence of the task and the detailed orbit parameters are shown in the following table:

[0172] Table 2 Summary of power sequence and orbit parameters

[0173]

[0174] The final transfer orbit is shown in the figure. Figure 3

[0175] The above specific description further describes the purpose, technical scheme and beneficial effects of the application. It should be understood that the above description is only a specific embodiment of the application and is not used to limit the protection scope of the application. Any modification, equivalent replacement, improvement, etc. within the spirit and principles of the application shall be included in the protection scope of the application.​

Claims

1. A method for optimizing multi-satellite transfer orbit in a planet system, characterized in that: Comprising the following steps, Step one, the initial orbit is a large elliptical orbit around the planet, the target satellite and the satellite's orbit height are determined, the sequence of the borrowed satellite that meets the detection requirements is selected, the transfer orbit period before the first borrowed satellite and the transfer orbit period between adjacent borrowed satellites in the sequence are taken as variables, only the maneuver when leaving the initial orbit and being captured by the target satellite is added, and the optimization target is the minimum total velocity increment, and the orbit optimization problem is constructed in the planet center inertia system; Step two, a maneuver is added at the apsis of the initial orbit, the spacecraft enters the transfer orbit and finally intersects with the borrowed satellite M1, the semi-major axis of the transfer orbit is calculated according to the period of the transfer orbit, and then the apsis velocities of the initial orbit and the transfer orbit are calculated respectively, and the size of the apsis maneuver is calculated according to the velocity difference; Step three, the velocity of the spacecraft relative to the borrowed satellite M1 is calculated according to the transfer orbit parameters before reaching the borrowed satellite M1 and the borrowed satellite orbit parameters; Step 4: Based on the spacecraft's arrival at satellite M i relative velocity V ∞,i (i = 1, 2, ..., n-1), the velocity state after leveraging the force is calculated for different leveraging angles using the leveraging model, and the mapping relationship between the orbital period and the radius of the pericenter after leveraging the force is obtained, i.e., Tr. p The image shows satellite M. i At relative velocity V ∞,i The corresponding Tr p The curve; using the period of the transfer orbit after leveraging the glide path as a variable, interpolation calculations are performed to obtain the pericentric radius of the transfer orbit, and then the semi-major axis and eccentricity of the transfer orbit are calculated; based on the parameters of the transfer orbit, the arrival time of the spacecraft at the next lever-path satellite M is determined. i+1 relative velocity V ∞,i+1 ; Step five, repeat step four until the relative velocity V to the target satellite is solved ∞,n ; Step six, solve the size of the insertion maneuver into the circular orbit according to the V ∞,n and the circular orbit height r h , solve the size of the insertion maneuver into the circular orbit Step seven, the global optimization algorithm is used to solve the orbit optimization problem, and the relative velocity V ∞,i and the corresponding transfer orbit period T i are obtained; without considering ephemeris, each segment transfer orbit is optimized through V ∞,i and T i .

2. The method of claim 1, wherein: Further comprising step eight, if the transfer of the orbit before and after the borrowing cannot be realized by one borrowing, the spacecraft realizes the transfer of the orbit before and after the borrowing by multiple resonance borrowings on the same borrowed satellite.

3. The method of claim 1 or 2, wherein: The implementation method of step one is, The initial orbit period is T0, and the pericenter radius is r p，0 ; the sequence of the transfer satellites is M1, M2, …, M n , and the sequence of the transfer satellites can have repeated transfer satellites, but adjacent transfer satellites are different; the target satellite is also included in the sequence of the transfer satellites but does not participate in the transfer, and is denoted as M n ; The periods of the transfer orbits as variables are T1, T2,..., T n where T1is the period of the transfer orbit between the initial orbit and the first M1, T2,..., T n are the periods of the transfer orbits between adjacent MIs in the sequence, respectively. Only two pulse maneuvers are added in the whole transfer orbit process, the first one is added at the apsis of the initial orbit to adjust the pericenter height, and then the spacecraft enters the transfer orbit before the borrowed satellite M1, this maneuver is recorded as PEM, and the second one is the insertion maneuver when entering the target satellite orbit, recorded as OIM; in summary, the total velocity increment of the transfer orbit is estimated as: Δv Seq = Δv PEM + Δv OIM (1) where Δv PEM is the telecentric point maneuver size, Δv OIM is the insertion maneuver size; Two constraints are needed in the optimization process. The first one is that the period T of each transfer orbit in the swing-by sequence i should be within the reachable range [T min,i ,T max,i ] after the satellite performs the swing-by maneuver. The other one is that the perigee radius r p,j of each transfer orbit should be higher than the safety radius R lim for safety consideration. The orbit optimization problem is modeled as follows in the planet center inertia system: x = (T1, T2,..., T n ) minf(x) = Δv Seq s.t.T min,i ≤T i ≤T max,i (i = 2, 3,..., n) r p,j ≥R lim (j = 1, 2,..., n).

4. The method of claim 3, wherein: The implementation method of step two is, The relationship between the semi-major axis and the period is shown in equation (2), In the formula, μ is the planetary gravitational constant; T is the orbit period; a is the orbit semi-major axis; According to equation (2), the initial orbit semi-major axis a0 and the transfer orbit semi-major axis a1 are calculated from T0 and T1 respectively; from the initial orbital semi-major axis a0and the pericenter radius r p，0 The apocenter radius r a,0 is found r a，0 ＝2a0-r p，0 (3) The transfer orbit apsis is the same as the initial orbit apsis, i.e. r a，1 = r a,0 ; The apsis velocity formula is shown in equation (4), where r a is the radius of the orbit's pericenter, V ra is the velocity of the orbit's apocenter; Substitute a0 and r a,0 and a1 and r a,1 into equation (4) respectively, and obtain the initial orbit parabolic point velocity V ra,0 and the transfer orbit parabolic point velocity V ra,1 ; The direction of the velocity increment applied at the apocenter is tangential to the velocity, so that Δv PEM : Δv PEM = |V ra,1 -V ra,0 | (5).

5. The method of claim 4, wherein: The implementation method of step three is, a1and the apocenter radius r a，1 The eccentricity e1is found: Based on Tisserand criterion, equation (7) is derived, and the eccentricity of the transfer orbit is calculated from the relative velocity state of the borrowed satellite; In formula (7), e is the eccentricity of the transfer orbit, a is the semi-major axis of the transfer orbit, R is the semi-major axis of the orbit of the piggyback satellite, V is the velocity of the spacecraft, V ∞ is the velocity of the spacecraft relative to the piggyback satellite, and a is the included angle between the direction of the V ∞ vector and the velocity V P of the piggyback satellite. Equation (8) is derived from equation (7): According to the relative relationship between the spacecraft velocity and the borrowed satellite velocity: In the formula, V P is the speed of the piggyback satellite around the central celestial body. Equation (10) is derived from equation (9): Equations (8) and (10) are combined to eliminate α to obtain: V = V0+ a1t 1,in : Wherein, R1 is the semi-major axis of the borrowed satellite M1 orbit; Equation (11) is obtained about the satellite M1 before borrowing: wherein V ∞,1 is the relative velocity of the spacecraft when it reaches the M1 satellite, V P,1 is the orbital velocity of the M1 satellite around the central celestial body; In formula (13), all quantities except V ∞,1 are known quantities, and V ∞,1 is obtained from formula (13) subject to the condition that 6. The method of claim 5, wherein: The implementation method of step four is, According to equation (9) the spacecraft passes the M i velocity V after the flyby i,out : Since V P,i and V ∞,i are known quantities, α out is the angle between the direction of the velocity vector of the spacecraft as it departs from the M i and the direction of the orbital velocity vector of the M i , and thus V i,out is expressed as a function of α out alone: V 1,out = g(a out ) (15) alpha out Any value in the theoretical range [0, 2pi]; solve different alpha out Corresponding to the boost satellite M i Period T of the post-transfer orbit i+1 And the perihelion radius r p,i+1 , and then get T i+1 -r p,i+1 Relationship curve; according to the relationship curve to determine T i+1 The range of values [T min,i+1 , T max,i+1 ] will be variable T i+1 Interpolation calculation to get r p,i+1 ; By T i+1 The semi-major axis a of the transfer orbit is obtained i+1 : r p,i+1 and a i+1 The transfer orbit eccentricity e i+1 is obtained: By a i+1 The velocity V of the spacecraft when it reaches the M i+1 maneuvering satellite is found by i+1,in : wherein R i+1 MEO satellite i+1 semi-major axis of the orbit Equation (19) is obtained from equation (13): In the formula, V P,i+1 MEO satellite i+1 The orbital speed around the central celestial body; All quantities in formula (19) except V ∞,i+1 are known quantities, and V ∞,i+1 is obtained from formula (19) subject to the condition 7. The method of claim 6, wherein: The implementation method of step six is, The radius of the circular orbit R circle is represented as: R circle = rad n + r h (20) wherein rad n represents the physical radius of the target satellite; Solving for the linear velocity v around the circular orbit according to equation (21) circle : where μ n is the gravitational constant of the target satellite; The velocity v of the spacecraft when it reaches the height of the circular orbit in is: The spacecraft insertion maneuver delta-v is found from equation (23) as the spacecraft enters a circular orbit OIM : Δv OIM = |v circle -v in | (23).

Images