Satellite constellation fuel filling task planning method based on improved ant colony algorithm

Abstract

The invention discloses a satellite constellation fuel filling task planning method based on an improved ant colony algorithm, and belongs to the technical field of satellite task planning. The problems that an optimal solution is difficult to search and the search efficiency is low in an existing method are solved. The method specifically comprises the following steps: step 1, initializing parameters of a DQN model, and training the initial DQN model to obtain a trained DQN model; 2, initializing parameters of an ant colony optimization algorithm, establishing a fitness function of the ant colony optimization algorithm, and establishing an inter-satellite transition probability equation based on the trained DQN model; and step 3, searching based on the established fitness function and transition probability equation to obtain a task planning result. According to the invention, the Q value of the DQN model is integrated into the ACO decision making process, so that the ant not only considers the pheromone trajectory and the heuristic distance, but also considers the learned Q value in the node selection process, and the movement decision of the ant is guided through the Q value. The method can be applied to satellite constellation fuel filling task planning.

Description

Technical Field

[0001] The present invention belongs to the technical field of satellite mission planning, and particularly relates to a satellite constellation fueling mission planning method based on an improved ant colony algorithm. Background Art

[0002] With the continuous progress of space technology, large-scale low Earth orbit (LEO) satellite constellations have become an indispensable infrastructure in the space field. These large constellation systems increase the complexity of mission planning requirements, including tasks such as fuel replenishment, maintenance operations, and orbital debris removal. To complete these tasks, service satellites usually need to visit other members within the constellation in a specific order, thus triggering complex combinatorial optimization problems involving orbital dynamics.

[0003] Theoretically, the optimal sequence search problem in a large-scale constellation can be abstracted as a three-dimensional traveling salesman problem (3D TSP) for ΔV. In the space field, ΔV refers to the velocity change required to perform one or more maneuvers to change the orbit of a spacecraft. ΔV is not only a basic index for evaluating mission design and orbital transfer efficiency, but also of great significance for ensuring the economy and feasibility of space missions.

[0004] However, in large-scale LEO constellations, the number of satellites involved in a mission usually ranges from several hundred to several thousand, making it very difficult to find an optimal or near-optimal access sequence within polynomial time. Traditional greedy algorithms have a high computational complexity when solving such problems.

[0005] In recent years, metaheuristic algorithms have received extensive attention due to their effectiveness in solving TSP problems. Murakami et al. successfully applied metaheuristic algorithms to the preliminary research of active debris removal (ADR) tasks. Inspired by this, Missel applied the genetic algorithm (GA) to orbital maintenance tasks, especially those involving the four-satellite launcher (4S) project. Similarly, Liu used GA to solve multi-objective optimization problems in space mission planning. Swarm intelligence algorithms, such as ant colony optimization (ACO) and particle swarm optimization (PSO), have shown high efficiency in solving problems with complex search spaces, thus gaining wide recognition in the field of mission planning. In particular, ACO, which mimics the foraging behavior of ants, has been proven to be an effective method for optimizing path selection. Variants and hybrid algorithms of GA have also been valued for their innovative adaptive solutions and ability to balance exploration and exploitation. However, as the problem scale expands, the search space grows exponentially, making it increasingly difficult to find the optimal solution, which poses a greater challenge to metaheuristic algorithms in searching for the global optimal solution. In addition, the populations generated by these algorithms in large-scale problems usually lack effective evolutionary dynamics, reducing the search efficiency. Summary of the Invention

[0006] The object of the present invention is to solve the problem that with the exponential growth of the search space, it is difficult for existing methods to search for the optimal solution and the search efficiency is low, and a satellite constellation fuel filling mission planning method based on an improved ant colony algorithm is proposed.

[0007] The technical solution adopted by the present invention to solve the above technical problems is: a satellite constellation fuel filling mission planning method based on an improved ant colony algorithm, and the method specifically includes the following steps:

[0008] Step 1, initialize the parameters of the DQN model, and after training the initial DQN model, obtain the trained DQN model;

[0009] Step 2, initialize the parameters of the ant colony optimization algorithm, establish the fitness function of the ant colony optimization algorithm, and establish the transition probability equation between satellites based on the trained DQN model;

[0010] The parameters of the ant colony optimization algorithm specifically include the number of ants m, the importance of pheromone α, the importance of heuristic β, the importance of Q value ξ, and the pheromone evaporation rate ρ;

[0011] Step 3, perform search based on the established fitness function and transition probability equation to obtain the task planning result.

[0012] Further, the specific process of training the initial DQN model is:

[0013] Step 11, randomly select a satellite x i as the starting point, and then calculate the costs required to go from the starting point to each of the other satellites in the constellation respectively. Use all the calculated costs to form the distance vector of satellite x i , and use the distance vector of satellite x i as the current state s;

[0014] Step 12, convert the current state s into a tensor and then input it into the DQN model; the DQN model selects the next satellite x j , x j based on the current state s, which is the action a output by the DQN model. The value network in the DQN model predicts the Q value of selecting the action a in the current state s, and denote the predicted Q value as Q(s, a; θ), where θ represents the parameters of the main value network;

[0015] Step 13, calculate the reward r obtained by executing the action a:

[0016]

[0017] where, represents the orbit transfer from satellite x i to satellite xj Velocity change;

