Satellite constellation fueling task planning method based on improved ant colony algorithm

Abstract

The application 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 application solves the problems that the existing method is difficult to search for an optimal solution and has low search efficiency.The application specifically comprises the following steps: step one, initializing parameters of a DQN model, training the initial DQN model, and obtaining a trained DQN model;step two, initializing parameters of an ant colony optimization algorithm, establishing a fitness function of the ant colony optimization algorithm, and establishing a transfer probability equation among satellites based on the trained DQN model;step three, searching based on the established fitness function and the transfer probability equation, and obtaining a task planning result.The application integrates Q values of the DQN model into an ACO decision process, so that the ants consider not only pheromone trails and heuristic distances but also learned Q values in a node selection process, and the movement decision of the ants is guided by the Q values.The application method can be applied to satellite constellation fuel filling task planning.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of satellite mission planning, and particularly relates to a satellite constellation fuel refueling mission planning method based on an improved ant colony algorithm. BACKGROUND

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

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

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

[0005] In recent years, metaheuristic algorithms have gained widespread attention due to their effectiveness in solving TSP problems. Murakami et al. successfully applied metaheuristic algorithms to the preliminary study of active debris removal (ADR) tasks. Inspired by this, Missel applied genetic algorithms (GA) to orbital maintenance tasks, particularly 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 performance in solving complex search space problems, thus gaining widespread recognition in the field of mission planning. In particular, ACO, which simulates 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 size expands, the search space grows exponentially, making it increasingly difficult to find the optimal solution, which poses greater challenges to metaheuristic algorithms in searching for global optimal solutions. In addition, the population generated by these algorithms in large-scale problems often lacks effective evolutionary dynamics, reducing search efficiency. SUMMARY

[0006] The application aims to solve the problem that the existing method is difficult to search for an optimal solution and has low search efficiency with the exponential growth of the search space, and proposes a satellite constellation fuel refueling task planning method based on an improved ant colony algorithm.

[0007] The technical scheme adopted by the application to solve the above technical problem is: a satellite constellation fuel refueling task planning method based on an improved ant colony algorithm, which specifically comprises the following steps:

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

[0009] Step two, 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 alpha, the importance of heuristic beta, the importance of Q value xi, and the evaporation rate of pheromone rho;

[0011] Step three, 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 one, randomly select a satellite x i as the starting point, and then calculate the cost of the starting point to each satellite in the constellation, use the calculated total cost 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 two, convert the current state s into a tensor and input it into the DQN model; the DQN model selects the next satellite x j based on the current state s, x j is the action a output by the DQN model, and the value network in the DQN model predicts the Q value of selecting action a under the current state s, and the predicted Q value is denoted as Q(s,a; theta), theta represents the parameters of the main value network;

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

[0016]

[0017] Wherein, represents the orbit transfer from satellite x i to satellite xj The speed change;

[0018] Step 1, Part 4: Calculate satellite x separately. j In the constellation, except for satellite X i The costs incurred by each of the other satellites, combined with the calculated total costs, form satellite x. j The distance vector of satellite x j The distance vector is used as the new state s′, and the new state s′ is converted into a tensor and then input into the DQN model;

[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' under the new state s'.

[0020] Step 16: Update the Q value of the selected 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 after choosing action a in the current state s, γ represents the discount factor, and Q(s′,a′; θ) - ) is the Q-value for choosing action a' in the new state s′, θ - Parameters representing the sub-value network;

[0023] Step 17: Calculate the loss function based on Q′(s,a) and Q(s,a;θ), and adjust the parameters of the DQN model through backpropagation using the loss function;

[0024] Step 18: Set the new state as the current state, return to Step 12, and continue training until the maximum number of iterations is reached.

[0025] Furthermore, the statement from satellite x i Perform orbital transfer to satellite x j speed change for:

