Centroidal Voronoi Tessellation and Model Predictive Control – Based Macro-Micro Trajectory Optimization of Microsatellite Swarm

Probabilistic swarm guidance enables autonomous microsatellites to generate their individual trajectories independently so that the entire swarm converges to the desired distribution shape. However, it is essential to avoid crowding for reducing the possibility of collisions between microsatellites. To determine the collision-free guidance trajectory of each microsatellite from the current position to the target space, a collision avoidance algorithm is necessary. A synthesis method is proposed that generate the collision avoidance trajectories. The idea is that the trajectory planning is divided into macro-planning and micro-planning; macro-planning guides where the microsatellites move step by step from the initial cube to the target cube by probabilistic swarm guidance with Centroidal Voronoi tessellation, while the micro-planning is to generate the optimal path for each step and ﬁ nally reach the speci ﬁ ed position in the target cube by model predictive control. Simulation results are presented for the collision-free guidance trajectory of microsatellites to verify the bene ﬁ ts of this planning scheme.


Introduction
As an emerging multi-satellite cooperative flight mode, microsatellite swarm has become an important future research issue for distributed space systems due to their advantages of low cost, rapid response, and collaborative decision-making. Recently, reference [1] develops hundreds or thousands of small femtosatellites and establishes effective swarm GNC strategies for potential synthetic aperture applications. The problem of shaping a large swarm (10 4 -10 5 ) of particles into a large reflective dish was addressed [2][3], with only a few (40) actuators in the boundary of the workspace. Due to their small size, the microsatellites have limited sensing, actuation, and computation capabilities, which require the guidance and control algorithms of the swarm to be both fuel and computationally efficient.
To address the coordination of swarms for autonomous agents, a probabilistic guidance approach was investigated [4][5][6], which contain sub-swarms with different mission objectives. The probabilistic guidance approach is based on designing a Markov chain such that its steady-state distribution corresponds to the desired swarm density. On the basis of the above method, a homogeneous Markov chain-based approach was presented for the probabilistic density control of swarm. The proposed approach specifies the time evolution of the probabilistic density distribution by using a Markov chain, which guides the swarm to a desired steady-state distribution, while satisfying the prescribed ergodicity, motion, and safety constraints [7][8][9][10][11]. Different from the above method based on homogeneous Markov chain is that a novel and generic distributed swarm guidance algorithm using inhomogeneous Markov chain that guarantees superior performance over existing homogeneous Markov chain based algorithms, when the feedback of the current swarm distribution is available [12][13][14][15][16]. In contrast with previous homogeneous or inhomogeneous Markov chain-based approaches, optimal transport based approach was presented which guarantees faster convergence, minimizes a given cost function, and reduces the number of transitions for achieving the desired formation [2,3,17]. Moreover, the optimal transport based approach was verified by the Robotarium platform [2], and Jet Propulsion Laboratory of NASA has developed a guidance algorithm for a swarm of assets deployed from the back shell of the Mars spacecraft for distributed science [18]. All these cases indirectly show that the probabilistic guidance method can be applied to the practical swarm satellite mission.
The solution of the probabilistic guidance algorithm establishes the target bin for each agent. Collision avoidance algorithms are necessary for generating collision-free guidance trajectories from the present location to the target bin. The previous homogeneous or inhomogeneous Markov chain-based approaches adopt safety constraints to avoid collision, and the optimal transport-based approach adopts a distributed collision avoidance algorithm for multi agent path planning by modifying the coverage control algorithm. However, these collision avoidance algorithms are only macro-planning, thinking that the agent is a particle, without considering their micro-dynamics model. Therefore, it is difficult to guarantee that collision avoidance between microsatellites in physical space and that the trajectories may not be optimal and will waste fuel.
The main contribution of this paper is the integration of probabilistic swarm guidance with a macro-micro trajectoryplanning algorithm. Model predictive control (MPC) has been a major research area in aerospace for over a decade [19]. To save fuel consumption and improve computational efficiency, MPC has been used in applications similar to swarm guidance, such as spacecraft landing [20], spacecraft rendezvous [21], and reconfiguration of a satellite constellation [22]. Additionally, some collision avoidance algorithms by using MPC for spacecraft swarm are proposed [23,24]. The guidance trajectories obtained from PDGS can be implemented using sequential convex programming and model predictive control [25], and the experiment verified the feasibility of MPC in path planning for swarm.
In this paper, we propose a synthesis method that Centroidal Voronoi tessellation (CVT) [26][27][28] and model predictive control based macro-micro trajectory optimization of microsatellite swarm, and it is different from [25]. The main contributions of this work are listed as follows: (1) Comparing with the existing swarm distribution guidance [4][5][6][7][8][9][10][11][12][13][14][15][16][17], the physical space over which the swarm agents are distributed from 2D to 3D, it is more suitable for the practical application of swarm (2) A synthesis method is proposed that generate the collision avoidance trajectories. The idea is that the trajectory planning is divided into macro-planning and micro-planning, macro-planning guides where the microsatellites move step by step from the initial cube to the target cube by probabilistic swarm guidance with CVT, while the micro-planning is to generate the optimal path for each step and finally reach the specified position in the target cube by MPC (3) Considering J 2 perturbation in the micro-modeling of microsatellite swarm, it is more corresponded to engineering practice The remainder of this paper is organized as follows. The problem formulation is introduced in Section 2. The macroplanning of microsatellite swarm is presented in Section 3, and the micro-planning of microsatellites is presented in Section 4. The effectiveness of the proposed scheme is validated by using a numerical example in Section 5. Some conclusions and future work are given in Section 6.

