Orbit Design Method for Large-Scale Inclination Change Using Multiple Gravity Assists

Abstract

The present invention discloses an orbital design method for large-range inclination change using multiple gravity assists. The orbital design method for large-range inclination change using multiple gravity assists has the following input information: assuming that the orbital inclination of the gravity assist small body around the planet is 0, at a certain moment, the spacecraft has achieved coincidence of its position with that of the gravity assist body through orbital control, the positions and velocities of the spacecraft and the small body relative to the planet are known, the orbital inclination of the spacecraft at this moment is i0, and the target inclination to be achieved is i f , then the inclination adjustment is achieved in the shortest time through multiple gravity assists. The present invention solves the problems in the prior art that there are too many variables to solve during the orbital inclination change process and large-range inclination change cannot be achieved.

Description

Technical Field

[0001] The present invention belongs to the technical field of aerospace navigation control, and particularly relates to an orbital design method for changing the inclination angle in a large range by using multiple gravity assists. Background Art

[0002] In the gravity assist approximation model, only the direction of the relative velocity can be changed. Therefore, the inclination angle change is coupled with the semi-major axis change. When the change amount of the gravity assist velocity direction in a single time is very small, it is necessary to accumulate it multiple times. During this process, the orbital period changes greatly, and the resonance ratio needs to be continuously changed until the orbital inclination control task is completed. How to plan the direction of each gravity assist and select the corresponding resonance ratio to enable the aircraft to complete the large-range inclination adjustment task in the shortest time is a complex optimization problem. The present invention proposes a mixed-integer programming model for changing the inclination angle by using multiple gravity assists, which can overall plan the gravity assist design variables and complete the large-range orbital transfer fastest under the premise of meeting various constraints. The method can quickly obtain the optimal gravity assist sequence and realize the large-range change of the detector's orbital inclination without using fuel. Summary of the Invention

[0003] The purpose of the present invention is to provide an orbital design method for changing the inclination angle in a large range by using multiple gravity assists, which solves the problems of too many variables to be solved and unable to realize the large-range inclination angle change in the prior art during the process of changing the orbital inclination angle.

[0004] The technical solution adopted by the present invention is an orbital design method for changing the inclination angle in a large range by using multiple gravity assists. The input information is as follows: assuming that the orbital inclination angle of the gravity assist small body orbiting the planet is 0, at a certain moment, the spacecraft has achieved the coincidence of the position with the gravity assist celestial body through orbital control, the position and velocity of the spacecraft and the small body relative to the planet are known, the orbital inclination angle of the spacecraft at this moment is i0, and the target inclination angle to be achieved is i f , then the inclination angle adjustment is realized in the shortest time by using multiple gravity assists, and the specific implementation steps are as follows:

[0005] Step 1: Calculate the angle of change of the relative velocity of the gravity assist according to the initial orbit, the target orbit elements, and the position and velocity.

[0006] Step 2: Calculate all available resonance period ratios during the transfer process according to the initial orbital period and the target orbital period.

[0007] Step 3: Establish a mixed-integer optimization model for multiple gravity assists.

[0008] Step 4: Use the differential evolution algorithm to optimize and obtain the result.

[0009] The characteristics of the present invention also lie in

[0010] Step 1 is specifically implemented according to the following steps:

[0011] According to the known gravity assist approximation model in the field, the gravity assist is an instantaneous process, and there is:

[0012]

[0013] where t is the gravity assist moment, t - , t + respectively represent the moments before and after the gravity assist instantaneously, r(t - ), r(t + ) respectively represent the position vectors of the spacecraft before and after the gravity assist instantaneously, r G (t) represents the position vector of the gravity assist celestial body, v(t - ), v(t + ) are respectively the velocity vectors of the spacecraft before and after the gravity assist instantaneously, v G (t) is the velocity vector of the gravity assist celestial body, respectively represent the relative velocity vectors of the spacecraft and the celestial body before and after the gravity assist instantaneously, v ∞ is the scalar of the relative velocity magnitude, δ is the angle between and caused by the gravity assist, and is calculated according to the following formula:

[0014]