[0026] Step 131: Place satellite x i The orbital parameters are denoted as D1={a1,e1,i1,Ω1,ω1,υ1}, and the satellite x j The orbital parameters are denoted as D2={a2,e2,i2,Ω2,ω2,υ2};

[0027] Where a1 is satellite x i The semi-major axis length of the elliptical orbit, e1 is the satellite x i The eccentricity of the orbit, i1 is the satellite x iinclination of the orbital plane with respect to the reference plane, Ω1 is the satellite x i right ascension of the ascending node of the orbit, υ1 is the satellite x i argument of the object's position on the orbit, ω1 is the satellite x i argument of the perigee of the orbit; a2 is the satellite x j length of the semi-major axis of the elliptical orbit, e2 is the satellite x j eccentricity of the orbit, i2 is the satellite x j inclination of the orbital plane with respect to the reference plane, Ω2 is the satellite x j right ascension of the ascending node of the orbit, υ2 is the satellite x j argument of the object's position on the orbit, ω2 is the satellite x j argument of the perigee of the orbit;

[0028] Step 1132, the satellite x i the instantaneous radius r1 on its own orbit is:

[0029]

[0030] the satellite x j the instantaneous radius r2 on its own orbit is:

[0031]

[0032] the satellite x i the instantaneous velocity v1 on its own orbit is:

[0033]

[0034] where μ is the gravitational parameter;

[0035] the satellite x j the instantaneous velocity v2 on its own orbit is:

[0036]

[0037] Step 1133, calculate the semi-major axis a of the Hohmann transfer orbit h :

[0038]

[0039] Step 1134, 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 the satellite x i is the velocity at the intersection of the Hohmann transfer orbit of the satellite x h2 is the velocity at the intersection of the Hohmann transfer orbit of the satellite x jthe speed of satellite x

[0042] Step 133, calculating the speed change ΔV1 required for satellite x i to transfer to the orbit of satellite x j Step 134, calculating the speed change ΔV2 required for satellite x

[0043] ΔV1 = |v1-v h1 |

[0044] ΔV2 = |v2-v h2 |

[0045] wherein ΔV1 is the speed change required for satellite x i to transfer to the orbit of satellite x j Step 135, calculating the speed change ΔV required for orbit plane adjustment plane :

[0046] ΔV = |v1-v2|

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

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

[0049]

[0050] Step 137, calculating the speed change ΔV required for satellite x i to transfer to the orbit of satellite x j :

[0051]

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

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

[0054] wherein L(θ) is the loss function, E[·] represents 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 primary value network is:

[0056] In the primary value network, the input sequentially passes through the input layer, the first gated recurrent unit, the second gated recurrent unit and the output layer, and the Q value is output through the output layer.

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

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

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

[0060] Wherein, 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 sigma is a sigmoid activation function.

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

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

[0063] Wherein, W z is the weight matrix of the update gate.

[0064] Step 3, calculating the candidate hidden state

[0065]

[0066] Wherein, represents element-by-element multiplication, tanh(·) represents an activation function, and W h represents the weight matrix of the candidate hidden state.

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

[0068]

[0069] Further, 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] Further, the fitness function is as follows:

[0071]

[0072] Wherein, x i represents the i-th satellite in the satellite constellation, and x i+1Xi+1 represents the i+1th satellite in the satellite constellation, and n+1 represents the total number of satellites in the satellite constellation, Xi represents the ith satellite in the satellite constellation, i Xi+1 represents the i+1th satellite in the satellite constellation, and n+1 represents the total number of satellites in the satellite constellation, i+1 minfit represents the transfer cost between the ith satellite x i and the i+1th satellite x i , minfit is the fitness function.

[0073] Further, the specific process of step three is as follows:

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

[0075] Step three, according to the established probability transfer equation, each ant respectively calculates the transfer probability of the current satellite to each satellite in the allowed transfer satellite set in the next step;

