A Mission Planning Method for In-Orbit Refueling in the Same Orbit Plane of GEO Orbit

Abstract

A method for planning in-orbit refueling tasks for coplanar GEO orbits. The present invention relates to a method for planning in-orbit refueling tasks for coplanar GEO orbits. The purpose of the present invention is to solve the problem that the existing in-orbit round-trip refueling for coplanar GEO orbits often adopts a single-layer or two-layer optimization strategy, which will lead to falling into local optimal solutions during the optimization process, and the convergence speed is not fast enough, making it difficult to find the optimal solution, resulting in the fuel consumption of the service satellite when traveling back and forth between the service station and the target satellite for refueling tasks not being able to reach the optimal value. The process is as follows: Step 1, establish a model for planning in-orbit refueling tasks for coplanar GEO orbits; Step 2, solve the model for planning in-orbit refueling tasks for coplanar GEO orbits according to the CGAPB three-layer optimization algorithm to obtain the optimal service sequence of the service satellite, the optimal service time sequence, and the optimal time node sequence for the service satellite to return to the service station. The present invention is used in the field of aerospace technology.

Description

Technical Field

[0001] The present invention relates to a method for planning an on-orbit refueling mission for a GEO orbit coplanarity. Background Art

[0002] As space technology becomes more and more mature, human space activities are increasing year by year, and the number of satellites in orbit is gradually increasing. According to data from the Celestrak website, as of April 30, 2022, there are 5,465 satellites in orbit, including 565 GEO satellites. In order to maintain their position, many GEO satellites need to perform orbital maneuvers, which consumes fuel. According to statistics1 (1 Zhai Guang, Zhang Jingrui, Zhou Zhicheng. Research Progress on In-Orbit Life Extension Technology for Geostationary Satellites [J]. Journal of Astronautics, 2012, 33(07): 849-859.), the number of GEO spacecraft that lost their function due to fuel exhaustion between 2008 and 2010 accounted for 85.4%, and with the improvement of the reliability design of hardware facilities, the number and proportion of fuel exhaustion will further increase. Due to the limited resources of the GEO orbit and the high launch cost of the GEO orbit, most of the GEO orbit satellites are high-value satellites. In order to prevent high-value GEO satellites from losing their use value prematurely and to ensure that they continue to perform their due functions, on-orbit refueling technology came into being.

[0003] Considering that it is too costly to launch a refueling spacecraft from the ground to the GEO orbit to refuel satellites, scholars have proposed the concept of a service station, that is, to deploy a service station in the GEO orbit in advance. The service satellites stationed in the service station can provide refueling services for GEO orbit satellites.

[0004] In order to use the fuel stored in the service station to replenish the GEO target satellite as much as possible, the fuel consumed by the orbital maneuver of the service satellite is required to be as little as possible. In application, when refueling a single target satellite, the proportion of fuel consumed by the orbital transfer of the service satellite is very high, so a service satellite often carries enough fuel to sequentially refuel multiple target satellites within a specified time. This requires planning the overall task, and planning the service order, corresponding service time, the time node for the service satellite to return to the fuel station, and the speed pulse of each orbit change under the constraints of time and the fuel carried by the service satellite. Taking into account the characteristics of the GEO orbit, all orbits are approximately considered to be circular orbits during the implementation of this patent. Summary of the invention

[0005] The purpose of the present invention is to solve the problem that the existing in-orbit round-trip refueling in the same GEO orbit plane often adopts a single-layer or two-layer optimization strategy, which may lead to getting stuck in local optimal solutions during the optimization process, and the convergence speed is not fast enough, making it difficult to find the optimal solution, resulting in the fuel consumption not being optimal when the service satellite travels back and forth between the service station and the target satellite for refueling tasks. Therefore, a method for planning in-orbit refueling tasks in the same GEO orbit plane is proposed.

[0006] The specific process of a method for planning in-orbit refueling tasks in the same GEO orbit plane is as follows:

[0007] Step 1: Establish a model for planning in-orbit refueling tasks in the same GEO orbit plane;

[0008] Step 2: Solve the model for planning in-orbit refueling tasks in the same GEO orbit plane according to the CGAPB three-layer optimization algorithm to obtain the optimal service sequence of the service satellite, the optimal service time sequence, and the optimal time node sequence for the service satellite to return to the service station.

[0009] The beneficial effects of the present invention are as follows:

[0010] The present invention proposes a CGAPB three-layer optimization algorithm for in-orbit round-trip refueling in the same GEO orbit plane. In the three-layer optimization algorithm, the first layer optimizes the service sequence X of the service satellite, the second layer optimizes the service time T of the service satellite, and the third layer optimizes the time node S for the service satellite to return to the service station, which speeds up the convergence rate, reduces the calculation amount, and is easy to find the optimal solution. It solves the problem that the current application of the two-layer algorithm to solve this problem has a slow convergence and it is difficult to obtain the optimal solution. Considering the characteristics of the CGAPB model, good results can also be achieved for in-orbit round-trip refueling in different GEO orbits. Description of the Drawings