[0018] Step 14. Calculate the costs that satellites other than satellite x in the constellation need to pay respectively. Use all the calculated costs to form the distance vector of satellite x j to each of the other satellites in the constellation except satellite x i . Use the calculated distance vector of satellite x j as the new state s′, and convert the new state s′ into a tensor and then input it into the DQN model; j

[0019] Step 15. The DQN model selects the next action a′ based on the new state s′. The value network in the DQN model predicts the Q value of selecting action a′ in the new state s′;

[0020] Step 16. Update the Q value of selecting action a in the current state s:

[0021] Q′(s,a) = r + γmax a′ Q(s′,a′;θ - )

[0022] where Q′(s,a) is the updated Q value of selecting action a in the current state s, γ represents the discount factor, and Q(s′,a′;θ - ) is the Q value of selecting action a′ in the new state s′, and θ - represents the parameters of the value network;

[0023] Step 17. Calculate the loss function according to Q′(s,a) and Q(s,a;θ), and adjust the parameters of the DQN model through backpropagation of the loss function;

[0024] Step 18. Take the new state as the current state, and return to execute Step 12 until the set maximum number of iterations is reached and the training stops.

[0025] Furthermore, the velocity change i when satellite x j performs an orbit transfer to reach satellite x is:

[0026] Step 131. Denote the orbital parameters of satellite x i as D 1 = {a 1 , e 1 , i 1 , Ω 1 , ω 1 , υ 1}, and denote the orbital parameters of satellite x j as D 2 = {a 2 , e 2 ​, i 2 , Ω 2 , ω 2 , υ 2};

[0027] Among them, a 1 is the semi-major axis length of the elliptical orbit of satellite x i , e 1 is the eccentricity of the orbit of satellite x i , i 1 is the inclination angle of the orbit plane of satellite x i relative to the reference plane, Ω 1 is the right ascension angle of the ascending node of the orbit of satellite x i , υ 1 is the object position angle on the orbit of satellite x i , ω 1 is the argument of perigee of the orbit of satellite x i ; a 2 is the semi-major axis length of the elliptical orbit of satellite x j , e 2 is the eccentricity of the orbit of satellite x j , i 2 is the inclination angle of the orbit plane of satellite x j relative to the reference plane, Ω 2 is the right ascension angle of the ascending node of the orbit of satellite x j , υ 2 is the object position angle on the orbit of satellite x j , ω 2 is the argument of perigee of the orbit of satellite x j ;

[0028] Step 132. The instantaneous radius r of satellite x i on its own orbit is: 1

[0029]

[0030] The instantaneous radius r of satellite x j on its own orbit is: 2

[0031]

[0032] The instantaneous velocity v of satellite x i on its own orbit is: 1

[0033]

[0034] Among them, μ is the gravitational parameter;

[0035] Satellite x j ​​​The instantaneous velocity v in its own orbit 2 is:

[0036]

[0037] Step 133: Calculate the semi-major axis a of the Hohmann transfer orbit h :

[0038]

[0039] Step 134: Calculate the velocity at the intersection of the Hohmann transfer orbit:

[0040]

[0041] where v h1 is the velocity at the intersection of the Hohmann transfer orbit of satellite x i and v h2 is the velocity at the intersection of the Hohmann transfer orbit of satellite x j ;

[0042] Step 135: Calculate the required velocity change for satellite x i orbit and the required velocity change for satellite x j orbit:

[0043] ΔV 1 = |v 1 - v h1 |

[0044] ΔV 2 = |v 2 - v h2 |

[0045] where ΔV 1 is the required velocity change for satellite x i orbit, and ΔV 2 is the required velocity change for satellite x j orbit;

[0046] Step 136: Calculate the required velocity change ΔV for orbit plane adjustment plane :

[0047] Δi = |i 1 - i 2 |, ΔΩ = |Ω 1 - Ω 2 |

[0048] cos(α) = cos(Δi) + sin(i 1 )sin(i 2 )(1 - cos(ΔΩ))

[0049]

[0050] Step 137, calculate from satellite x i Perform an orbital transfer to reach satellite x j of the velocity change

[0051]

[0052] Further, in the said Step 17, the loss function is:

[0053] L(θ) = E[(r + γmax a′ Q(s′, a′; θ - ) - Q(s, a; θ)) 2

[0054] where L(θ) is the loss function, E[·] represents taking the expectation, θ - represents the parameters of the secondary value network, and θ represents the parameters of the primary value network.

[0055] Further, the working process of the said primary value network is:

[0056] Inside the primary value network, the input passes through the input layer, the first gated recurrent unit, the second gated recurrent unit, and the output layer in sequence, and the Q value is output through the output layer;

[0057] The first gated recurrent unit includes a reset gate and an update gate, and the working process of the first gated recurrent unit is:

[0058] Step 1, calculate the output r of the reset gate t :

[0059] r t = σ(W r · [h t-1 , x t )

[0060] where x t is the current input, h t-1 is the previous hidden state, W r represents the weight matrix of the reset gate, and σ is the sigmoid activation function;

[0061] Step 2, calculate the output z of the update gate t :

[0062] z t = σ(W z · [h t-1 , x t )

[0063] where W z is the weight matrix of the update gate;​

[0064] Step 3, calculate the candidate hidden state

[0065]

[0066] where ⊙ represents element-wise multiplication, tanh(·) represents the activation function, and W h represents the weight matrix of the candidate hidden state;

[0067] Step 4, update the hidden state h t :

[0068]

[0069] Furthermore, the working process of the main value network is the same as that of the slave value network, and the working process of the first gated recurrent unit is the same as that of the second gated recurrent unit.

[0070] Furthermore, the fitness function is:

[0071]

[0072] where x i represents the i-th satellite in the satellite constellation, x i+1 represents the (i + 1)-th satellite in the satellite constellation, n + 1 represents the total number of satellites in the satellite constellation, represents the transfer cost between continuously accessing the i-th satellite x i and the (i + 1)-th satellite x i+1 and minfit i is the fitness function.

[0073] Furthermore, the specific process of Step 3 is as follows:

[0074] Step 3-1, each ant randomly selects a satellite from the satellite constellation as the starting point;

[0075] Step 3-2, according to the established probability transfer equation, each ant calculates the transfer probability from the current satellite to each satellite in the set of satellites allowed to transfer to the next step;

[0076] Step 3-3, each ant selects the next satellite according to the maximum transfer probability corresponding to itself and updates the pheromone concentration on the path;

[0077] Step 3-4, each ant takes the satellite selected by itself in Step 3-3 as the current satellite and returns to execute Step 3-2 until the end of the entire search process, and respectively obtains the satellite access order sequence searched by each ant;

[0078] Step 35: Use the satellite access sequence corresponding to the minimum fitness function value as the final task planning result.

[0079] Furthermore, the transition probability equation is as follows:

[0080]

[0081] where is the transition probability of the k-th ant from satellite x i to satellite x j , represents the set of satellites that the k-th ant is allowed to transfer to next, τ i,j (t) is the pheromone concentration on the path between satellite xi and satellite xj before update, Q i,j (t) represents the Q value predicted by the trained DQN model for transferring from satellite x i to satellite x j , η i,j (t) is the heuristic value for transferring from satellite x i to satellite x j .

[0082] Even further, the pheromone concentration update method is as follows:

[0083]

[0084] where τ i,j (t + 1) is the pheromone concentration on the path between satellite x i and satellite x j after update, τ i,j (t) is the pheromone concentration on the path between satellite x i and satellite x j before update, represents the amount of pheromone left by the k-th ant on the path between satellite x i and satellite x j at step t.

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

[0086] The present invention proposes a data-driven method that combines ACO and DQN models, integrates the Q value of the DQN model into the ACO decision-making process. In this way, when ants select nodes, they not only consider pheromone trails and heuristic distances but also the learned Q value. By guiding the movement decisions of ants through the Q value, it is possible to enhance the exploration efficiency of the satellite constellation member access sequence while preventing premature convergence to suboptimal solutions and ultimately searching for the optimal solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0087] Figure 1It is the flowchart of a method for satellite constellation fuel filling mission planning based on an improved ant colony algorithm of the present invention;

[0088] Figure 2 It is the architecture diagram of the deep network for Q-value estimation;

[0089] Figure 3 It is the graph of the change in the ΔV value required for the transfer between satellites in Dataset I;

[0090] In the figure, Rendezvous point represents the rendezvous point;

[0091] Figure 4 It is the graph of the change in the ΔV value required for the transfer between satellites in Dataset II;

[0092] Figure 5 It is the graph of the change in the fitness and average fitness value of the method of the present invention on Dataset I;

[0093] In the figure, Iteration represents the number of iterations, Fitness represents the fitness value, global best represents the global optimal solution, and population average represents the average fitness value;

[0094] Figure 6 It is the graph of the change in the fitness and average fitness value of the method of the present invention on Dataset II;

[0095] Figure 7 It is the T-SNE visualization of the task planning result of the method of the present invention on Dataset I;

[0096] Figure 8 It is the T-SNE visualization of the task planning result of the method of the present invention on Dataset II;

[0097] Figure 9 It is log 500 The box plot of the iterative optimal solution at the scale. Detailed implementation manners

[0098] Detailed implementation manner 1: Combine Figure 1 To illustrate this implementation manner. A method for satellite constellation fuel filling mission planning based on an improved ant colony algorithm described in this implementation manner, the method specifically includes the following steps:

[0099] Step 1: Initialize the parameters of the DQN model, and after training the initial DQN model, obtain the trained DQN model;

[0100] Step 2: Initialize the parameters of the ant colony optimization algorithm, establish the fitness function of the ant colony optimization algorithm, and establish the transfer probability equation between satellites based on the trained DQN model;

[0101] The parameters of the ant colony optimization algorithm specifically include the number of ants m, the importance of pheromone α, the importance of heuristic β, the importance of Q value ξ, and the pheromone evaporation rate ρ;

[0102] Step 3: Search based on the established fitness function and transition probability equation to obtain the task planning result.

[0103] The purpose of the improved ant colony optimization algorithm (DDEACO) of the present invention is to find a solution close to or optimal within a short time. The DDEACO algorithm framework combines deep learning technology and traditional ACO to improve the quality and efficiency of the solution. Specifically, the present invention uses a deep Q-network (DQN) to preliminarily explore the solution space to effectively learn high-quality task planning solutions, and uses the ΔV transfer cost from the current point to the target point as a reward mechanism. In addition, a gated recurrent unit (GRU) is used to capture the time dependence in the optimal sequence, which is further used to determine the Q value of the subsequent rendezvous point. The pre-trained GRU network is integrated into the ACO algorithm, enabling the ACO algorithm to have rich heuristic path selection capabilities, thereby improving the efficiency of identifying the optimal solution.

[0104] Specific Embodiment 2: The difference between this embodiment and Specific Embodiment 1 is that the specific process of training the initial DQN model is as follows:

[0105] Step 11: Randomly select a satellite x i as the starting point, and then calculate the costs required for the starting point to reach each of the other satellites in the constellation. Use all the calculated costs to form the distance vector of satellite x i and use the distance vector of satellite x i as the current state s;

[0106] Step 12: Convert the current state s into a tensor and input it into the DQN model; the DQN model selects the next satellite x j , x j based on the current state s, where x i is the action a output by the DQN model. The value network in the DQN model predicts the Q value of selecting action a in the current state s, and denote the predicted Q value as Q(s, a; θ), where θ represents the parameters of the main value network;

[0107] Step 13: Calculate the reward r obtained by executing action a:

[0108]

[0109] where, represents the velocity change for orbital transfer from satellite x i to satellite x j ;

[0110] Step 14: Calculate the costs that each satellite in the constellation except satellite x needs to pay respectively, and use all the calculated costs to form the distance vector of satellite x. Take the distance vector of satellite x as the new state s′, and convert the new state s′ into a tensor and then input it into the DQN model; j to each of the other satellites in the constellation except satellite x i The calculated total cost forms the distance vector of satellite x j , and use the distance vector of satellite x j as the new state s′, and after converting the new state s′ into a tensor, input it into the DQN model;

[0111] Step 15: The DQN model selects the next action a′ based on the new state s′, and the value network in the DQN model predicts the Q value of selecting the action a′ in the new state s′;

[0112] Step 16: Update the Q value of selecting the action a in the current state s:

[0113] Q′(s,a) = r + γmax a′ Q(s′,a′;θ - )

[0114] where Q′(s,a) is the updated Q value of selecting the action a in the current state s, γ represents the discount factor, Q(s′,a′;θ - ) is the Q value of selecting the action a′ in the new state s′, and θ - represents the parameters of the value network;

[0115] Step 17: Calculate the loss function according to Q′(s,a) and Q(s,a;θ), and adjust the parameters of the DQN model through the backpropagation of the loss function;

[0116] Step 18: Take the new state as the current state, return to execute Step 12, and stop training until the set maximum number of iterations is reached.

[0117] Other steps and parameters are the same as those in the first specific implementation manner.

[0118] Specific implementation manner three: The difference between this implementation manner and the first or second specific implementation manner is that it is assumed that the satellite orbit has zero eccentricity and runs along a circular orbit. The speed change i for the satellite x j to perform an orbit transfer to reach satellite x is:

[0119] Step 131: Denote the orbital parameters of satellite x i as D 1 = {a 1 , e 1 , i 1 , Ω 1 , ω 1 , υ 1}, and denote the satellite x jThe orbital parameters are denoted as D 2 ={a 2 , e 2 , i 2 , Ω 2 , ω 2 , υ 2};

[0120] Among them, a 1 is the semi-major axis length of the elliptical orbit of satellite x i , e 1 is the eccentricity of the orbit of satellite x i (for a circular orbit, e 1 =0; for an elliptical orbit, 0 < e 1 <1), i 1 is the inclination angle of the orbital plane of satellite x i relative to the reference plane (usually relative to the Earth's equatorial plane or the ecliptic plane), Ω 1 is the right ascension angle of the ascending node of the orbit of satellite x i (i.e., the angle between the intersection point of the orbital plane and the reference plane (ascending node) and the reference direction (such as the vernal equinox point)), υ 1 is the position angle of the object on the orbit of satellite x i (indicating the position of the object on the orbit relative to the perigee), ω 1 is the argument of perigee of the orbit of satellite x i ; a 2 is the semi-major axis length of the elliptical orbit of satellite x j , e 2 is the eccentricity of the orbit of satellite x j (for a circular orbit, e 2 =0; for an elliptical orbit, 0 < e 2 <1), i 2 is the inclination angle of the orbital plane of satellite x j relative to the reference plane (usually relative to the Earth's equatorial plane or the ecliptic plane), Ω 2 is the right ascension angle of the ascending node of the orbit of satellite x j (i.e., the angle between the intersection point of the orbital plane and the reference plane (ascending node) and the reference direction (such as the vernal equinox point)), υ 2 is the position angle of the object on the orbit of satellite x j (indicating the position of the object on the orbit relative to the perigee), ω 2 is the argument of perigee of the orbit of satellite x j ;

[0121] Step 132. The instantaneous radius r of satellite x i on its own orbit 1is:

[0122]

[0123] Satellite x j The instantaneous radius r in its own orbit 2 is:

[0124]

[0125] Satellite x i The instantaneous velocity v in its own orbit 1 is:

[0126]

[0127] where μ is the gravitational parameter;

[0128] Satellite x j The instantaneous velocity v in its own orbit 2 is:

[0129]

[0130] Step 133. Calculate the semi-major axis a of the Hohmann transfer orbit h :

[0131]

[0132] Step 134. Calculate the velocity at the intersection of the Hohmann transfer orbit:

[0133]

[0134] where v h1 is the velocity at the intersection of the Hohmann transfer orbit of satellite x i and v h2 is the velocity at the intersection of the Hohmann transfer orbit of satellite x j ;

[0135] Step 135. Calculate the required velocity change of satellite x i orbit and the required velocity change of satellite x j orbit:

[0136] ΔV 1 = |v 1 - v h1 |

[0137] ΔV 2 = |v 2 - v h2 |

[0138] where ΔV 1 is the satellite x iVelocity change required for the orbit, ΔV 2 is the velocity change required for the orbit of satellite x j ;

[0139] Step 136: Calculate the velocity change ΔV required for orbit plane adjustment plane :

[0140] Δi = |i 1 - i 2 |, ΔΩ = |Ω 1 - Ω 2 |

[0141] cos(α) = cos(Δi) + sin(i 1 )sin(i 2 )(1 - cos(ΔΩ))

[0142]

[0143] Step 137: Calculate the velocity change for satellite x i to perform an orbit transfer to reach satellite x j ;

[0144]

[0145] Other steps and parameters are the same as those in the first or second specific implementation manner.

[0146] Specific implementation manner four: The difference between this implementation manner and one of the first to third specific implementation manners is that in step 17, the loss function is:

[0147] L(θ) = E[(r + γmax a′ Q(s′, a′; θ - ) - Q(s, a; θ)) 2

[0148] where L(θ) is the loss function, E[·] represents taking the expectation, θ - represents the parameters of the secondary value network, and θ represents the parameters of the primary value network.

[0149] Other steps and parameters are the same as those in one of the first to third specific implementation manners.

[0150] The parameters of the secondary value network are synchronized with the primary value network regularly. During the training process, ensure that the Q value of the action in a given state is recognized. The higher the Q value, the better the path selection.

[0151] Specific implementation manner five: Combine Figure 2 to illustrate this implementation manner. The difference between this implementation manner and one of the first to fourth specific implementation manners is that the working process of the primary value network is:​

[0152] Within the main value network, the input sequentially passes through the input layer, the first gated recurrent unit (GRU), the second gated recurrent unit, and the output layer, and the Q value is output through the output layer;

[0153] The first gated recurrent unit includes a reset gate and an update gate, and the working process of the first gated recurrent unit is as follows:

[0154] Step 1: Calculate the output r of the reset gate t :

[0155] r t = σ(W r · [h t-1 , x t )

[0156] where x t is the current input, h t-1 is the previous hidden state, W r represents the weight matrix of the reset gate, and σ is the sigmoid activation function;

[0157] Step 2: Calculate the output z of the update gate t :

[0158] z t = σ(W z · [h t-1 , x t )

[0159] where W z is the weight matrix of the update gate;

[0160] Step 3: Calculate the candidate hidden state

[0161]

[0162] where ⊙ represents element-wise multiplication, tanh(·) represents the activation function, and W h represents the weight matrix of the candidate hidden state;

[0163] Step 4: Update the hidden state h t :

[0164]

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

[0166] Embodiment Six: The difference between this embodiment and any one of Embodiments One to Five is that the working process of the main value network is the same as that of the secondary value network, and the working process of the first gated recurrent unit is the same as that of the second gated recurrent unit.

[0167] Other steps and parameters are the same as those in any one of Embodiments One to Five.

[0168] GRU implements a gating mechanism that can capture long-term dependencies, thus alleviating the problems of gradient vanishing and explosion.

[0169] Embodiment Seven: The difference between this embodiment and any one of Embodiments One to Six is that the fitness function is as follows:

[0170]

[0171] where x i represents the i-th satellite in the satellite constellation, x i+1 represents the (i + 1)-th satellite in the satellite constellation, n + 1 represents the total number of satellites in the satellite constellation, represents the transfer cost between continuously accessing the i-th satellite x i and the (i + 1)-th satellite x i+1 minfit i is the fitness function, and the smaller the fitness function value, the better.

[0172] Other steps and parameters are the same as those in any one of Embodiments One to Six.

[0173] Embodiment Eight: The difference between this embodiment and any one of Embodiments One to Seven is that the specific process of Step Three is as follows:

[0174] Step 3-1: Each ant randomly selects a satellite from the satellite constellation as the starting point;

[0175] Step 3-2: According to the established probability transfer equation, each ant calculates the transfer probability from the current satellite to each satellite in the set of satellites allowed to be transferred in the next step (removing the satellites that have been selected from the set of satellites allowed to be transferred);

[0176] Step 3-3: Each ant selects the next satellite according to the maximum transfer probability corresponding to itself and updates the pheromone concentration on the path;

[0177] Step 3-4: Each ant takes the satellite selected by itself in Step 3-3 as the current satellite and returns to execute Step 3-2 until the end of the entire search process, respectively obtaining the satellite access order sequence searched by each ant;

[0178] Step 35: Use the satellite access sequence corresponding to the minimum fitness function value as the final mission planning result.

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

[0180] Specific Embodiment 9: The difference between this embodiment and any one of the first to eighth specific embodiments is that the transition probability equation is:

[0181]

[0182] where is the transition probability of the k-th ant from satellite x i to satellite x j , represents the set of the next allowed transition satellites of the k-th ant, τ i,j (t) is the pheromone concentration before update on the path between satellite xi and satellite xj, Q i,j (t) represents the Q value predicted by the trained DQN model for the transition from satellite x i to satellite x j , η i,j (t) is the heuristic value for the transition from satellite x i to satellite x j , and η i,j (t) is usually the reciprocal of the distance or time delay between satellite x i and satellite x j .

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

[0184] The present invention comprehensively considers the influence of pheromone concentration (τ i,j (t)), heuristic information (η i,j (t)), and Q value (Q i,j (t)) on the transition probability. The coefficients α, β, and ξ respectively adjust the influence of pheromone, heuristic information, and Q value to ensure that all factors are considered in balance when selecting a path. The path search benefits from the pheromone-guided exploration of ACO and the prediction ability of DQN, promoting a more effective search for the optimal satellite access sequence.

[0185] Specific Embodiment 10: The difference between this embodiment and any one of the first to ninth specific embodiments is that the pheromone concentration update method is:

[0186]

[0187] where τ i,j (t + 1) is between satellite x i and satellite x jThe updated pheromone concentration τ on the path between i,j (t) is satellite x i and satellite x j The pheromone concentration before update on the path between represents the amount of pheromone left by the k-th ant at the t-th step on the path from satellite x i to satellite x j on the path between.

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

[0189] Experimental part

[0190] To verify the effectiveness of the algorithm proposed in the present invention, two datasets of different scales were generated using STK (Satellite Tool Kit) for numerical analysis, namely datasets containing 30 satellites and 200 satellites. Table 1 shows the detailed information of the datasets. The satellite orbits in all datasets have nearly zero eccentricity, indicating mainly circular orbits. The key orbital parameters affecting the calculation of ΔV required for orbital transfer include the semi-major axis (a), orbital inclination (i), and right ascension of the ascending node (Ω). These parameters define the spatial orientation and shape of the orbit and affect the energy required for transfer maneuvers.

[0191] Table 1 Details of the datasets of the present invention

[0192]

[0193]

[0194] The transfer cost matrix of the dataset was calculated through ΔV, and 3D visualization was used to assist understanding, Figure 3 and Figure 4 shows the ΔV required for transfer between satellites, with darker colors indicating higher ΔV values.

[0195] A variety of benchmark methods were selected for comparison, including Deep Q-Network (DQN), Asynchronous Advantage Actor-Critic (A3C), Monte Carlo Tree Search (MCTS), Ant Colony Optimization (ACO) and its variant Max-Min Ant System (MMAS), as well as Genetic Algorithm (GA), Particle Swarm Optimization (PSO), Quantum Particle Swarm Optimization (QPSO), Adaptive Firefly Algorithm (AFSA), Artificial Bee Colony (ABC), and Sparrow Search Algorithm (SSA).

[0196] In the DDEACO framework of the present invention, the maximum number of iterations of the ACO algorithm is set to 300, and the maximum number of iterations of DQN is set to 100000. The GRU unit is adjusted according to the size of the data set, with a learning rate of 0.01 and a discount factor of 0.99. The number of ants in the ACO algorithm is 25, the pheromone importance (α) is 2, the heuristic information factor (β) is 4, the Q value importance (ξ) is 1, the pheromone evaporation rate (ρ) is 0.3, and the exploration factor (Q) is 1. The population size of all algorithms is consistent with 30. All reinforcement learning models use the Adam optimizer. The algorithm is implemented using Python 3.10 and PyTorch, and the calculation process is performed on the Intel i9-12900k CPU and RTX 4090 GPU on the Ubuntu system.

[0197] 1. Results of DDEACO on Datasets I and II

[0198] The effectiveness of DDEACO is demonstrated through computational experiments on two datasets, focusing on the fundamental mission planning challenge of optimizing satellite constellation access sequences. The core evaluation metric is the total ΔV required for path transfer, which is used as a fitness function to identify the optimal sequence.

[0199] Figure 5 and Figure 6 The fitness and average fitness value changes of DDEACO in 300 iterations are shown, showing the performance dynamics of DDEACO. The fitness curve shows that DDEACO combines the fast convergence characteristics of the traditional ACO algorithm with the evolutionary ability of reinforcement learning, ensuring its effectiveness in handling complex optimization tasks with high-dimensional search spaces.

[0200] The best solution of DDEACO on datasets I and II is 80.17 and 84.00 respectively. For detailed analysis, the t-distributed stochastic neighbor embedding (T-SNE) method is used for visualization, as shown in Figure 7 and Figure 8 As shown in Figure 1, T-SNE reduces the data dimension to two dimensions while retaining the original distance, making it easier to understand the performance of the algorithm in the search space.

[0201] 2. Compare the performance of DDEACO and the benchmark method in terms of ΔV optimal solution and computational complexity. The specific results are shown in Table 2.

[0202] Table 2 Comparison of ΔV optimal solutions and computational complexity of different algorithms

[0203]

[0204] In Table 2, "Time complexity" represents time complexity, "Space complexity" represents space complexity, "Reinforcement Learning" represents reinforcement learning, "Metaheuristic" represents metaheuristic algorithm, "RandomSearch" represents random search, "Data-driven" represents data-driven, and "EnhancedMetaheuristic" represents enhanced metaheuristic algorithm. The comparative analysis results are as follows:

[0205] 1. DDEACO algorithm: The DDEACO algorithm performs significantly better than other benchmark methods on datasets of different scales. Especially on Dataset II, the optimal ΔV of DDEACO is 84.0 m / s, demonstrating the effectiveness of DDEACO.

[0206] 2. Reinforcement learning methods: The ΔV of DQN and A3C on Dataset II are 3994.79 m / s and 3753.04 m / s respectively, and their time and space complexities are both O(N + N 2 ), and their performance is significantly inferior to DDEACO.

[0207] 3. Random search method: The ΔV of MCTS on Dataset II is 10962.0 m / s, indicating its limited adaptability in large-scale constellation mission planning.

[0208] 4. Ant colony optimization-based methods: The performance difference between ACO and MMAS is small. The optimal ΔV on Dataset II are 108.54 m / s and 161.47 m / s respectively, and they perform better among other benchmark algorithms, but their capabilities are still inferior to DDEACO.

[0209] 5. Other metaheuristic methods:

[0210] GA and QPSO perform similarly on Dataset II, and their ΔV are 3470.82 m / s and 3582.46 m / s respectively. Although their time complexities are both O(N), their optimal solutions are significantly inferior to DDEACO. The performance of AFSA, PSO, ABC, and SSA drops significantly in a larger solution space and may not be suitable for large-scale swarm mission planning.

[0211] In addition, we conducted ten independent optimization iterations to statistically analyze the robustness of the algorithms. Due to the inherent variability of population initialization and network weight positions, the results of each iteration may be different. Figure 9 The box plot shown presents the distribution of the optimal solutions after each iteration. To enhance the visualization of the differences, the data is presented on a log 500 scale.

[0212] FromFigure 9 It can be clearly seen that DDEACO consistently shows excellent accuracy in multiple random initializations, and the box plot of DDEACO is more compact, proving the strong robustness of the algorithm. In contrast, the benchmark algorithm is inferior to DDEACO in both accuracy and robustness.

[0213] The above numerical examples of the present invention are only for illustrating in detail the calculation model and calculation process of the present invention, rather than limiting the implementation manners of the present invention. For those of ordinary skill in the art, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation manners here. Any obvious changes or modifications derived from the technical solutions of the present invention still fall within the protection scope of the present invention.

Claims

1. A satellite constellation fuel filling mission planning method based on improved ant colony algorithm, characterized in that: The method specifically comprises the following steps: Step 1: Initialize the parameters of the DQN model, and after training the initial DQN model, obtain the trained DQN model; Step 2: Initialize the parameters of the ant colony optimization algorithm, establish the fitness function of the ant colony optimization algorithm, and establish the transition probability equation between satellites based on the trained DQN model; The parameters of the ant colony optimization algorithm specifically include the number of ants m, pheromone importance α, heuristic importance β, Q value importance ξ and pheromone evaporation rate ρ; Step 3: Search based on the established fitness function and transition probability equation to obtain the task planning result.

2. According to claim 1, a satellite constellation fuel filling mission planning method based on improved ant colony algorithm is characterized in that: The specific process of training the initial DQN model is as follows: Step 1: Randomly select a satellite x i As the starting point, calculate the cost of each satellite in the constellation from the starting point, and use the calculated total cost to form the satellite x i The distance vector of satellite x i The distance vector is taken as the current state s; Step 1 and 2: Convert the current state s into a tensor and input it into the DQN model; the DQN model selects the next satellite x based on the current state s j , x j That is, action a is output by the DQN model. The value network in the DQN model predicts the Q value of selecting action a under the current state s. The predicted Q value is recorded as Q(s,a;θ), where θ represents the parameters of the main value network. Step 13: Calculate the reward r obtained by executing action a: r=-ΔV xi,xj Where, ΔV xi,xj Represents the satellite x i Perform orbit transfer to reach satellite x j Speed ​​changes; Step 14: Calculate the satellite x separately j Remove satellite x from the constellation i The cost of each satellite other than , and the total cost calculated to form satellite x j The distance vector of satellite x j The distance vector is taken as the new state s′, and the new state s′ is converted into a tensor and input into the DQN model; Step 15: The DQN model selects the next action a′ based on the new state s′. The value network in the DQN model predicts the Q value of selecting action a′ under the new state s′. Step 16: Update the Q value of action a in the current state s: Q′(s,a)=r+γmax a′ Q(s′,a′;θ - ) Where Q′(s,a) is the updated Q value of action a in the current state s, γ represents the discount factor, Q(s′,a′; θ - ) is the Q value of action a′ in the new state s′, θ - Represents the parameters of the sub-value network; Step 17: Calculate the loss function based on Q′(s,a) and Q(s,a;θ), and adjust the parameters of the DQN model through back propagation of the loss function; Step 18: Use the new state as the current state and return to step 12 to stop training until the maximum number of iterations is reached.

3. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 2 is characterized in that: The satellite x i Perform orbit transfer to reach satellite x j Speed ​​change for: Step 131: Satellite x i The orbital parameters of the satellite x are recorded as D1 = {a1, e1, i1, Ω1, ω1, υ1}. j The orbital parameters are recorded as D2 = {a2, e2, i2, Ω2, ω2, υ2}; Where a1 is the satellite x i The semi-major axis length of the elliptical orbit, e1 is the length of the satellite x i The eccentricity of the orbit, i1 is the satellite x i The inclination of the orbital plane relative to the reference plane, Ω1 is the satellite x i The right ascension angle of the orbital ascending node, υ1, is the right ascension angle of the orbital ascending node. i The position angle of the object in orbit, ω1 is the satellite x i The orbit perigee argument; a2 is the satellite x j The semi-major axis length of the elliptical orbit, e2 is the length of the satellite x j The eccentricity of the orbit, i2 is the satellite x j The inclination of the orbital plane relative to the reference plane, Ω2 is the satellite x j The right ascension angle of the orbital ascending node, υ2, is the right ascension angle of the orbital ascending node. j The position angle of the object in orbit, ω2 is the satellite x j The orbital perigee argument; Step 132, Satellite x i The instantaneous radius r1 on its own orbit is: Satellite x j The instantaneous radius r2 on its own orbit is: Satellite x i The instantaneous speed v1 on its own track is: Among them, μ is the gravitational parameter; Satellite x j The instantaneous speed v2 on its own track is: Step 133: Calculate the semi-major axis a of the Hohmann transfer orbit h : Step 1, 3, 4. Calculate the velocity of the intersection point of the Hohmann transfer orbit: Among them, v h1 It is on satellite x i The velocity of the Hohmann transfer orbit intersection, v h2 It is on satellite x j The velocity of the intersection point of the Hohmann transfer orbit; Step 135: Calculate satellite x i The required velocity change of the orbit and the satellite x j Required velocity change of track: ΔV1=|v1-v h1 | ΔV2=|v2-v h2 | Where ΔV1 is the satellite x i The velocity change required for the orbit, ΔV2, is the satellite x j the required speed change of the track; Step 136. Calculate the speed change ΔV required for track plane adjustment plane : Δi=|i1-i2|, ΔΩ=|Ω1-Ω2| cos(α)=cos(Δi)+sin(i1)sin(i2)(1-cos(ΔΩ)) Step 137. Calculate the x from the satellite i Perform orbit transfer to reach satellite x j Speed ​​change 4. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 3 is characterized in that: In step 17, the loss function is: L(θ)=E[(r+γmax a′ Q(s′,a′;θ - )-Q(s,a;θ)) 2 ] Among them, L(θ) is the loss function, E[·] represents the expectation, and θ - represents the parameters of the secondary value network, and θ represents the parameters of the primary value network.

5. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 4 is characterized in that: The working process of the main value network is as follows: In the main value network, the input passes through the input layer, the first gated recurrent unit, the second gated recurrent unit and the output layer in sequence, and the Q value is output through the output layer; The first gated recurrent unit includes a reset gate and an update gate. The working process of the first gated recurrent unit is as follows: Step 1: Calculate the output r of the reset gate t : r t =σ(W r ·[h t-1 ,x t ]) Among them, x t is the current input, h t-1 is the previous hidden state, W r represents the weight matrix of the reset gate, σ is the sigmoid activation function; Step 2: Calculate the output z of the update gate t : z t =σ(W z ·[h t-1 ,x t ]) Among them, W z is the weight matrix of the update gate; Step 3: Calculate candidate hidden states Among them, ⊙ represents element-by-element multiplication, tanh(·) represents the activation function, and W h A weight matrix representing candidate hidden states; Step 4: Update the hidden state h t :

6. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 5 is characterized in that: The working process of the main value network is the same as the working process of the slave value network, and the working process of the first gated loop unit is the same as the working process of the second gated loop unit.

7. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 6 is characterized in that: The fitness function is: Among them, x i represents the i-th satellite in the satellite constellation, x i+1 represents the i+1th satellite in the satellite constellation, n+1 represents the total number of satellites in the satellite constellation, Indicates continuous visit to the i-th satellite x i and the i+1th satellite x i+1 The transfer cost between i is the fitness function.

8. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 7 is characterized in that: The specific process of step three is: Step 31. Each ant randomly selects a satellite from the satellite constellation as the starting point; Step 32: According to the established probability transfer equation, each ant calculates the transfer probability from the current satellite to each satellite in the next allowed transfer satellite set; Step 3. Each ant selects the next satellite according to its corresponding maximum transfer probability and updates the pheromone concentration on the path; Step 34: Each ant uses the satellite selected by itself in step 33 as the current satellite, and returns to execute step 32 until the entire search process is completed, and the satellite access sequence searched by each ant is obtained; Step 35: The satellite access sequence corresponding to the minimum fitness function value is taken as the final mission planning result.

9. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 8, characterized in that: The transition probability equation is: in, is the kth ant from satellite x i Transfer to Satellite x j The transition probability, represents the next set of satellites that the k-th ant is allowed to transfer to, τ i,j (t) is the pheromone concentration before updating on the path between satellite xi and satellite xj, Q i,j (t) represents the predicted value from satellite x by the trained DQN model i Transfer to Satellite x j Q value, η i,j (t) is the distance from satellite x i Transfer to Satellite x j The inspiration value of .

10. The satellite constellation fuel filling mission planning method based on improved ant colony algorithm according to claim 9, characterized in that: The pheromone concentration updating method is: Among them, τ i,j (t+1) is the updated pheromone concentration on the path between satellite xi and satellite xj, τ i,j (t) is the distance between satellite xi and satellite x j The pheromone concentration before the update on the path between them, Indicates that at step t, the kth ant is on satellite x i To Satellite x j The amount of pheromone left on the path between them.

Images