[0076] Step three, 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 three, each ant takes the satellite selected in step three as the current satellite, and returns to step three to perform step three, until the end of the whole search process, and respectively obtains the satellite access order sequence searched by each ant;

[0078] Step three, the satellite access order sequence corresponding to the minimum fitness function value is taken as the final task planning result.

[0079] Further, the transfer probability equation is:

[0080]

[0081] wherein, is the transfer probability of the kth ant from satellite x i to satellite x j , S i,j (t) represents the next allowed transfer satellite set of the kth ant, τ i,j (t) is the pheromone concentration on the path between satellite xi and satellite xj before updating, Q i (t) represents the Q value predicted by the trained DQN model from satellite x j to satellite x i,j , η i (t) is the heuristic value from satellite x j to satellite x i,j .

[0082] Further, the pheromone concentration updating method is:

[0083]

[0084] Where, τ i,j (t+1) is satellite x i With satellite x j The updated pheromone concentration, τ, along the path between them i,j (t) is satellite x i With satellite x j The pheromone concentration before the update along the path between them This indicates that at step t, the k-th ant is on satellite x. i To satellite x j The amount of pheromone left on the path between them.

[0085] The beneficial effects of this invention are:

[0086] This invention proposes a data-driven method that combines ACO with DQN models. By integrating the Q-value of the DQN model into the ACO decision-making process, ants not only consider pheromone trajectories and heuristic distances during node selection, but also the learned Q-value. By using the Q-value to guide the ants' movement decisions, this method can enhance the efficiency of exploring satellite constellation member access sequences, prevent premature convergence to suboptimal solutions, and ultimately find the optimal solution. Attached Figure Description

[0087] Figure 1 This is a flowchart of a satellite constellation refueling mission planning method based on an improved ant colony algorithm, according to the present invention.

[0088] Figure 2 This is a diagram of a deep network architecture used for Q-value estimation;

[0089] Figure 3 This is a graph showing the change in the ΔV values ​​required for transfer between satellites in dataset I;

[0090] In the diagram, the Rendezvous point represents the meeting point;

[0091] Figure 4 This is a graph showing the change in ΔV values ​​required for transfer between satellites in Dataset II;

[0092] Figure 5 This is a graph showing the changes in fitness and average fitness values ​​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 This is a graph showing the changes in fitness and average fitness values ​​of the method of the present invention on dataset II;

[0095] Figure 7 This is a visualization of the task planning results (T-SNE) of the method of the present invention on dataset I;

[0096] Figure 8 This is a visualization of the task planning results of the method of the present invention on dataset II using T-SNE.

[0097] Figure 9 It is log 500 Box plot of the iterative optimal solution at different scales. Detailed Implementation

[0098] Specific implementation method one: Combining Figure 1 This embodiment describes a satellite constellation refueling mission planning method based on an improved ant colony algorithm. The method specifically includes the following steps:

[0099] Step 1: Initialize the parameters of the DQN model and train the initial DQN model to 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 transition 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, pheromone importance α, heuristic importance β, Q-value importance ξ, and pheromone evaporation rate ρ.

[0102] Step 3: Based on the established fitness function and transition probability equation, perform a search to obtain the task planning results.

[0103] The improved Ant Colony Optimization (DDEACO) algorithm of this invention aims to find near-optimal or optimal solutions in a short time. The DDEACO algorithm framework combines deep learning techniques with traditional ACO to improve the quality and efficiency of solutions. Specifically, this invention employs a Deep Q-Network (DQN) to initially explore the solution space to effectively learn high-quality task planning schemes, and uses the ΔV transition cost from the current point to the target point as a reward mechanism. Furthermore, a Gated Recurrent Unit (GRU) is used to capture the temporal dependencies in the optimal sequence, which is further used to determine the Q-value of subsequent rendezvous points. Integrating a pre-trained GRU network into the ACO algorithm endows it with rich heuristic path selection capabilities, thereby improving the efficiency of identifying the optimal solution.