[0011] Figure 1 It is the flowchart of the three-layer optimization algorithm of the present invention;

[0012] Figure 2 It is the crossover operation diagram;

[0013] Figure 3 It is the duplicate removal operation diagram;

[0014] Figure 4 It is the mutation operation diagram;

[0015] Figure 5 It is the phasing transfer diagram;

[0016] Figure 6 It is the phasing transfer characteristic curve diagram;

[0017] Figure 7 It is the refueling model diagram for traveling back and forth to the service station;

[0018] Figure 8It is the population evolution curve graph of the two - layer optimization algorithm;

[0019] Figure 9 It is the population evolution curve graph of the CGAPB algorithm. Specific implementation manners

[0020] Specific implementation manner one: The specific process of a method for planning the in - orbit refueling task of co - planar GEO orbits in this implementation manner is as follows:

[0021] Task scenario

[0022] There are currently multiple co - planar satellites in the GEO orbit. Due to years of service, they lack fuel and need to be refueled within a specified time. There is a service station on the same orbit, and a certain amount of fuel is stored in the service station. There is a service satellite in the service station that shuttles between the service station and the target satellites to perform the refueling task. Since the fuel stored in the service station is extremely precious, when meeting the fuel requirements of all target satellites, it is required that the fuel consumed by the transfer orbit of the service satellite is the least. It is necessary to give the service sequence for this task, the service time for each target satellite, when the service satellite returns to the service station for resupply, and the mass of fuel carried by the service satellite each time it departs from the service station. This patent is to give a scheme with the least fuel consumption through an algorithm when dealing with this kind of task.

[0023] Key steps of the invention patent

[0024] For the GEO co - planar in - orbit refueling problem faced by this patent, it is very difficult to obtain the optimal solution by using traditional single intelligent algorithms such as genetic algorithms and particle swarm algorithms, or various double - layer combined intelligent optimization algorithms. Moreover, the convergence speed is very slow and the computational amount is very large. This patent proposes a three - layer optimization algorithm based on chaotic genetic particle swarm branch - and - bound (CGAPB), which well solves this problem. In the three - layer optimization algorithm, the first layer optimizes the service sequence X of the service satellite, the second layer optimizes the service time T of the service satellite, and the third layer optimizes the time node S when the service satellite returns to the service station, which speeds up the convergence rate and reduces the computational amount.

[0025] Step 1: Establish a planning model for the in - orbit refueling task of co - planar GEO orbits;

[0026] Step 2: Solve the planning model for the in - orbit refueling task of co - planar GEO orbits according to the CGAPB three - layer optimization algorithm to obtain the optimal service sequence sequence of the service satellite, the optimal service time sequence, and the optimal time node sequence when the service satellite returns to the service station.

[0027] Specific implementation manner two: The difference between this implementation manner and specific implementation manner one is that in step 1, a planning model for the in - orbit refueling task of co - planar GEO orbits is established; the specific process is as follows:

[0028] Establish a mission planning model for in-orbit refueling with GEO orbits coplanar;

[0029]

[0030] Among them, n is the total number of target stars, and x i and x j are the numbers of the i-th and j-th target stars respectively;

[0031] t i is the time consumed for the service star to refuel the i-th target star;

[0032] s i is the decision flag on whether the service star returns to the service station after refueling the i-th target star, 0 means not returning, and 1 means returning;

[0033] s n is the decision flag for the service star to return to the service station after refueling the n-th target star;

[0034] f i is the mass of fuel carried by the service star when departing from the service station to refuel the i-th target star, in kg. If it does not depart from the service station to refuel the i-th target star, it is 0;

[0035] f is the mass of fuel consumed by the service star for the entire refueling mission, in kg;

[0036] C is the upper limit of fuel carried by the service star, in kg;

[0037] T max is the constraint on the total mission duration.

[0038] x i is i; x j is j.

[0039] Other steps and parameters are the same as those in the specific implementation method 1.

[0040] Specific implementation method 3: The difference between this implementation method and the first or second specific implementation method is that in step 2, the GEO orbit coplanar in-orbit refueling mission planning model is solved according to the CGAPB three-layer optimization algorithm to obtain the optimal service sequence of the service star, the optimal service time sequence, and the optimal time node sequence for the service star to return to the service station; the specific process is as follows:

[0041] Step 2-1. Convert the GEO orbit coplanar in-orbit refueling mission planning model into the following:

[0042] find X best =[x1,x2,…,x n

[0043] ​

[0044] where n is the total number of target stars, and x i , x j is the number of a specific target star;

[0045] X best is the optimal service order sequence of the service star;

[0046] find means to search for;

[0047] T best is the optimal service time sequence corresponding to the service sequence X best of the service star, with the unit of day;

[0048] t i is the time consumed by the service star for refueling the i-th target star;

[0049] T max is the total mission duration constraint, with the unit of day;

[0050] S best is the optimal time node sequence for the service star to return to the service station, where 0 means not to return and 1 means to return;

[0051] s i is the decision flag indicating whether the service star returns to the service station after refueling the i-th target star, where 0 means not to return and 1 means to return;

[0052] s n is the decision flag for the service star to return to the service station after refueling the n-th target star;

[0053] f i is the fuel mass carried by the service star when departing from the service station to refuel the i-th target star, with the unit of kg. If not departing from the service station to refuel the i-th target star, it is 0;

[0054] f is the fuel mass consumed by the service star for the entire refueling mission, with the unit of kg;

[0055] C is the upper limit of the fuel carried by the service star, with the unit of kg;

[0056] Step 2: Solve the on-orbit refueling mission planning model for GEO orbit coplanar according to the CGAPB three-layer optimization algorithm to obtain the optimal service order sequence of the service star, the optimal service time sequence, and the optimal time node sequence for the service star to return to the service station.

[0057] Other steps and parameters are the same as those in the specific implementation method 1 or 2.

[0058] Embodiment 4: The difference between this embodiment and any one of Embodiments 1 to 3 is that in Step 22, the CGAPB three-layer optimization algorithm is used to solve the GEO orbit coplanar in-orbit refueling mission planning model to obtain the optimal service order sequence of the service satellites, the optimal service time sequence, and the optimal time node sequence for the service satellites to return to the service station; the specific process is as follows:

[0059] Step 221: Use the chaotic genetic algorithm to optimize and solve the service sequence;

[0060] Step 222: Use the particle swarm optimization (PSO) algorithm to optimize the service time;

[0061] Step 223: Use the branch and bound (B&B) algorithm to optimize the time node for returning to the service station.

[0062] Other steps and parameters are the same as those in any one of Embodiments 1 to 3.

[0063] Embodiment 5: The difference between this embodiment and any one of Embodiments 1 to 4 is that in Step 221, the chaotic genetic algorithm is used to optimize and solve the service sequence; the specific process is as follows:

[0064] Step 2211: Initialize the population; the specific process is as follows:

[0065] Use the mapping method from integer permutation to decimal to perform chaotic initialization on the service order sequence set X of the service satellites;

[0066] X = [X1, X2, …, X i′ , …, X NP ; the optimal service order sequence X best in Step 1 is one of the sequences in this set;

[0067] X i′ = [x1, x2, …, x i , …, x n ;

[0068] Among them, NP is the population size (selected artificially), that is, the number of service order sequences of the service satellites in the service order sequence set of the service satellites; X i′ is the i'-th service order sequence in the service order sequence set of the service satellites; x i is the number of the i-th target satellite;

[0069] Step 2212: Calculate the fitness value; the specific process is as follows:

[0070] Assign each X i′ to Step 222 to obtain the service time T i′ corresponding to each X i′, the time node sequence S of the service stars returning to the service station i′ and the fitness value F i′ ;

[0071] Compare the fitness values to obtain the optimal individual in this generation, update the optimal individual, and record it as the service order sequence, service time, time node sequence of the service stars returning to the service station, and fitness value;

[0072] Step 2213. Perform a selection operation on the sequence X = [X1, X2, …, X i′ , …, X NP : Adopt the roulette wheel strategy;

[0073] Step 2214. Perform a crossover operation. As Figure 2 shown, randomly select the crossover position, and crossover the genes after the selected position to obtain Figure 2 the offspring;

[0074] Then perform a de-duplication operation: from left to right, replace the repeated genes with those that have not appeared, as Figure 3 shown;

[0075] Step 2215. Mutation operation: Randomly select two genes in the chromosome and exchange them, as Figure 4 shown;

[0076] Step 2216. Perform chaotic perturbation;

[0077] Step 2217. Determine whether the iteration requirement is met. If it is met, end and obtain the optimal service order sequence X best , the corresponding optimal service time T best , the optimal time node sequence S of the service stars returning to the service station best and the optimal fitness value F best ; If the iteration requirement is not met, execute Step 2212.

[0078] Other steps and parameters are the same as those in any one of the specific embodiments 1 to 4.

[0079] Specific embodiment 6: The difference between this embodiment and any one of the specific embodiments 1 to 5 is that the chaotic initialization service star service order sequence set X is performed by using the mapping method from integer arrangement to decimal;

[0080] X = [X1, X2, …, X i′ , …, X NP ; The service order sequence X best of the optimal service star in Step 1 is one of the set;

[0081] X i′ = [x1, x2, …, xi , …, x n ;

[0082] Among them, NP is the population size (selected artificially), that is, the number of service star service order sequences in the service star service order sequence set; X i′ is the i'-th service order sequence in the service star service order sequence set; x i is the number of the i-th target star;

[0083] The specific method is as follows:

[0084] Suppose there are n target stars in total, and there are n! permutation ways;

[0085] Set z0 = 0.2 to obtain the initial value;

[0086] Select the motion equation z i″+1 = 4z i″ (1 - z i″ ), i″ = 1, …, n, …, 3NP. Map the decimal z i″ to the permutation X j , j = 1, 2, …, n! (3NP X j constitute X 3NP ). The specific mapping method is as follows:

[0087] Arrange the n! permutation ways in ascending order, then z i represents the round(z i × n!) -th permutation way;

[0088] where round is rounding;

[0089] Randomly select the time series T 3NP corresponding to the 3NP service sequence populations X 3NP , and execute step two two four to obtain the corresponding node S 3NP returning to the service station and the corresponding fitness value F 3NP ; Take the first NP service order sequences X = [X1, X2, …, X i′ , …, X NP corresponding to the smaller fitness values as the initial population.

[0090] Other steps and parameters are the same as one of the specific implementation manners one to five.

[0091] Specific implementation manner seven: The difference between this implementation manner and one of the specific implementation manners one to six is that chaotic perturbation is performed in step two two one six; The specific process is as follows:

[0092] First, X i′The decimal z mapped to 0 to 1 (including 0 and 1) i′ , and the specific mapping method is as follows:

[0093] X after the sequence undergoes the mutation operation of step two-two-one-five i′ Arrange them in ascending order as n! ones, where the N Xi′ th (the N Xi′ th is one of the n! ones) arranged and mapped decimal is

[0094] Set the chaotic perturbation coefficient k is the number of iterations;

[0095] Set the perturbation amount Δz = rand × α × (-1) randi ×p;

[0096] rand is a random number from 0 to 1 (0 ≤ rand ≤ 1), randi is 0 or 1, and p needs to be determined according to the length of n;

[0097] Add the perturbation amount Δz to z i′ to obtain the new decimal z i′ ';

[0098] According to the mapping method of step two-two-one-one, map to the sequence X i′ ', assign the sequence X i′ ' to step two-two-two to obtain the service time T i′ ' corresponding to each X i′ ', the time node sequence S i′ ' when the corresponding service star returns to the service station, and the corresponding fitness value F i′ ';

[0099] If a smaller fitness value is obtained, then replace X i′ ' with X i′ to obtain the new sequence X = [X1, X2,..., X i′ ',..., X NP .

[0100] Other steps and parameters are the same as those in any one of the specific embodiments one to six.

[0101] Specific embodiment eight: The difference between this embodiment and any one of the specific embodiments one to seven is that in step two-two-two, the particle swarm optimization (PSO) algorithm is used to optimize the service time T best ; The specific process is as follows:

[0102] Step two-two-one-one, initialize the particle swarm parameters; The specific process is as follows:

[0103] Randomly generate a particle population T = [T1, T2,..., T i′,…,T NP , T i′ =[t1, t2, …, t i , …, t n , NP is the number of particle populations (selected manually), satisfying the constraint t i ∈[1, T max , sum(t i ) ≤ T max ;

[0104] Randomly generate the particle velocity V;

[0105] Step Two Two Two Two: Combine T1, T2, …, T i′ , …, T NP with the service order in Step Two Two One (X i′ or X i′ ′) and assign them to Step Two Two Three, calculate the time node sequence when the service star returns to the service station and the corresponding fitness value (fuel consumption value);

[0106] Record the time series corresponding to the minimum fitness value as the optimal service time;

[0107] Step Two Two Two Three: Update the velocity V and the randomly generated particle population (position) T NP , where the inertia weight:

[0108]

[0109] where iter_max is the maximum number of iterations, w min represents the minimum inertia weight, w max represents the maximum inertia weight, and iter represents the current number of iterations;

[0110] Step Two Two Two Four: If the number of iterations is satisfied, return F min , and the corresponding time series; if not, execute Step Two Two Two Two.

[0111] Other steps and parameters are the same as those in any one of the first to seventh specific embodiments.

[0112] Specific Embodiment Nine: The difference between this embodiment and any one of the first to eighth specific embodiments is that in Step Two Two Three, the branch and bound (B&B) algorithm is used to optimize the time node for returning to the service station; the specific process is as follows:

[0113] Step Two Two Three One: Initialize the parameters of the branch and bound (B&B) algorithm; the specific process is as follows:

[0114] Initialize the matrices ind0 = [0] and ind1 = [1], and the combined matrix Ind = [ind0; ind1]. Since the last digit of the return service station node must be 1, only the first n - 1 digits need to be solved, where n is the number of target stars;

[0115] Step 2232: According to T i′ and the service order in Step 221, combine with (X i′ or X i′ ′) to optimize the time nodes of the return service station and the corresponding fitness values.

[0116] Other steps and parameters are the same as those in any one of Embodiments 1 to 8.

[0117] Embodiment 10: The difference between this embodiment and any one of Embodiments 1 to 9 is that in Step 2232, according to T i′ and the service order in Step 221, combine with (X i′ or X i′ ′) to optimize the time nodes of the return service station and the corresponding fitness values; the specific steps are as follows:

[0118] Step 1: Obtain the number of rows and columns of Ind (i.e., obtain which target has been processed); if the number of columns is less than n - 1, execute Step 2, otherwise execute Step 3;

[0119] Step 2: Fill 0 after each row in the matrix ind0 and fill 1 after each row in the matrix ind1 respectively. According to T i′ and the service order in Step 221, combine with (X i′ or X i′ ′) to calculate the fitness value of each row in the matrix ind1 (the fitness value is obtained from Formulas (9), (10), (11), (12), and in Formula (9) can be obtained from Formula (6), and Formula (6) is obtained from Formulas (1), (2), (3), (4), (5)), compare the calculated fitness values of each row, and the row with the minimum fitness value among all rows is the partial optimal solution in ind1 (row vector, i.e., the time node of the return service station) I opt , and update Ind = [ind0; I opt , and execute Step 1;

[0120] Step 3: Fill 1 after each row in the matrix ind1. According to T i′ and the service order in Step 221, combine with (X i′ or X i′ ′) to calculate the fitness value of each row in the matrix ind1 (the fitness value is obtained from Formulas (9), (10), (11), (12), and in Formula (9) It can be obtained from formula (6), and formula (6) is obtained from formulas (1), (2), (3), (4), and (5). Compare the fitness values calculated for each row. The row with the smallest fitness value among all rows is the optimal complete solution in ind1 (a row vector, i.e., the time node for returning to the service station).

[0121] Step4: Obtain the optimal complete solution and the corresponding fitness value.

[0122] Other steps and parameters are the same as those in any one of the specific implementation manners one to nine.

[0123] In summary, the flow chart of the CGAPB three-layer optimization algorithm Figure 1 is shown as follows.

[0124] Related technologies

[0125] Phase modulation transfer and its characteristics

[0126] Phase modulation transfer can adjust the phase angle within the same orbital plane and change the position of the spacecraft. For phase angle adjustment using a double-pulse method within the same orbital plane, phase modulation transfer is more fuel-efficient than Hohmann transfer and Lambert transfer. Theoretically, phase modulation transfer does not have the constraints of the transfer time of half an orbital period required by Hohmann transfer and the transfer time of Lambert, and has more advantages in coplanar phase modulation.

[0127] The schematic diagram of phase modulation transfer is as Figure 5 shown. For two spacecraft at different positions on the same orbit to conduct rendezvous, B is the service satellite, and A and C are the target satellites. Initially, all three satellites are on the same circular orbit ① with a semi-major axis of r1. Then, the periods, angular velocities, and velocities of spacecraft A, B, and C are respectively:

[0128]

[0129] where G represents the gravitational constant and M represents the mass of the Earth;

[0130] Assume that it is required for the service satellite B to rendezvous with the spacecraft A at the initial point of the service satellite B after N orbits. Then, the transfer time required for the service satellite B is:

[0131]

[0132] where θ1 represents the initial phase difference between the service satellite and the target satellite;

[0133] Then, the transfer orbit period and semi-major axis after the maneuver of spacecraft B are:

[0134]

[0135] Thus, the eccentricity and the perigee length of the transfer orbit of spacecraft B can be further obtained as follows:

[0136]

[0137] Then, the magnitude of the velocity at the apogee (i.e., the initial point of B) of the transfer orbit of spacecraft B is:

[0138]

[0139] Then, the magnitudes of the two pulses when spacecraft B rendezvouses with A are respectively:

[0140]

[0141] As known from formula (6), the magnitude of the velocity pulse required for phasing transfer is related to the number of transfer loops and the phase angle to be adjusted. The longer the transfer duration and the smaller the phase angle to be adjusted, the less fuel is consumed. The specific relationship is as Figure 6 shown.

[0142] As Figure 6 shown, the phasing parameters with phase angle differences of 40 degrees, 80 degrees, 120 degrees, and 160 degrees are respectively selected for simulation from 3 loops (3 days) to 28 loops (28 days). It can be seen that as the number of phasing loops increases, the fuel consumed during phasing transfer gradually decreases regardless of the phase angle difference, and the decreasing amplitude changes from fast to slow. In theory, when the phasing time is long enough, the required fuel tends to 0, but in practice, there must be a time limit.

[0143] Refueling mode of round-trip service station

[0144] As Figure 7 shown, the current target stars numbered 1 - 8 lack fuel and need on-orbit refueling. The service star departs from the service station for sequential service. Since the fuel carried by the service star is limited, it is necessary to return to the service station to replenish fuel. Figure 7 There are 3 service paths in total as shown in

[0145] After the service star carries enough fuel at the service station and starts to serve stars ④ and ② in the order of path 1, due to insufficient fuel carried, it needs to return to the fuel station for replenishment. Then, it refuels stars ⑦, ⑥, ⑤, ③, and ① in sequence according to path 2. After completion, it returns to the fuel station to replenish fuel, then departs to refuel star ③, and finally returns to the service station.

[0146] Fuel consumption model

[0147] Key Elements and Establishment of the Model

[0148] The optimization objective of the model is to minimize the fuel consumed for the orbital transfer of the servicing satellite after servicing all target stars, that is, when meeting the fuel requirements of all target stars, the total mass of fuel replenished by the servicing satellite at the servicing station is minimized:

[0149]

[0150] Among them, n is the total number of target stars, and f i is the mass of fuel carried by the servicing satellite when departing from the servicing station to refuel the i-th target star, with the unit of kg. If it does not depart from the servicing station to refuel the i-th target star, it is 0; f is the mass of fuel consumed by the servicing satellite for the entire refueling mission, with the unit of kg;

[0151] The design variable X = [x1, x2, …, x n is the servicing sequence, T = [t1, t2, …, t n is the corresponding servicing time, and S = [s1, s2, …, s n is the decision sequence for whether to return to the servicing station, where 0 means not returning and 1 means returning. Then f i = F(X, T, S), and the optimization index becomes:

[0152]

[0153] Among them, T max is the constraint on the total mission duration, and C is the upper limit of the fuel carried by the servicing satellite each time.

[0154] Calculation of the Objective Function Value

[0155] Assume that the servicing satellite leaves the servicing station for the -th time and will refuel n target stars. Then a total of 2(n + 1) velocity pulses are required. Let be the velocity increment required for the servicing satellite to leave the servicing station for the -th time to visit the i-th target star x i , which can be obtained from formula (6); be the remaining mass of the servicing satellite after refueling the i-th target star x when leaving the servicing station for the i -th time. Then there is:

[0156]

[0157] Among them, is the total mass of the servicing satellite when leaving the servicing station for the -th time, I sp is the specific impulse of the satellite, and g0 is the acceleration of gravity on the earth. For the target star x i The mass of the required fuel, After the service star completes this refueling mission and returns to the service station, there is:

[0158]

[0159] To minimize the fuel consumption of the mission, assume that the fuel is just consumed when the service star returns to the service station after completing the mission; let the dry weight of the satellite be M dry , there is ( is When i takes n in, the mass of the service star when it returns to the service station for the last time after completing the mission);

[0160] Then, after determining the service order and service time, the formula for calculating the initial mass of the service star can be obtained:

[0161]

[0162] Then the mass of the fuel that the service star needs to carry is:

[0163]

[0164] At the same time, the mass of the fuel carried each time needs to meet the constraint condition f i < C.

[0165] The following embodiments are used to verify the beneficial effects of the present invention:

[0166] Embodiment 1:

[0167] Relevant parameters and simulation results

[0168] Relevant parameters

[0169] Select 10 currently orbiting GEO satellites according to the data on the UCS Satellite Database website, and give the mass of the required fuel. The data is shown in Table 1.

[0170] Table 1 Orbital parameters of the target star to be refueled

[0171]

[0172] The algorithm parameters are shown in Table 2.

[0173] Table 2 Algorithm parameters

[0174]

[0175] To illustrate the effectiveness of the proposed CGAPB algorithm of the present invention, a widely used two-layer optimization algorithm is designed for comparison. The genetic algorithm is used in the outer layer to optimize the service sequence X and the corresponding time series Y, and the B&B algorithm is used in the inner layer to optimize the time node S for returning to the service station. The algorithm parameters are shown in Table 3.

[0176] Table 3 Algorithm Parameters

[0177]

[0178]

[0179] Other satellite parameters are the same as the data in Table 1.

[0180] Scheme Results

[0181] Through calculation, the final scheme is as follows: The results obtained by the two-layer and CGAPB three-layer optimization algorithms are shown in Table 4. The population change process of the two-layer algorithm is as Figure 8 shown, and the population change process of the CGAPB algorithm is as Figure 9 shown.

[0182] Table 4 Results of Two-Layer and CGAPB Optimization Algorithms

[0183] Variable symbol Two-layer optimization result Three-layer optimization result X [1,8,4,3,6,5,7,2,10,9] [1,2,10,9,7,8,6,5,4,3] T [6,12,9,12,10,7,12,13,5,14] [8,14,7,5,15,8,8,15,9,11] S [0,1,0,1,0,0,1,0,0,1] [1,0,0,0,1,0,0,1,0,1] Sf [729.553,824.884,853.976,721.505] [293.059,992.966,998.511,826.904] f 3129.918 3111.440

[0184] It can be seen from the above results that in the process of population evolution, the two-layer optimization algorithm is obviously trapped in a local optimum, and its convergence rate is not fast enough. The CGAPB three-layer optimization algorithm adopted in this patent converges faster than the two-layer optimization algorithm, and the obtained solution is better. It can be seen from Figure 9 that when the population evolves to the 31st generation, all the population individuals have evolved into the optimal individuals, indicating that the convergence speed of this algorithm is very fast and the optimal solution is found, which is very effective in solving the problem of on-orbit refueling mission planning for GEO orbit coplanarity.

[0185] The present invention can also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art can make various corresponding changes and deformations according to the present invention, but these corresponding changes and deformations should all fall within the protection scope of the appended claims of the present invention.

Claims

1. A method for planning in-orbit refueling tasks for GEO orbit coplanarity, characterized in that: The specific process of the method is as follows: Step 1: Establish a mission planning model for in-orbit refueling with GEO orbits coplanar; Step 2: Solve the mission planning model for in-orbit refueling with GEO orbits coplanar according to the CGAPB three-layer optimization algorithm to obtain the optimal service order sequence of the service satellites, the optimal service time sequence, and the optimal time node sequence for the service satellites to return to the service station; In Step 1, the mission planning model for in-orbit refueling with GEO orbits coplanar is established. The specific process is as follows: Establish a mission planning model for in-orbit refueling with GEO orbits coplanar; where n is the total number of target stars, and x i , x j are the numbers of the i-th and j-th target stars respectively; t i The time consumed for refueling the i-th target star by the service star; s i It is the decision flag on whether to return to the service station after refueling the i-th target star by the service star. 0 means not to return, and 1 means to return; s n It is the decision flag returned to the service station after refueling the nth target star by the service star; f i The mass of the fuel carried by the service satellite for refueling the i-th target satellite when departing from the service station, with the unit of kg. If it does not depart from the service station for refueling the i-th target satellite, it is 0. f is the mass of fuel consumed by the service satellites for the entire refueling mission, with the unit of kg; C is the upper limit of fuel carried by the service satellites, with the unit of kg; T max is the total task duration constraint; In Step 2, the mission planning model for in-orbit refueling with GEO orbits coplanar is solved according to the CGAPB three-layer optimization algorithm to obtain the optimal service order sequence of the service satellites, the optimal service time sequence, and the optimal time node sequence for the service satellites to return to the service station. The specific process is as follows: Step 2-1: Convert the mission planning model for in-orbit refueling with GEO orbits coplanar into the following: find X best =[x1,x2,…,x n ​ where n is the total number of target stars, and x i , x j is the number of a specific target star; X best is the service order sequence of the optimal service stars; find means to search for; T best is the optimal service time sequence corresponding to the service star according to the service sequence X best in days; t i Time consumed for refueling the i-th target star by the service star; T max is the total task duration constraint, in days; S best It is the optimal time node sequence for the service star to return to the service station, where 0 means not returning and 1 means returning; s i It is the decision flag on whether to return to the service station after refueling the i-th target star by the service star. 0 means not to return, and 1 means to return; s n The service star returns to the service station decision flag after refueling the nth target star; f i The mass of the fuel carried by the service satellite for refueling the i-th target satellite when departing from the service station, with the unit of kg. If it does not depart from the service station for refueling the i-th target satellite, it is 0; f is the mass of fuel consumed by the service satellites for the entire refueling mission, with the unit of kg; C is the upper limit of fuel carried by the service satellites, with the unit of kg; Step 2-2: Solve the mission planning model for in-orbit refueling with GEO orbits coplanar according to the CGAPB three-layer optimization algorithm to obtain the optimal service order sequence of the service satellites, the optimal service time sequence, and the optimal time node sequence for the service satellites to return to the service station; In Step 2-2, the mission planning model for in-orbit refueling with GEO orbits coplanar is solved according to the CGAPB three-layer optimization algorithm to obtain the optimal service order sequence of the service satellites, the optimal service time sequence, and the optimal time node sequence for the service satellites to return to the service station. The specific process is as follows: Step 2-2-1: Use the chaotic genetic algorithm to optimize and solve the service sequence; Step 2-2-2: Use the particle swarm algorithm to optimize the service time; Step 2-2-3: Use the branch and bound algorithm to optimize the time nodes for returning to the service station; In Step 2-2-1, the chaotic genetic algorithm is used to optimize and solve the service sequence. The specific process is as follows: Step 2-2-1-1: Initialize the population. The specific process is as follows: Use the mapping method from integer permutation to decimal to perform chaotic initialization on the set X of service order sequences of the service satellites; X = [X1, X2, …, X i′ , …, X NP ; X i′ = [x1, x2, …, x i , …, x n ; where NP is the population size, that is, the number of service star service order sequences in the service star service order sequence set; X i′ is the i'-th service order sequence in the service star service order sequence set; x i is the number of the i-th target star; Step 2-2-1-2: Calculate the fitness value. The specific process is as follows: Assign each X i′ to step two two two, obtaining each X i′ corresponding service time T i′ , time node sequence S when the service star returns to the service station i′ and fitness value F i′ ; Compare the fitness values to obtain the optimal individual in this generation, update the optimal individual, and record it as the service order sequence, service time, time node sequence for the service satellites to return to the service station, and fitness value; Step 213: Perform a selection operation on the sequence X = [X1, X2, …, X i′ , …, X NP : Adopt the roulette wheel strategy; Step 2-2-1-4: Perform the crossover operation; Then perform the de-duplication operation: Replace the repeated genes from left to right with those that have not appeared; Step 2-2-1-5: Mutation operation: Randomly select two genes in the chromosome and exchange them; Step 2-2-1-6: Perform chaotic perturbation; Step 2217: Determine whether the iteration requirements are met. If so, end and obtain the optimal service order sequence X best , the corresponding optimal service time T best , the time node sequence S when the optimal service star returns to the service station best and the optimal fitness value F best ; if the iteration requirements are not met, execute Step 2212; In Step 2-2-1-6, the chaotic perturbation is performed. The specific process is as follows: First, map the X after the mutation operation in Step 2215 i′ to the decimal number z between 0 and 1 i′ , and the specific mapping method is as follows: X after the sequence undergoes the mutation operation of step 2215 i′ They are arranged in ascending order as n!, among which the decimal number mapped by the Set the chaotic perturbation coefficient k is the number of iterations; Set the perturbation amount Δz = rand × α × (-1) randi ×p; rand is a random number from 0 to 1, randi is 0 or 1, and p needs to be determined according to the length of n; At z i′ After adding the perturbation amount Δz, a new decimal z i′ ′; Map to sequence X according to the mapping method in step 2211 i′ ′, assign sequence X i′ ′ to step 222 to obtain the service time T i′ ′ corresponding to each X i′ ′, the time node sequence S i′ ′ when the corresponding service star returns to the service station, and the corresponding fitness value F i′ ′; If a smaller fitness value is obtained, then X i′ ' is replaced with X i′ , resulting in a new sequence X = [X1, X2, …, X i′ ', …, X NP ​ In the second step 222, the particle swarm optimization algorithm is used to optimize the service time T best ; The specific process is as follows: Step 2-2-2-1: Initialize the particle swarm parameters. The specific process is as follows: Randomly generate a particle swarm \(T = [T_1, T_2, \ldots, T i′ , \ldots, T NP \), where \(T i′ = [t_1, t_2, \ldots, t i , \ldots, t n \), \(NP\) is the number of particle swarms, satisfying the constraint \(t i \in [1, T max \), \(\sum(t i )\leq T max ;\) Randomly generate the particle velocity V; Step Two Two Two Two: Combine T1, T2, …, T i′ , …, T NP with the service order in Step Two Two One respectively and assign them to Step Two Two Three to calculate the time node sequence for the service stars to return to the service station and the corresponding fitness values; Record the time series corresponding to the minimum fitness value as the optimal service time; Step Two Two Two Three: Update the velocity V and randomly generate the particle population T using the particle swarm method NP , where the inertia weight is: where iter_max is the maximum number of iterations, w min represents the minimum inertia weight, w max represents the maximum inertia weight, and iter represents the current iteration number; Step Two Two Two Four, return F if the iteration count is met min , and the corresponding time series; otherwise, execute Step Two Two Two Two; In step 223, the branch and bound algorithm is used to optimize the time node for returning to the service station; the specific process is as follows: Step 2231: Initialize the parameters of the branch and bound algorithm; the specific process is as follows: Initialize the matrix ind0 = [0] and ind1 = [1], and their combined matrix Ind = [ind0; ind1]. Since the last digit of the node returning to the service station must be 1, only the first n - 1 digits need to be solved, where n is the number of target stars; Step 2232: Based on T i′ and the service order in Step 221, optimize and return the time nodes of the service station and the corresponding fitness values.

2. The method for planning an in-orbit refueling mission coplanar with the GEO orbit according to claim 1, wherein: The chaotic initialization service star service order sequence set X is carried out by using the mapping method from integer permutation to decimal; X = [X1, X2, …, X i′ , …, X NP ; X i′ = [x1, x2, …, x i , …, x n ; where NP is the population size, i.e., the number of service star service order sequences in the service star service order sequence set; X i′ is the i'-th service order sequence in the service star service order sequence set; x i is the number of the i-th target star; The specific method is as follows: Suppose there are n target stars, and there are n! permutation methods; Set z0 = 0.2 to obtain the initial value; Select the motion equation z i″+1 = 4z i″ (1 - z i″ ), i″ = 1, … n, … 3NP, map the decimal z i″ to the permutation X j , j = 1, 2, … n!; The specific mapping method is as follows: Arrange the \(n!\) permutations in ascending order, then \(z\) i represents the \(round(z\) i \(\times n!)\)-th permutation; where round is rounding; Obtain 3NP service sequence populations X 3NP , and the corresponding time series T 3NP Randomly select and execute Step Two Two Three to obtain the corresponding node S of the return service station 3NP and the corresponding fitness value F 3NP ; Take the first NP service order sequences X = [X1, X2, …, X i′ , …, X NP with smaller fitness values as the initial population.

3. A method for planning in-orbit refueling tasks for GEO orbit coplanarity according to claim 2, characterized in that: In step 2232, according to T i′ and the service order in step 221, the time node for optimizing the return to the service station and the corresponding fitness value are combined; the specific steps are as follows: Step 1: Obtain the number of rows and columns of Ind; if the number of columns is less than n - 1, execute Step 2, otherwise execute Step 3; Step 2: Fill 0 after each row in matrix ind0 and fill 1 after each row in matrix ind1. According to T i′ and the service order of Step 221, calculate the fitness value of each row in matrix ind1, compare the calculated fitness values of each row, and the row with the smallest fitness value among all rows is the partially optimal solution I in ind1 opt , and update Ind = [ind0; I opt , and execute Step 1; Step 3: Fill 1 after each row in matrix ind1, and calculate the fitness value of each row in matrix ind1 according to T i′ and the service order in Step 221. Compare the calculated fitness values of each row. The row with the minimum fitness value among all rows is the optimal complete solution in ind1; Step 4: Obtain the optimal complete solution and the corresponding fitness value.

Images