Problem Formulation
2.1. Probabilistic Swarm Guidance. According to the probabilistic swarm guidance method described in literature [4], the agents can move in 2D plane. In this paper, we will investigate the transfer of swarm microsatellites in 3D space. The physical space over which the swarm microsatellites are distributed is denoted as R 3 ; it is assumed that the space are referred to as cubes.
Let rðtÞ be the position of a microsatellite at time index t. Let xðtÞ be a probability vector so that x½Ω c ðtÞ is the probability of the event that the microsatellite will be in cube R Ω c at time t, Consider a swarm consisting of N independently acting microsatellites so that (1) holds for N independent events: where r k ðtÞ denotes the position of the k-th microsatellite at time t. Here, x½Ω c ðtÞ is called the swarm distribution. The distribution guidance problem is defined as follows: Given any initial probability distribution, the swarm must be guided to a specified steady-state probability distribution v: Subject to motion constraints, specified by the adjacency matrix A a as The idea behind probabilistic swarm guidance is to control the propagation of probability vectors x½Ω c rather than the individual micorsatellite positions fr k ðtÞg N k=1 . The probability density vector xðtÞ over R 3 evolves as The probabilistic guidance algorithm is implemented by providing the Markov matrix M to each microsatellites. The probabilistic guidance algorithm can be referenced [4,5].

Safety Analysis of Collision Avoidance.
In the method of probabilistic swarm guidance, it is essential to avoid crowding for reducing the possibility of collisions between microsatellites. To determine the collision-free guidance trajectory of each microsatellite from the current position to the target space, a collision avoidance algorithm is needed. However, with high-level coordination that uses the macroscopic models, collision-free trajectories are very hard to generate. In this paper, a synthesis method was presented that the collision avoidance trajectories are generated. The idea is that the trajectory planning is divided into macro-planning and micro-planning, as shown in Figure 1, where macro-planning guides the microsatellites move step by step from the initial cube to the target cube, while the micro-planning is to generate the optimal path for each step and finally reach the specified position in the target cube.
Collision avoidance needs to be considered during trajectory generation. The minimum collision avoidance distance between all microsatellites is required, represented by R col . Collision avoidance analysis is based on finding the lower bound R col of the minimum distance between all microsatellites at any time. Suppose the distance between two microsatellites is jjη ij jj, so the difference between it and the minimum safe distance is γ ij ðη ij Þ = R col − kη ij k 2 , then The positive invariance condition can be expressed as When the distance between any microsatellites is less than the lower bound of R col , the time derivative of the dis-tance function can maintain or increase the distance, so collision avoidance can be guaranteed.

Macro-Planning of Microsatellite Swarm
is defined as the set of n distinct points in set Q. The Voronoi cell V i corresponding to the generation point p i is defined as follows: where k·k denotes the Euclidean distance in ℝ d . From the above definition, it can be seen that all points in Q whose Euclidean distance to the generation point p i is less than to the generation point p j≠i and the points are included in the Except for the boundary, all Voronoi cells are disjoint, and Q = S N i=1 V i is true, that is, the set V is a partition of the set Q.