[0104] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that the specific process of training the initial DQN model is as follows:

[0105] Step one, randomly select a satellite x i As a starting point, the cost that each satellite in the constellation needs to pay to satellite x i is calculated respectively, and the distance vector of satellite x i is composed of all the calculated costs, and the distance vector of satellite x j is taken as the current state s.

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

[0107] Step three, the reward r obtained by executing action a is calculated:

[0108]

[0109] Wherein, represents the change of velocity of satellite x j when performing orbit transfer to satellite x j ;

[0110] Step four, the cost that each satellite in the constellation needs to pay to satellite x i is calculated respectively, and the distance vector of satellite x j is composed of all the calculated costs, and the distance vector of satellite x j is taken as the new state s', and the new state s' is converted into a tensor and input into the DQN model.

[0111] Step five, 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 action a' under the new state s'.

[0112] Step six, update the Q value of selecting action a under the current state s:

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

[0114] Wherein, Q'(s, a) is the updated Q value of selecting action a under the current state s, γ represents the discount factor, Q(s', a'; θ - ) is the Q value of selecting action a' under the new state s', and θ - represents the parameters of the secondary value network.

[0115] Step seventeen, calculate the loss function according to Q'(s, a) and Q(s, a; θ), and adjust the parameters of the DQN model by back propagation of the loss function;

[0116] Step eighteen, take the new state as the current state, and return to step twelve until the maximum number of iterations is reached.

[0117] The other steps and parameters are the same as in the first embodiment.

[0118] The third embodiment is different from the first or second embodiment in that the satellite orbit is assumed to have zero eccentricity and run along a circular orbit. The orbit parameters of the satellite x i The velocity change j for orbit transfer to the satellite x is:

[0119] Step thirteen, record the orbit parameters of the satellite x i as D1 = {a1, e1, i1, Ω1, ω1, υ1}, and the orbit parameters of the satellite x j as D2 = {a2, e2, i2, Ω2, ω2, υ2};

[0120] where a1 is the semi-major axis length of the elliptical orbit of the satellite x i , e1 is the eccentricity of the orbit of the satellite x i (0 < e1 < 1 for an elliptical orbit), i1 is the inclination of the orbit plane of the satellite x i relative to the reference plane (usually relative to the Earth's equatorial plane or ecliptic plane), Ω1 is the right ascension of the ascending node of the orbit of the satellite x i (i.e., the angle between the intersection of the orbit plane and the reference plane (ascending node) and the reference direction (e.g., the vernal equinox)), and υ1 is the argument of latitude of an object on the orbit of the satellite x i (indicating the position of the object on the orbit relative to the perigee), and ω1 is the argument of perigee of the satellite x i ; a2 is the semi-major axis length of the elliptical orbit of the satellite x j , e2 is the eccentricity of the orbit of the satellite x j (0 < e2 < 1 for an elliptical orbit), i2 is the inclination of the orbit plane of the satellite x j relative to the reference plane (usually relative to the Earth's equatorial plane or ecliptic plane), Ω2 is the right ascension of the ascending node of the orbit of the satellite x j (i.e., the angle between the intersection of the orbit plane and the reference plane (ascending node) and the reference direction (e.g., the vernal equinox)), and υ2 is the argument of latitude of an object on the orbit of the satellite x jThe object position angle (denoting the position of the object on the orbit relative to the perigee), ω2 is the satellite x j The orbit perigee argument (Argument of Perigee);

[0121] Step 132, the satellite x i The instantaneous radius r1 on the orbit of the satellite x

[0122]

[0123] The satellite x j The instantaneous radius r2 on the orbit of the satellite x

[0124]

[0125] The satellite x i The instantaneous velocity v1 on the orbit of the satellite x

[0126]

[0127] Wherein, μ is the gravitational parameter;

[0128] The satellite x j The instantaneous velocity v2 on the orbit of the satellite x

[0129]