[0015] where μ G is the gravitational constant of the gravity assist celestial body, h p is the altitude at which the spacecraft flies over the celestial body, R G is the radius of the celestial body, h p is a variable, usually h p has the following constraint:

[0016] h p <h p min (3)

[0017] h p min is the minimum value of the altitude at which the spacecraft flies over the celestial body. From equation (2), it can be known that when other inputs remain unchanged, h p min constrains the maximum value δ max of δ, that is:

[0018]

[0019] Given the initial orbit and the target orbit of the spacecraft, the corresponding spacecraft velocity can be obtained according to the known methods in the field. Assume that the position coordinates of the spacecraft and the gravity assist celestial body in the inertial system are rG =(0, -a G , 0), a G represents the semi-major axis of the orbit of the gravity assist celestial body. Since the gravity assist celestial body is approximately in a circular orbit with an inclination of 0, thus the modulus of r G is always equal to a G , when r G takes other values, through the coordinate rotation between inertial systems, make r G equal to (0, -a G , 0) in the new coordinate system,

[0020] Let the velocity of the gravity assist celestial body be where μ J is the planetary gravitational constant. Let the inclination of the initial orbit be i0 and the velocity vector be Then the relative velocity can be obtained from Equation (5) The velocity of the target orbit is written as Then the initial relative velocity of the spacecraft can be obtained and the included angle Δ between the relative velocity of the target orbit :

[0021]

[0022] where v G (t) is the velocity of the gravity assist celestial body, Δ is much larger than the maximum angle δ changed by a single gravity assist max , so it is necessary to realize step by step the rotation of to through multiple gravity assists. The total angle turned is Δ, and v ∞ is the modulus of

[0023] In step 2,

[0024] Given that the initial orbit period is T0, the target orbit period is T f , and the orbit period of the gravity assist celestial body is T G , then obviously the resonance ratio after each gravity assist is between and . Therefore, it is necessary to find all possible resonance ratios M, N are natural numbers, N ≤ N max , N max is the given threshold, indicating that the time interval between two gravity assists is less than N max T G , for selection in the next step of optimization

[0025] Step 2 is specifically implemented according to the following steps:

[0026] The specific process of obtaining all available resonance period ratios is as follows: Start traversing backward from N = 1. For a given N and all M within the range of 1 ≤ M ≤ N, retain the M that satisfies the constraint If T f <T0, then the constraint is taken as The i-th value that satisfies the constraint is denoted as M N,i , and the number of M N,i is denoted as n N , obtaining the corresponding sequence as Then increment N by 1 and continue to search for the M and sequence that satisfy the constraint until N > N max when it ends,

[0027] Finally, merge all the sequences and sort them by size, denoted as

[0028] Step 3 is specifically implemented according to the following steps:

[0029] According to the following formula:

[0030]

[0031] where v = v ∞ + v G is the spacecraft velocity vector corresponding to any given v ∞ , α is the semi-major axis of the spacecraft orbit. When α is constant, v ∞ ·v G is constant, then is also constant, and the corresponding semi-major axis of the orbit remains unchanged. And the orbital period is uniquely determined by the semi-major axis, indicating that for a certain resonance period ratio the corresponding α is obtained by back-calculating through the spacecraft orbital period i ,

[0032]

[0033] where, T i is the spacecraft period corresponding to the resonance period ratio , a and v are intermediate variables in the operation; then the v ∞ before and after each gravity assist is on the conical bottom surface corresponding to a certain α i ; where, α i represents the angle between the before the i-th gravity assist and v G , β i represents the angle between the projection of on the conical bottom surface corresponding to α i and the z-axis;

[0034] In summary, the Denoted as n, the optimization variables include {x i} and {β i}, both {x i} and {β i} are 2n-dimensional. Among them, the first n x i are 0-1 variables, indicating whether the optimal transfer trajectory passes through the corresponding periodic orbit, where 0 means not passing through and 1 means passing through; the last n β i are real variables, used to locate the when the spacecraft leaves the current orbit period and transfers to the next orbit period. Let β0 be the initial f angle between the projection of the cone bottom corresponding to α0 on the z-axis, and β be the target f angle between the projection of the cone bottom corresponding to α i on the z-axis. Then, for a given set of {x i} and {β

[0035] }, the total transfer time is calculated according to the following process: i Traverse backward from x1 to find the i for which the x value is 1, and the corresponding i α i and β i α i and β ∞ can directly calculate