Centroidal Voronoi Tessellation.
When the function of the target space is set, its generation point is the centroid of each Voronoi cell (the density function is used to calculate the centroid); it is the Centroidal Voronoi tessellation.
Assuming that the density function of space Q is defined as ρðqÞ ≥ 0, then the centroid of each Voronoi cell is positioned as When the position of the generation point is the centroid, namely: is Centroidal Voronoi tessellation. When the cost function is less than a certain threshold, it can be judged that the Centroidal Voronoi tessellation tends to be stable. Given the density function, the cost function is The cost function estimates the error between the generation point p i and the centroid of the Voronoi tessellation by the density function of the space.

CVT-Based Trajectory Planning for Microsatellite
Swarm. The target position of each microsatellite is determined by the centroid generated through CVT algorithm, and all microsatellites move to the corresponding centroid until the algorithm converges. In the present work, a code library of CVT based on the probabilistic Lloyd's method is incorporated to produce the CVT generators. According to the location of the centroid, the final distribution of the microsatellite swarm in the space is obtained. The main algorithm step is as follows: Determine 3D space Q, density function ρðqÞ ≥ 0, and define integer k microsatellites.
(1) Select the k points fp i g k i=1 randomly as generated points, namely, the initial position of the microsatellites (2) The probabilistic guidance algorithm [4,5] was executed to determine the cube where each microsatellite was located The reference frame is the local-vertical-local-horizontal (LVLH) denoted by F L . The origin of LVLH coordinate system is the centroid of the central microsatellite, the x-axis points to the centroid of the central microsatellite from the center of the earth, the z-axis is perpendicular to the orbital plane and points to the direction of the orbital angular momentum, and the y-axis is perpendicular to the x and z axes and is determined by the right hand rule.
The relative dynamics equation of the microsatellite considering the J 2 perturbation is where ½x y z T is the position vector of the microsatellite in the LVLH coordinate system of the reference satellite, ½u x u y u z T is the control torque imposed on microsatellite, n = ffiffiffiffiffiffiffiffi μ/r 3 p is the mean orbital motion of the microsatellite, r is geocentric distance, and μ is the Earth's gravitational constant.
The orbital elements are described as ða, e, ι, Ω, ω, νÞ, a is semi-major axis, e is eccentricity, ι is orbit inclination, Ω is the right ascension of ascending node, ω is argument of perigee, and ν is true anomaly. ½ f x f y f z T is the perturbing acceleration caused by J 2 : where the coefficient κ = −3/2μJ 2 ðR e 2 /r 4 Þ, R e is the radius of the Earth, J 2 is the second zonal harmonic coefficient of the Earth, and u = ω + ν is argument of latitude.
The relative dynamics model of microsatellite swarm is expressed in state space: where N is the number of microsatellites, and the state of the j-th microsatellite at time t is x j ðtÞ, x j = ½x j y j z j _ x j _ y j _ z j T , and u j = ½u jx u jy u jz T . Meanwhile, To satisfy the conditions of solving convex optimization problem, the dynamics model was discretized and obtained from the continuous dynamics model by zero-order equal interval sampling: Then, the dynamics model of discretization is

Convexification of Collision Avoidance Constraints.
To avoid two members of microsatellite swarm moving to the same location at the same time, collision trajectory between members should be considered in the process of swarm reconfiguration. Then, the anti-collision avoidance constraint can be expressed as where and R com is the maximum communication distance. Equation (21) is linearized at the discrete points corresponding to the two reference trajectories of microsatellite i and j, and the affine function about the state of the two adjacent microsatellites can be obtained, so as to realize the convexity of collision avoidance constraints between microsatellites: where x j is the nominal trajectory of each trajectory point. The accurate collision avoidance constraint requires only that the distance between two adjacent microsatellites is no less than the safe distance. However, after convexization, as shown in Figure 2, the constraint becomes that the two adjacent microsatellites cannot simultaneously be located in the belt region with width R col , which is perpendicular to the connection line of the discrete reference points of the two adjacent microsatellites, and the position of the space is not fixed, so it can move along the connection direction of the reference points.