[0130] Step 133, the semi-major axis a of the Hohmann transfer orbit is calculated h :

[0131]

[0132] Step 134, the velocity of the Hohmann transfer orbit intersection point is calculated

[0133]

[0134] Wherein, v h1 is the velocity of the Hohmann transfer orbit intersection point of the satellite x i , v h2 is the velocity of the Hohmann transfer orbit intersection point of the satellite x j ;

[0135] Step 135, the velocity change required by the satellite x i orbit and the velocity change required by the satellite x j orbit are calculated

[0136] ΔV1 = |v1-v h1 |

[0137] ΔV2 = |v2-v h2 |

[0138] where ΔV1 is the velocity change required for satellite x i to reach the orbit of satellite x j ;

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

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

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

[0142]

[0143] Step 137, calculate the velocity change ΔV required for satellite x i to reach the orbit of satellite x j from satellite x

[0144]

[0145] The other steps and parameters are the same as those in embodiment one or two.

[0146] Embodiment four: the difference between this embodiment and one of the embodiments one to three 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 the expectation, θ - represents the parameters of the secondary value network, and θ represents the parameters of the primary value network.

[0149] The other steps and parameters are the same as those in one of the embodiments one to three.

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

[0151] Embodiment five: in combination Figure 2 with the above embodiments, this embodiment is described. The difference between this embodiment and one of the embodiments one to four is that, the working process of the primary value network is:

[0152] In 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 outputs the Q value 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] Wherein, x t is the current input, h t-1 is the previous hidden state, W r is the weight matrix of the reset gate, and sigma 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] Wherein, W z is the weight matrix of the update gate.

[0160] Step 3, calculate the candidate hidden state

[0161]

[0162] Wherein, represents element-by-element multiplication, tanh(·) represents an 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] The other steps and parameters are the same as one of the first to fourth embodiments.

[0166] Embodiment six: The difference between this embodiment and one of the first to fifth embodiments is that 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.

[0167] The other steps and parameters are the same as one of embodiments one to five.

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

[0169] Embodiment seven: different from one of embodiments one to six, the fitness function is:

[0170]

[0171] Wherein, x i represents the i-th satellite in the satellite constellation, x i+1 represents the i+1-th satellite in the satellite constellation, and n+1 represents the total number of satellites in the satellite constellation, represents the transfer cost between 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 is, the better.

[0172] The other steps and parameters are the same as one of embodiments one to six.

[0173] Embodiment eight: different from one of embodiments one to seven, the specific process of step three is:

[0174] Step three one, each ant randomly selects a satellite from the satellite constellation as the starting point;

[0175] Step three two, each ant calculates the transfer probability of the current satellite to each satellite in the set of allowed transfer satellites (remove the satellite that has been selected from the set of allowed transfer satellites) according to the established probability transfer equation;

[0176] Step three three, 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 three four, each ant takes the satellite selected in step three three as the current satellite, and returns to step three two until the end of the search process, and respectively obtains the satellite access order sequence searched by each ant;

[0178] Step three five, the satellite access order sequence corresponding to the smallest fitness function value is taken as the final task planning result.

[0179] The other steps and parameters are the same as one of embodiments one to seven.

[0180] Specific implementation nine: different from one of the specific implementations one to eight, the transfer probability equation is:

[0181]

[0182] wherein, is the transfer probability of the kth ant from satellite x i to satellite x j , denotes the next allowed transfer satellite set of the kth ant, τ i,j (t) is the pheromone concentration on the path between satellite x i,j i and satellite x i j before updating, Q j (t) denotes the Q value predicted by the trained DQN model from satellite x i,j to satellite x i , j η i,j (t) is the heuristic value from satellite x i to satellite x j ,

[0183] Other steps and parameters are the same as one of the specific implementations one to eight.