[0036]

[0037] That is, the first gravity assist represented by x i and β i changes the relative velocity of the spacecraft from to The obtained angle value is:

[0038]

[0039] Then, according to the fact that the maximum angle of changing the relative velocity by a single gravity assist is a fixed value δ max , check the constraint. Obviously, when |α i - α0| > δ max , the gravity assist cannot be achieved, and a penalty term needs to be added to the objective function; when δ max ≤ |α i - α0| ≤ Δ i , a single gravity assist can be achieved, and the spacecraft intersects with the position of the gravity assist celestial body again after N i T i time; and when Δi > δ max When it is, the orbit transfer process is split into multiple gravity assists without changing the period. That is, first use one gravity assist to change the resonance period ratio of the spacecraft orbit relative to the gravity assist celestial body orbit, and then keep the resonance period ratio unchanged and use multiple gravity assists to make the relative velocity vector of the spacecraft reach

[0040] Calculate the included angle δ with by solving the following system of equations max , and the included angle with the x-axis is α i of the vector, denoted as

[0041]

[0042]

[0043] Then find corresponding β i , then continue to transfer to The number of gravity assists is The integer obtained by rounding up, denoted as m i , then from Change to The required time is (m i +1)N i T i ;

[0044] Then continue to traverse backward to find the next x j value is 1, and the corresponding α j and β j ; Continue to calculate the time required to change from to in the same way, as well as the constraint penalty term;

[0045] Until the traversal ends at n. Assume that there are k non-zero variables in x i , i = 1,... n, corresponding to k + 1 times of changes. The total duration and constraint penalty term calculated are used as the optimization objective function:

[0046]

[0047] Among them, Δt j is the time required for each change, k p is the number of times the constraint is not satisfied, c p is the constant penalty factor;

[0048] In summary, {x is obtained.i} and {β i} are optimization variables, and a mixed-integer optimization model with J as the objective function.

[0049] The beneficial effect of the present invention is that an orbital design method for large-range inclination change using multiple gravity assists can perform a large-range adjustment of the orbital inclination through multiple gravity assists of small celestial bodies without consuming fuel, and through optimization, a multiple gravity assist trajectory with the shortest orbital transfer time can be obtained. The method has high solution efficiency and can provide a useful reference for the design and analysis of such tasks. Brief Description of the Drawings

[0050] Figure 1 is a schematic diagram of the relative velocity between the spacecraft and the gravity assist celestial body;

[0051] Figure 2(a) shows the change in the resonance ratio in the single-step gravity assist type;

[0052] Figure 2(b) shows that the resonance ratio remains unchanged in the single-step gravity assist type;

[0053] Figure 3 is the process of calculating the objective function according to the optimization variables;

[0054] Figure 4 is the flow chart of the present invention;

[0055] Figure 5(a) shows the X-Y plane coordinate system in the optimal trajectory of the spacecraft obtained in the embodiment; Figure 5(b) shows the three-dimensional coordinate system in the optimal trajectory of the spacecraft obtained in the embodiment;

[0056] Figure 6 is the path of the change in the velocity vector of the spacecraft relative to Ganymede obtained in the embodiment. Detailed Embodiment

[0057] The present invention will be described in detail below in conjunction with the drawings and specific embodiments.

[0058] The orbital design method of the present invention for large-range inclination change using multiple gravity assists has the following input information: assuming that the orbital inclination of the gravity assist small celestial body around the planet is 0, at a certain moment, the spacecraft has achieved the coincidence of the position with the gravity assist celestial body through orbital control, the position and velocity of the spacecraft and the small celestial body relative to the planet are known, the orbital inclination of the spacecraft at this moment is i0, and the target inclination to be achieved is i f , then the inclination adjustment is achieved in the shortest time through multiple gravity assists, and the flow chart is as Figure 4 shown, and the specific implementation is carried out according to the following steps:

[0059] Step 1: Calculate the change angle of the relative velocity by the gravity assist according to the initial orbit, target orbit elements, and position and velocity;