MPC-Based Trajectory Optimization for Microsatellites.
To implement real-time trajectory planning, model predictive control is introduced; the MPC implementation uses a receding horizon to update the optimal trajectories based on the current state information. Assuming that the prediction horizon and the control horizon are both p, according to the principle of model predictive control, it can be obtained wherex k is initial state, and The cost function is defined as where Q, R, and F are the weight matrix, x k+i is the predicted state, x k+p is the predicted terminal state, Q and R are usually the diagonal matrix, and F is the terminal cost weight matrix obtained by solving the algebraic Riccarti equation. where In the whole reconfiguration stage, it is assumed that microsatellite members of swarm have limited communica-tion capacity, so the members can only communicate with their neighboring microsatellites, so each member microsatellite needs to consider anti-collision constraints. So the trajectory planning problem of microsatellite swarm considering collision constraints can be described as

Results
In this section, numerical simulation is carried out to verify the proposed macro-micro trajectory planning method of microsatellite swarm. In this paper, a virtual central microsatellite is given, and a large-scale (300) microsatellite swarm is designed with an omnidirectional flight configuration. The initial orbital element of the central satellite is (7106.14 km, 0.05°, 98.3°, 270°, 0°, 0°). The earth's gravitational constant is 3:986 × 10 5 km 3 s −2 ; the radius of the earth is 6378.140 km. In LVLH coordinates, the position vectors of 300 microsatellites relative to the virtual center microsatellite are randomly generated, the azimuth angle is randomly selected, and the initial distance is randomly selected within 100 km.
In the micro-planning, we do not have a reference trajectory only a reference final position. Hence, a terminal  Space: Science & Technology constraint is added to the MPC. The state error cost is kept low because, we do not have a state trajectory to follow only the final destination. Meanwhile, the control cost is kept high to emphasize the optimization of fuel consumption to be minimum to the reference which is taken to be zero. The sampling time is 60s, the prediction horizon is 50. The Q is 1 × 10 −11 * I 6×6 , R is 1 * I 3×3 , and F is 1 × 10 −7 * I 3×3 .
The overall configuration transformation of microsatellite swarm is shown in Figures 3 and 4. The initial configuration generated is shown in Figures 3(a) and 3(d); the desired configuration design is of type "0". Here, the numbers of microsatellites in each sub-swarm have been specifi-cally chosen so that the final configuration is uniform over the shape of the "0". With the increase of iteration times, microsatellite swarm gradually forms the desired configuration under the action of macro and micro trajectory planning. It can be seen that the desired configuration can be formed quickly. After the 10 iterations, the desired configuration is basically formed and maintained in a stable state; as shown in Figure 4, there is no significant change in distance scale. In the meanwhile, the total variation in Figure 5 indicates that convergence is achieved at roughly t = 10 steps.
In this paper, we use the CVT algorithm to divide regions, and so as to determine the position of the   Space: Science & Technology microsatellites to be transferred at the next moment. As shown in Figure 6, we selected one of the cubes in the transfer process and perform CVT on it to determine the transfer position of the microsatellite. After 50 iterations, a stable configuration is obtained, and the position where the microsatellite moves at the next moment is determined. Due to the large scale of the microsatellite swarm, the process of achieving the final configuration requires many transitions. To verify the proposed trajectory optimization based on model predictive control, one of the microsatellites is selected from the initial point to the next desired target point at a certain moment, as shown in Figure 7. The individual microsatellites can reach the desired point well. After the desired point is reached, the next iteration will be carried out, and due to the influence of orbital dynamics, the microsatellite may not remain the target point without control constraints. Therefore, the trajectory curve after T = 6000 s of the simulation in Figure 7 is not a reference. To make the mission of microsatellite swarm more practical, we used MPC in micro-planning to improve the performance of microsatellite swarm in terms of fuel consumption and resource utilization.

Conclusion
This paper developed a synthesis trajectory planning method for microsatellite swarm, which include macro-planning and micro-planning. The former used Centroidal Voronoi tessellation, while the latter adopted MPC to generate the optimal trajectories. The proposed method can not only realize collision avoidance of microsatellite swarm maneuvering at the macrolevel, but also provide optimal trajectories for each microsatellite of swarm individuals at the micro-level. Macro-planning provides guidance for micro-planning, and micro-planning implements the commands of macro-planning; the proposed method accords well with engineering practice. Future work will focus on the extending the presented scheme to achieve a faster, more fuel-efficient trajectory optimization method. Moreover, the validation experiments will be carried out.

Data Availability
The simulation date and the program files in this paper cannot be shared with others as they are our basis for next research.