[0184] The present application comprehensively considers the influence of pheromone concentration (τ i,j (t)), heuristic information (η i,j (t)) and Q value (Q i,j (t)) on transfer probability. Coefficients α, β and ξ adjust the influence of pheromone, heuristic information and Q value respectively, ensuring balanced consideration of all factors when selecting paths. Path search benefits from the pheromone-guided exploration of ACO and the prediction ability of DQN, promoting more effective search for the optimal satellite access sequence.

[0185] Specific implementation ten: different from one of the specific implementations one to nine, the pheromone concentration updating method is:

[0186]

[0187] wherein, τ i,j (t+1) is the updated pheromone concentration on the path between satellite x i and satellite x j , τ i,j (t) is the pheromone concentration on the path between satellite x i and satellite x j before updating, denotes the kth ant at step t in satellite x ithe amount of pheromone left on the path between satellite x j

[0188] The other steps and parameters are the same as one of the first nine embodiments.

[0189] Experimental Section

[0190] To verify the effectiveness of the proposed algorithm, two different size datasets were generated using STK (Satellite Tool Kit) for numerical analysis, containing 30 and 200 satellites, respectively. Table 1 shows the detailed information of the datasets. The orbits of all satellites in the datasets have near zero eccentricity, representing mainly circular orbits. The key orbital parameters that affect the required ΔV for inter-orbit transfers include 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, affecting the energy required for transfer maneuvers.

[0191] Table 1 Detailed information of the datasets of the present application

[0192]

[0193]

[0194] The transfer cost matrix of the datasets was obtained by ΔV calculation, and three-dimensional visualization was used to help understand, Figure 3 and Figure 4 The required ΔV for inter-satellite transfers is shown, with darker colors representing 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 application, the maximum number of iterations for the ACO algorithm is set to 300, and the maximum number of iterations for the DQN is set to 100,000. 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. In the ACO algorithm, the number of ants is 25, the importance of pheromone (a) is 2, the heuristic information factor (b) is 4, the Q value importance (z) is 1, the pheromone evaporation rate (p) is 0.3, and the exploration factor (Q) is 1. The population size of all algorithms is consistent at 30. All reinforcement learning models use the Adam optimizer. The algorithms are implemented using Python 3.10 and PyTorch, and the calculation process is performed on an Intel i9-12900k CPU and an RTX 4090 GPU on an Ubuntu system.

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

[0198] The present application demonstrates the effectiveness of DDEACO through computational experiments on two data sets, focusing on the basic mission planning challenge of optimizing satellite constellation access sequences. The core evaluation metric is the total AV required for path transfers, which is used as the 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 traditional ACO algorithms with the evolutionary capabilities of reinforcement learning, ensuring effectiveness in handling complex optimization tasks with high-dimensional search spaces.

[0200] The optimal solutions of DDEACO on data sets I and II are 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 T-SNE reduces the data to two dimensions while preserving the original distances, making it easier to understand the algorithm's performance in the search space.

[0201] 2. Comparison of DDEACO and benchmark methods in terms of AV optimal solution and computational complexity, as shown in Table 2.

[0202] Table 2 Comparison of AV optimal solution 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, Random Search represents random search, Data-driven represents data-driven, Enhanced Metaheuristic represents enhanced metaheuristic algorithm, and the comparative analysis results are as follows:

[0205] 1. DDEACO algorithm: The performance of the DDEACO algorithm on different sizes of data sets is significantly better than other benchmark methods. Especially on data set II, the optimal AV of DDEACO is 84.0 m / s, which proves the effectiveness of DDEACO.

[0206] 2. Reinforcement learning method: The AV of DQN and A3C on data set II is 3994.79 m / s and 3753.04 m / s, respectively, and the time and space complexity is O(N+N 2 ), which is obviously not as good as DDEACO.

[0207] 3. Random search method: The AV of MCTS on data set II is 10962.0 m / s, indicating that its adaptability in large-scale constellation task planning is limited.