[0060] Step 1 is specifically implemented according to the following steps:

[0061] According to the known gravity assist approximation model in the field, the gravity assist is an instantaneous process, and there is:

[0062]

[0063] where t is the gravity assist moment, t - , t + respectively represent the moments before and after the gravity assist instantaneously, r(t - ), r(t + ) respectively represent the position vectors of the spacecraft before and after the gravity assist instantaneously, r G (t) represents the position vector of the gravity assist celestial body, v(t - ), v(t + ) are respectively the velocity vectors of the spacecraft before and after the gravity assist instantaneously, v G (t) is the velocity vector of the gravity assist celestial body, respectively represent the relative velocity vectors of the spacecraft and the celestial body before and after the gravity assist instantaneously, v ∞ is the scalar of the relative velocity magnitude, δ is the angle between and caused by the gravity assist, and is calculated according to the following formula:

[0064]

[0065] where μ G is the gravitational constant of the gravity assist celestial body, h p is the height at which the spacecraft flies over the celestial body, R G is the radius of the celestial body, h p is a variable, usually h p has the constraint:

[0066] h p <h p min (3)

[0067] h p min is the minimum value of the height at which the spacecraft flies over the celestial body. From formula (2), it can be known that when other inputs remain unchanged, h p min constrains the maximum value δ max of δ, that is:

[0068]

[0069] Given the initial orbit and target orbit of the spacecraft, the corresponding spacecraft velocity can be obtained according to the known methods in the field. Assuming that the position coordinates of the spacecraft and the gravity assist celestial body in the inertial system are rG =(0, -a G , 0), a G represents the semi-major axis of the orbit of the gravity assist celestial body. Since the gravity assist celestial body is approximately in a circular orbit with an inclination of 0, thus the modulus of r G is always equal to a G , when r G takes other values, through the coordinate rotation between inertial systems, make r G equal to (0, -a G , 0) in the new coordinate system, which does not affect the subsequent calculations, so this assumption is not loss of generality.

[0070] Let the velocity of the gravity assist celestial body be where μ J is the planetary gravitational constant. Let the inclination of the initial orbit be i0 and the velocity vector be Then the relative velocity can be obtained from Equation (5) According to Figure 1 , the velocity of the target orbit is written as Then the initial relative velocity of the spacecraft can be obtained and the relative velocity of the target orbit The included angle Δ is

[0071]

[0072] where v G (t) is the velocity of the gravity assist celestial body, Δ is much larger than the maximum angle δ changed by a single gravity assist max , so it is necessary to achieve step-by-step rotation of to through multiple gravity assists. The total angle turned is Δ, v ∞ is the modulus of Figure 1 as shown.

[0073] Step 2: Calculate all available resonance period ratios during the transfer process according to the initial orbit period and the target orbit period; in Step 2, according to Equations (5) and (4), it can be seen that at least Δ / δ max times of gravity assist are required to achieve the inclination adjustment. However, in the actual gravity assist orbit adjustment process, the period of the spacecraft also needs to be considered so that when it returns to the position again after several orbits, the gravity assist celestial body also passes through this position. That is, the ratio of the period of the spacecraft to the period of the gravity assist celestial body orbiting the planet should be an integer ratio so that the next gravity assist can be continued.

[0074] Given that the initial orbit period is T0, the target orbit period is T f , and the orbit period of the gravity assist celestial body is T G , then obviously the resonance ratio after each gravity assist is between and Therefore, it is necessary to find all possible resonance ratios where M and N are natural numbers, N ≤ N max , N max is a given threshold, indicating that the time interval between two gravity assists is less than N max T G , for selection in the next step of optimization.

[0075] Step 2 is specifically implemented according to the following steps:

[0076] The specific process of obtaining all available resonance period ratios is as follows: Traverse backward starting from N = 1. For a given N and all M in the range 1 ≤ M ≤ N, keep the M that satisfies the constraint . If T f < T0, then the constraint is taken as which does not affect subsequent calculations. The i-th value that satisfies the constraint is denoted as M N,i , M N,i The number of them is denoted as n N , and the corresponding sequence is obtained as Then increment N by 1 and continue to find the M and sequence that satisfy the constraint until N > N max when it ends.