[0208] 4. Ant colony optimization-based method: ACO and MMAS have small performance differences, with optimal AVs of 108.54 m / s and 161.47 m / s on data set II, respectively. They perform better than other benchmark algorithms, but their ability is still not as good as DDEACO.

[0209] 5. Other metaheuristic methods:

[0210] GA and QPSO have similar performance on data set II, with AVs of 3470.82 m / s and 3582.46 m / s, respectively. Although the time complexity of both is O(N), their optimal solutions are significantly worse than DDEACO. AFSA, PSO, ABC, and SSA have significantly decreased performance in larger solution spaces, which may not be suitable for large-scale group task planning.

[0211] In addition, we performed ten independent optimization iterations to calculate the robustness of the algorithm. Due to the inherent variability of population initialization and network weight positions, the results of each iteration may differ. Figure 9 The box plots shown show the distribution of optimal solutions after each iteration. To enhance the visualization of differences, data is shown in log 500 scale.

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

[0213] The above examples of the present application are only to illustrate the calculation model and calculation process of the present application, and are not limited to the embodiments of the present application. Based on the above description, other different forms of changes or variations can be made by those of ordinary skill in the art, and it is impossible to exhaust all the embodiments here. Any obvious changes or variations derived from the technical solutions of the present application are still within the protection scope of the present application.

Claims

1. A satellite constellation fuel refueling task planning method based on an improved ant colony algorithm, characterized in that, The method specifically comprises the following steps: Step one, initializing the parameters of the DQN model, and training the initial DQN model to obtain a trained DQN model; The specific process of training the initial DQN model is as follows: Step 11: Randomly select a satellite Starting from this point, the costs required to reach each of the other satellites in the constellation are calculated separately. All calculated costs are then used to assemble the satellite constellation. The distance vector of the satellite The distance vector is used as the current state. ; Step one: convert current state into tensor Input DQN model after conversion; DQN model based on current state Select next satellite , That is the action output by DQN model The value network in DQN model predicts the Q value of selecting action Under current state The predicted Q value is denoted as , Denote parameters of main value network Step one three, calculating to perform action Obtained reward : wherein, representing a satellite orbit transfer to a satellite velocity change; Step one four, the cost paid by each satellite in the constellation except the satellite is calculated respectively, and the distance vector of the satellite is composed of all the calculated costs, the distance vector of the satellite is taken as the new state , and the new state is converted into a tensor and input into the DQN model ; Step one five, the DQN model based on the new state selecting the next action , the value network in the DQN model predicts the Q value of selecting the action under the new state ; Step one six, update current state Next action to choose Q-value of wherein, is the current state the selected action the updated Q-value, denotes a discount factor, is the new state the selected action the Q-value, denotes the parameters of the value network; Step seventeen, according to and computing a loss function, adjusting parameters of the DQN model by backpropagation through the loss function; Step one eight, taking the new state as the current state, returning to step two, and stopping training until the maximum number of iterations is reached; Step two, initializing the parameters of the ant colony optimization algorithm, establishing the fitness function of the ant colony optimization algorithm, and establishing the transfer probability equation between satellites based on the trained DQN model; The fitness function is: wherein, denotes the i-th satellite of a satellite constellation, denotes the i-th satellite of a satellite constellation, denotes the i-th satellite of a satellite constellation, denotes the i-th satellite of a satellite constellation, denotes the total number of satellites of a satellite constellation, denotes the transfer cost between the i-th satellite and the i-th satellite and the i-th satellite and the i-th satellite and the i-th satellite is a fitness function; The transfer probability equation is: wherein, is the first only ant from satellite is transferred to satellite with a transfer probability, represents the next allowed transfer satellite set for the first only ant, is the pheromone concentration on the path between satellite and satellite before updating, represents the Q value predicted by the trained DQN model for transferring from satellite to satellite , is the heuristic value for transferring from satellite to satellite . The parameters of the ant colony optimization algorithm specifically include an ant quantity m, pheromone importance , heuristic importance , Q value importance , and pheromone evaporation rate ρ; Step three, searching based on the established fitness function and transfer probability equation to obtain a task planning result.