[0077] Finally, merge all the sequences and sort them by size, denoted as

[0078] Step 3: Establish a mixed-integer optimization model for multiple gravity assists;

[0079] Combined with Figure 3 , Step 3 is specifically implemented according to the following steps:

[0080] From Figure 1 it can be seen that when the magnitude of v ∞ remains unchanged, the vector v ∞ can be described by two angular variables. One is the angle α between v ∞ and v G ; the other is the angle β of the projection of v ∞ onto the bottom surface of the cone with v G as the axis and α as the semi-aperture angle relative to the xz plane. According to the following formula:

[0081]

[0082] where v = v ∞ + v G is the spacecraft velocity vector corresponding to any given v ∞ , and α is the semi-major axis of the spacecraft orbit. It can be seen that when α is constant, v∞ ·v G If it remains unchanged, then it also remains unchanged, and the corresponding semi-major axis of the orbit also remains unchanged. Since the orbital period is uniquely determined by the semi-major axis, it shows that for a certain resonance period ratio the corresponding α is obtained by back-calculating from the spacecraft orbital period i ,

[0083]

[0084] where T i is the resonance period ratio corresponding to the spacecraft period, and a and v are intermediate variables in the operation. Then the v before and after each gravity assist ∞ lies on the conical base corresponding to a certain α i ; Therefore, to optimize the resonance orbit of the gravity assist, it is only necessary to optimize the angular change of v ∞ inside the same circle and the jump time between different circles, as shown in Figures 2(a) and 2(b). Figure 2(a) shows that the spacecraft orbital period and the resonance period ratio have changed; Figure 2(b) shows that the spacecraft orbital period and the resonance period ratio remain unchanged. Among them, α i represents the angle between before the i-th gravity assist and v G , and β i represents the angle between the projection of i on the conical base corresponding to α and the z-axis;

[0085] In summary, denoting what is obtained in step 3 as n, the optimization variables include {x } and {β i}, both {x i} and {β i} are 2n-dimensional. Among them, the first n x i are 0-1 variables (can only take 0 or 1), indicating whether the optimal transfer trajectory passes through i the periodic orbit corresponding to, where 0 means not passing through and 1 means passing through; the last n β are real variables, used to locate the time when the spacecraft leaves the current orbital period and transfers to the next orbital period i Let β0 be the angle between the projection of the initial on the conical base corresponding to α0 and the z-axis, and β f be the angle between the projection of the target on the conical base corresponding to α and the z-axis. Then the total transfer time corresponding to a given set of {x f} and {β i} is calculated according to the following process: i} is calculated according to the following process:

[0086] Traverse backward from x1 to find x i i with a value of 1, and the corresponding α i and β i , through α i and β i and the invariant v ∞ can be directly calculated to obtain

[0087]

[0088] That is, x i and β i indicating that the first gravity assist changes the relative velocity of the spacecraft from to The obtained angle value is:

[0089]

[0090] Then, according to the fact that the maximum angle of changing the relative velocity by a single gravity assist is a fixed value δ max , check the constraint. Obviously, when |α i - α0| > δ max , the gravity assist cannot be achieved, and a penalty term needs to be added to the objective function; when δ max ≤ |α i - α0| ≤ Δ i , a single gravity assist can be achieved, and the spacecraft intersects with the position of the gravity assist celestial body again after N i T i time; while when Δ i > δ max , although a single gravity assist cannot be achieved, the orbit transfer process can be split into multiple gravity assists without changing the period. That is, first use a gravity assist to change the resonance period ratio of the spacecraft orbit relative to the gravity assist celestial body orbit, as shown in Fig. 2(a), and then keep the resonance period ratio unchanged and use multiple gravity assists to make the relative velocity vector of the spacecraft reach First transfer to the intermediate vector in the way of Fig. 2(a) Then transfer to

[0091] times in the way of Fig. 2(b) by solving the following system of equations to calculate the vector with an included angle of δ with the x-axis and an included angle of α max with the x-axis, denoted as i , denoted as

[0092]

[0093]