2. The method according to claim 1, wherein, The satellite Orbit transfer to the satellite Velocity change Is: Step one three one, record the orbit parameters of the satellite as Step one three two, record the orbit parameters of the satellite as ; in, It is a satellite The length of the semi-major axis of the elliptical orbit. It is a satellite The eccentricity of the orbit, It is a satellite The inclination angle of the orbital plane relative to the reference plane. It is a satellite The right ascension angle of the ascending node of the orbit. On satellite The position angle of an object on the track, It is a satellite Perigee angle of the orbit; It is a satellite The length of the semi-major axis of the elliptical orbit. It is a satellite The eccentricity of the orbit, It is a satellite The inclination angle of the orbital plane relative to the reference plane. It is a satellite The right ascension angle of the ascending node of the orbit. On satellite The position angle of an object on the track, It is a satellite Perigee angle of the orbit; Step one thirty two, satellite Instantaneous radius on its own orbit Is: Satellite Instantaneous radius on its own orbit is: Satellite Instantaneous velocity in its own orbit is: wherein is a gravitational parameter; Satellite Instantaneous velocity in its own orbit is: Step 133, calculate the semi-major axis of the Hohmann transfer orbit : Step one three four, calculating the speed of the Hohmann transfer orbit intersection point: wherein, is the velocity at the satellite's Hohmann transfer orbit intersection, is the velocity at the satellite's Hohmann transfer orbit intersection; Step one thirty-five, calculating satellite Velocity changes required for orbit and satellite Velocity changes required for orbit: wherein is a satellite the change in velocity required for the orbit, is a satellite the change in velocity required for the orbit; Step one thirty-six, calculate velocity change needed for orbital plane adjustment : Step one 37, calculating from satellite Orbit transfer to satellite Velocity change : 。 3. The method according to claim 2, wherein, In step one seven, the loss function is: wherein, is a loss function, denotes taking the expectation, denotes parameters of the secondary value network, denotes parameters of the primary value network.

4. The method according to claim 3, wherein, The working process of the main value network is as follows: In the main value network, the input successively passes through the input layer, the first gated recurrent unit, the second gated recurrent unit and the output layer, and the Q value is output through the output layer; The first gated recurrent unit comprises a reset gate and an update gate, and the working process of the first gated recurrent unit is as follows: Step 1, calculating the output of the reset gate : wherein, is the current input, is the previous hidden state, denotes the weight matrix of the reset gate, is the sigmoid activation function; Step 2, compute the output of the update gate : wherein, is an updated weight matrix for the gate; Step 3, calculating candidate hidden states : wherein, represents element-wise multiplication, represents an activation function, represents a weight matrix of the candidate hidden state; Step 4, update the hidden state : 。 5. The method according to claim 4, wherein, 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.

6. The method according to claim 5, wherein, The specific process of step three is as follows: Step three one, each ant randomly selects a satellite from the satellite constellation as a starting point; Step three two, each ant calculates the transfer probability of the current satellite to each satellite in the set of allowed transfer satellites according to the established probability transfer equation; Step three three, each ant selects the next satellite according to the maximum transfer probability corresponding to itself, and updates the pheromone concentration on the path; Step three four, each ant takes the satellite selected in step three three as the current satellite, and returns to step three two until the end of the whole search process, and the satellite access order sequence searched by each ant is obtained respectively; Step three five, the satellite access order sequence corresponding to the minimum fitness function value is taken as the final task planning result.

7. The method according to claim 6, wherein, The pheromone concentration updating mode is as follows: wherein is the satellite with the satellite updated pheromone concentration on the path between the satellite is the satellite with the satellite updated pheromone concentration on the path between the satellite denotes the number of the step amount of pheromone left by the scout ant on the path between the satellite and the satellite .

Images