[0094] Then calculate the corresponding β i , and then continue to transfer to The number of gravity assist times is the integer obtained by rounding up, denoted as m i , then from change to The required time is (m i + 1)N i T i ;

[0095] Then continue to traverse backward to find the next x j with a value of 1, and the corresponding α j and β j ; Continue to calculate the time required to change from to in the same way, as well as the constraint penalty term;

[0096] Until the traversal ends at n, assuming that there are k non-zero variables in x i , i = 1,...n, corresponding to k + 1 times of changes, and the total duration and constraint penalty term calculated are used as the optimization objective function:

[0097]

[0098] where Δt j is the time required for each change, k p is the number of times the constraint is not satisfied, and c p is the constant penalty factor;

[0099] In summary, a mixed-integer optimization model with {x i} and {β i} as optimization variables and J as the objective function is obtained.

[0100] Step 4: Use the differential evolution algorithm to optimize and obtain the results.

[0101] Step 4 is specifically implemented according to the following steps:

[0102] The differential evolution algorithm is a commonly used public method in the field. It is used to quickly obtain a set of variable values that optimize the objective function, that is, the optimal solution, by means of a principle similar to population evolution when given optimization variables and an optimization metric function. The process of the differential evolution algorithm is as follows: First, randomly initialize a population containing multiple sets of values. Then, calculate the objective function of the population, randomly select individuals in the population and apply evolutionary operators such as crossover and mutation. Then, recalculate the objective function and continue to apply the evolutionary operators until the maximum number of generations is reached or the optimal solution is obtained.

[0103] Embodiment

[0104] Suppose the spacecraft is in the gravitational field of Jupiter. At the initial moment, its position intersects with Jupiter's moon Ganymede (abbreviated as Ganymede). The position vector is [0, -1070587.469, 0] km. The initial velocity vector of the spacecraft is [0, 0, 3.937909859367544] km / s, and the velocity vector of Ganymede is [10.878, 0, 0] km / s. The gravitational constant of Jupiter is μ J = 126686534.922 km 3 / s 2 and the gravitational constant of Ganymede is μ G = 9887.834 km 3 / s 2 Its own radius is r G = 2634.0 km, the semi-major axis of its orbit around Jupiter is a G = 1070587.469 km, the period is 7.16 days, and the lowest altitude allowed for gravity assist is h p = 50 km. The target orbital velocity vector is [6.73206, 10.80055, 0] km / s. The orbital elements corresponding to the initial orbit and the target orbit of the spacecraft are shown in Table 1.

[0105] Table 1 Initial orbit and target orbit of the spacecraft

[0106] orbit a (km) e i Ω ω M initial 572827.0345 0.8689 90 270 180 180 target 1696094.69611 0.87077 0 0 0 340.253

[0107] The solution steps are as follows

[0108] (1) First, calculate the initial target The modulus of the two is v ∞ = 11.569 km / s, the included angle between the two is Δ = 70.308 degrees, while the maximum angle corresponding to a single gravity assist is only δ max = 3.070 degrees.

[0109] (2) Calculate all available resonance period ratios during the transfer process. The period ratio of the initial spacecraft to Ganymede is 0.3914, and the period ratio of the target orbit to Ganymede is 1.9941. Then, according to the method in step two, the possible sequence of resonance period ratios is obtained as follows: Dimension is 22.

[0110] (3) Establish a 44 - dimensional mixed - integer optimization model.

[0111] (4) Use the differential evolution algorithm to solve and obtain the optimal solution. The total transfer time is 109 times the period of Ganymede, that is, J = 780.119 days, which can transfer the spacecraft from an orbit with an inclination of 90 degrees to an orbit with an inclination of 0 degrees. Figure 5 shows the resonance orbits with different periods during the transfer process. Figures 5(a) and 5(b) show the change process of the relative velocity from to during the gravity - assist process.

[0112] The time consumption of the optimization calculation of the present invention is only 60 seconds; moreover, considering the constraint of the maximum rotation angle of a single gravity - assist, the obtained trajectory can be used as a reference trajectory for mission analysis.

Claims

1. An orbital design method for making large - scale inclination changes using multiple gravity assists, characterized in that, The input information is as follows: Assume that the orbital inclination of the gravity-assist small body around the planet is 0. At a certain moment, the spacecraft has achieved the coincidence of its position with that of the gravity-assist body through orbit control. The position and velocity of the spacecraft and the small body relative to the planet are known. At this moment, the orbital inclination of the spacecraft is i0, and the target inclination to be achieved is i. f , then the inclination adjustment is achieved in the shortest time through multiple gravity assists, and the specific implementation steps are as follows: Step 1: Calculate the change angle of the relative velocity of the gravity assist according to the initial orbit, the target orbit elements, and the position and velocity. Step 2: Calculate all available resonance period ratios during the transfer process according to the initial orbit period and the target orbit period. Step 3: Establish a mixed-integer optimization model for multiple gravity assists. Step 4: Use the differential evolution algorithm to optimize and obtain the results.

2. The orbital design method for large - scale inclination change using multiple gravitational assists according to claim 1, characterized in that, The specific implementation of Step 1 is as follows: According to the known gravity assist approximation model in the field, the gravity assist is an instantaneous process, and there is: where \(t\) is the gravity assist moment, \(t\) - - , \(t\) + + represent the moments before and after the gravity assist instant respectively, \(r(t\) - - ), \(r(t\) + + ) represent the position vectors of the spacecraft before and after the gravity assist instant respectively, \(r\) G G (t) represents the position vector of the gravity assist celestial body, \(v(t\) - - ), \(v(t\) + + ) are the velocity vectors of the spacecraft before and after the gravity assist instant respectively, \(v\) G G (t) is the velocity vector of the gravity assist celestial body, represent the relative velocity vectors of the spacecraft and the celestial body before and after the gravity assist instant respectively, \(v\) ∞ ∞ is the scalar of the relative velocity magnitude, \(\delta\) is the angle between and calculated according to the following formula: where μ G is the gravitational constant of the gravity assist celestial body, h p is the altitude at which the spacecraft flies over the celestial body, R G is the radius of the celestial body, h p is a variable, usually h p has the constraint: h p <h p min (3) h p min is the minimum value of the altitude for flying over the celestial body. From Equation (2), when other inputs remain unchanged, h p min constrains the maximum value δ of δ max , that is: Given the initial orbit and the target orbit of a spacecraft, the corresponding spacecraft velocity can be obtained according to the known methods in the field. Assume that the position coordinates of the spacecraft and the gravity assist celestial body in the inertial system are r G =(0, -a G , 0), where a G represents the semi-major axis of the orbit of the gravity assist celestial body. Since the gravity assist celestial body is approximately in a circular orbit and the inclination is 0, the modulus of r G is always equal to a G . When r G takes other values, the coordinate rotation between inertial systems is used to make r G equal to (0, -a G , 0) in the new coordinate system. Set the velocity of the gravity assist celestial body where μ J is the planetary gravitational constant. Let the inclination of the initial orbit be i0 and the velocity vector be Then the relative velocity can be obtained from Equation (5) The velocity of the target orbit is written as The initial relative velocity of the spacecraft can be obtained and the relative velocity of the target orbit The included angle Δ is where v G (t) is the velocity of the gravity assist celestial body, and Δ is much larger than the maximum angle δ changed by a single gravity assist max , so it is necessary to achieve the step-by-step rotation of to through multiple gravity assists. The total angle turned is Δ, and v ∞ is the modulus of.

3. The orbital design method for large-range inclination change using multiple gravity assists according to claim 2, wherein In Step 2, Given that the initial orbital period is \(T_0\) and the target orbital period is \(T\). f The orbital period of the gravity assist celestial body is \(T\). G Then it is obvious that the resonance ratio after each gravity assist is between and Therefore, it is necessary to find all possible resonance ratios where \(M\) and \(N\) are natural numbers, \(N\leq N\). max \(N\). max is a given threshold, indicating that the time interval between two gravity assists is less than \(N\). max \(T\). G for selection in the next step of optimization.

4. The orbital design method for large-scale inclination change using multiple gravity assists according to claim 3, characterized in that The specific implementation of Step 2 is as follows: The specific process for obtaining all available resonance period ratios is as follows: Start traversing backward from N = 1. For a given N and all M in the range 1 ≤ M ≤ N, retain the M that satisfies the constraint If f T < T0, then the constraint is taken as The i-th value that satisfies the constraint is denoted as M N,i , and the number of M N,i is denoted as n N , obtaining the corresponding sequence as Then increment N by 1 and continue to search for the M that satisfies the constraint and the sequence; until N > N max when it ends. Finally, merge all the sequences and sort them by size, denoted as 5. The orbital design method for large-scale inclination change using multiple gravity assists according to claim 4, characterized in that, The specific implementation of Step 3 is as follows: According to the following formula: where \(v = v\) ∞ + v G is the spacecraft velocity vector corresponding to any given \(v\), \(\alpha\) is the semi-major axis of the spacecraft orbit. When \(\alpha\) is constant, \(v\) ∞ · \(v\) ∞ is invariant, then G is also invariant, and the corresponding semi-major axis of the orbit also remains unchanged. The orbital period is uniquely determined by the semi-major axis, indicating that for a certain resonance period ratio the corresponding \(\alpha\) can be obtained by back-calculating from the spacecraft orbital period , i ​ Among them, T i is the resonance period ratio corresponding to the spacecraft period, and a and v are intermediate variables for calculation; then the v before and after each gravity assist ∞ is on the conical bottom surface corresponding to a certain α i ; among them, α i represents the angle between before the i-th gravity assist and v G , and β i represents the angle between the projection of on the conical bottom surface corresponding to α i and the z-axis; In summary, denote the result obtained in step 3 as n. Then the optimization variables include {x } and {β i}, both {x i} and {β i} are 2n-dimensional. Among them, the first n x i are 0-1 variables, indicating whether the optimal transfer trajectory passes through the corresponding periodic orbit, where 0 means not passing through and 1 means passing through; the last n β i are real variables, used to locate the i when the spacecraft leaves the current orbit period and transfers to the next orbit period. Let β0 be the initial the angle between the projection of the cone bottom corresponding to α0 on the z-axis, and β f is the target the angle between the projection of the cone bottom corresponding to α f on the z-axis. Then, for a given set of {x i} and {β i}, the total transfer time is calculated according to the following process: Traverse backward from x1 to find x i i with a value of 1, and the corresponding α i and β i , through α i and β i and the invariant v ∞ can be directly calculated to obtain i.e., x i and β i The first gravity assist represented changes the relative velocity of the spacecraft from to The obtained angle value is: Then, the maximum angle for changing the relative velocity by a single gravity assist is a fixed value δ max , and the constraint is verified. Obviously, when |α i - α0| > δ max , the gravity assist cannot be achieved, and a penalty term needs to be added to the objective function; when δ max ≤ |α i - α0| ≤ Δ i , a single gravity assist can be achieved, and the spacecraft intersects the position of the gravity assist celestial body again after N i T i time; and when Δ i > δ max , the orbit transfer process is split into multiple gravity assists without changing the period. That is, first, a single gravity assist is used to change the resonance period ratio of the spacecraft's orbit relative to the orbit of the gravity assist celestial body, and then, while keeping the resonance period ratio unchanged, multiple gravity assists are used to make the relative velocity vector of the spacecraft reach Calculate the included angle δ with by solving the following system of equations max , and the vector with an included angle of α with the x-axis i is denoted as Then calculate the corresponding β i , then continue to transfer to The number of gravity assist times is The integer obtained by rounding up, denoted as m i , then from Change to The required time is (m i + 1)N i T i ; Then continue to traverse backward to find the next x j with a value of 1, and the corresponding α j and β j ; continue to calculate in the same way the time required from changing to as well as the constraint penalty term; Until the traversal ends at n, assume x i , among i = 1,... n, there are k non-zero variables, corresponding to k + 1 times of change. The total duration and the constraint penalty term calculated are used as the optimization objective function: where, Δt j is the time required for each change, k p is the number of times the constraint is not satisfied, and c p is a constant penalty factor; In summary, a mixed-integer optimization model with {x i} and {β i} as optimization variables and J as the objective function is obtained.

Images