Satellite Cluster Formation Reconfiguration Based on the Bifurcating Potential Field

Abstract

The satellite cluster formation reconfiguration has received considerable attention in recent years. However, the traditional centralized control methods are challenging to apply to satellite clusters because of the enormous fuel consumption, and few studies have addressed the mathematical characterization of satellite clusters. This research aims to propose a mathematical characterization method for satellite clusters and investigate the formation reconfiguration control of satellite clusters. This study provided the five-element characterization method to represent the cluster characteristics and internal correlation characteristics of orbiting satellite clusters. In addition, a control method based on bifurcating potential fields was proposed to realize satellite cluster formation’s dynamic migration and rapid reconfiguration. A cluster with 50 satellites was simulated to verify the feasibility and effectiveness of the proposed formation control algorithm. The results show that various formation topologies were achieved by simply changing the bifurcation parameter and configuration adjustment parameters. The five descriptive elements of the satellite cluster can intuitively and effectively reflect the running state of the satellite cluster.

1. Introduction

In recent years, the application of low-cost micro-nano satellite clusters with good expansibility and flexibility to perform space missions has received considerable attention. A satellite cluster is a satellite system composed of several micro-nano satellites to complete space missions. The satellites in the cluster operate loosely in close proximity within a bounded space area and do not need to maintain strict spatial geometry configuration. Most satellite clusters are heterogeneous, relying on local information interaction to maintain the relative motion bounded [1,2,3], which are loose clusters that accomplish space tasks through autonomous cooperation. At present, the National Aeronautics and Space Administration (NASA), European Space Agency (ESA), and other space agencies have developed or planned several satellite cluster systems for detection, remote sensing, communication, and surveillance. Among them, the KickSat project, the Flock earth observing constellation, the Starlink constellation, the BlackJack program, the Starling mission, the SWIFT project, the ANTS mission, the OLFAR project, the Tiantuo-3 mission, and the SULFRO mission are more representative [4,5,6,7,8,9,10,11,12,13].

Mathematical characterization is fundamental to the study of satellite clusters. The basic formation configurations in traditional satellite formation flying (SFF), such as accompanying formation and flying-around formation, can be accurately represented by relative motion equations. Satellite cluster is derived from the traditional SFF, and its operation is based on the relative movement of satellites within the cluster. However, unlike the conventional SFF, satellite cluster does not need to maintain strict spatial geometry configuration. The cluster has some characteristics, which are different from the individual. Therefore, the description of a large-scale satellite cluster should not be limited to the relative motion between two satellites. At present, an effective mathematical characterization method is needed to describe the overall state of the satellite cluster and the relationships among satellites within the cluster.

Formation reconfiguration is the crucial technique for the application of satellite clusters. Given the volume and power limitation of micro-nano satellites [14] and the loose characteristic of satellite clusters, some traditional centralized control methods are challenging to apply to satellite clusters because of the enormous fuel consumption. Although distributed control methods, such as the artificial potential field (APF) method and the behavior-based control method, have achieved good results in robots and Unmanned Aerial Vehicles (UAVs), these methods cannot be directly applied to satellite formation control. It is relatively easy to control the UAV in three axes. However, the satellite operates in an exceptional space environment, so it is necessary to consider the difference in fuel consumption between in-plane and out-of-plane orbit controls.

The primary purpose of this research is to propose a mathematical characterization method for satellite clusters and investigate the formation reconfiguration control of satellite clusters based on the bifurcating potential field.

The APF method was first introduced by Khatib for mobile robot path planning and obstacle avoidance [15]. The basic principle of the APF method is to build artificial potential functions according to the surrounding environment of the controlled object while making the expected state of the managed object located at the global minimum point of the system’s potential field. By setting the control law, the managed object can move along the negative gradient direction of the potential field until it converges to the expected state. The APF method is very close to the sliding mode control and the Lyapunov function-based feedback control. All of them can guide the system to be stable in the desired state and have strong robustness. However, the APF method focuses more on path planning, which can easily maintain the relationship between satellites, as well as between satellites and the environment, such as collision avoidance and obstacle avoidance. At present, the APF method has been successfully applied to the multi-robot system and multi-UAV system and is gradually used for spacecraft rendezvous, in-orbit self-assembly, and spacecraft obstacle avoidance [16,17,18,19,20,21,22,23,24]. The name “bifurcation” was first introduced by Henri Poincaré and later commonly applied to the mathematical study of dynamical systems [25]. A bifurcation occurs when a small smooth change made to the bifurcation parameter values causes a sudden qualitative or topological change in its behavior [26]. The mathematical expression of the bifurcation function is generally straightforward. However, the system properties can be changed by the simple change of bifurcation parameters, which fit in with the formation reconfiguration of satellite clusters. Therefore, we consider combining the APF method with bifurcation theory to achieve satellite cluster formation reconfiguration.

The paper proceeds as follows. In Section 2, the five-element characterization method is proposed to describe the characteristics of orbiting satellite clusters. The satellite cluster is modeled, explaining the APF method and bifurcation theory. Section 3 shows the simulation results of satellite cluster formation reconfiguration, using the novel control law based on bifurcating potential field, and analyzes the regulating effect of the parameters in the formation control algorithm. Section 4 presents the conclusions and applications of this research.

2. Materials and Methods

2.1. The Five Descriptive Elements of Satellite Cluster

Satellite cluster formation control needs to consider constraints such as collision avoidance, obstacle avoidance, and maximum communication range. Since the satellites in the cluster operate nearby, the collision probability is closely related to the distribution within the cluster. The maximum communication range of the cluster is closely related to the scale and outer boundary of the cluster. The inter-satellite relationship within the cluster and the overall motion state of the whole cluster should be considered. This paper proposes the five-element characterization method to describe the characteristics of orbiting satellite clusters. In this method, a set of five basic characteristic parameters with specific geometric significance are used to intuitively characterize the satellite cluster from the aspects of cluster size, cluster movement, spatial cluster structure, and internal cluster distribution.

2.1.1. Cluster Size

A satellite cluster is a satellite system composed of several satellites, which operate in a bounded space area to complete space missions. Similar to the characteristics of the biological cluster, satellite clusters of different sizes show different clustering capabilities. Therefore, the cluster size parameter is defined as the number of member satellites in the satellite cluster from a macro perspective, denoted as $N$, which reflects the size of the satellite cluster to a certain extent.

2.1.2. Cluster Movement

Each member satellite in the cluster operates in its orbit, but the satellite cluster as a whole appears as a bounded cluster. Each member satellite’s running state determines the group motion state of the cluster. The satellite cluster can be regarded as a sizeable virtual satellite, as shown in Figure 1. The satellite’s position vector and velocity vector in the earth-centered inertial frame correspond to the six orbital elements, which can accurately describe the spatial geometric characteristics of the satellite’s orbit and the satellite’s current position in orbit. Therefore, the cluster movement parameter is defined as the average value of all member satellites’ position vector and velocity vector in the earth-centered inertial frame, denoted as $\left({\mathit{r}}_o,{\mathit{v}}_o\right)$, which reflects the current spatial orientation and orbital operation of the whole satellite cluster. The cluster movement parameter is mathematically defined as:

$${\mathit{r}}_o=\frac{{\displaystyle \sum _1^n{\mathit{r}}_i}}{N},$$

(1)

$${\mathit{v}}_o=\frac{{\displaystyle \sum _1^n{\mathit{v}}_i}}{N},$$

(2)

where ${\mathit{r}}_i$ and ${\mathit{v}}_i$ are the position and velocity of any member satellite in the earth-centered inertial frame in units of km and km/s.

2.1.3. Outer Boundary of the Cluster

The ${\mathit{r}}_o$ in the cluster movement parameter describes the geometric center position of the cluster. On this basis, it is also necessary to determine the outer boundary of the satellite cluster, which can be understood as determining the space volume of the large virtual satellite, since the maximum communication range constraint is closely related to the outer boundary of the cluster, which is determined by the satellite farthest from the geometric center of the cluster. The outer boundary of the cluster is defined as the maximum distance between each member satellite and the geometric center of the satellite cluster, denoted as $d_{\max }$, which reflects the space range occupied by the cluster. The member satellites remain in the spherical region with this parameter as the radius. A schematic diagram of the outer boundary is shown in Figure 2. This parameter can intuitively judge whether the cluster meets the expected communication range constraint. The outer boundary of the cluster is mathematically defined as:

$$d_{\max }=\begin{array}{cc}\max (d_{io})& (1\le i\le n),\end{array}$$

(3)

where $d_{io}$ is the distance from any member satellite to the geometric center of the satellite cluster in units of km.

2.1.4. Inner Boundary of the Cluster

The description of satellite cluster configuration needs to consider the overall behavior of the cluster and the individual behavior within the cluster. Due to the member satellites being close to each other, it is necessary to pay attention to the inner boundary to make collision avoidance warnings in time and ensure the safe operation of the satellite cluster. The inner boundary of the cluster is defined as the minimum distance between any two satellites in the satellite cluster, denoted as $d_{\min }$, which reflects the distance relationship between satellites in the cluster. A schematic diagram of the inner boundary is shown in Figure 3. This parameter can intuitively judge whether the cluster meets the expected collision avoidance constraint. The inner boundary of the cluster is mathematically defined as:

$$d_{\min }=\begin{array}{cc}\min (d_{ij}-a_i-a_j)& (1\le i\le N,1\le j\le N,i\ne j),\end{array}$$

(4)

where $d_{ij}$ is the distance between the centroids of any two member satellites in the satellite cluster, in units of km, and $a_i,a_j$ are the space radius of satellite $i$ and satellite $j$ respectively, in units of km.

2.1.5. Distribution Uniformity of the Cluster

In the non-mission state, the satellite cluster runs loosely within a specific space like bees or birds. In the mission state, the satellite cluster forms unique formation configurations to improve the efficiency of mission execution. Under different formation configurations of the same satellite cluster, the distribution of member satellites is different and has specific rules to follow. For example, when the outer boundary of the cluster remains constant, the distribution uniformity of the linear configuration is generally lower than that of the disk configuration. Distribution uniformity of the cluster, denoted as $L$, is equal to the total volume of each member satellite’s monopolized sphere divided by the volumes of the cluster’s monopolized sphere. $L$ is mainly defined to judge the formation configuration of the satellite cluster based on its current distribution uniformity, representing the uniformity of the cluster in most situations. The distribution uniformity of the cluster is mathematically defined as:

$$L=\frac{{\displaystyle \sum _1^n{V}_i}}{V},$$

(5)

where $V_i$ is the volume of any member satellite’s monopolized sphere, and $V$ is the volume of a circumscribed sphere of all satellites’ monopolized spheres.

In this research, the satellite’s monopolized sphere is defined as follows. Denote the distance between satellite $i$ and its nearest neighbor as $d_i$. Then, the center of satellite $i$’s monopolized sphere is itself, and the radius of that monopolized sphere is $d_i/2$. A schematic diagram of the monopolized sphere of each satellite in the satellite cluster is shown in Figure 4.

The specific calculation steps of $L$ are as follows:

• The radius of each member satellite’s monopolized sphere is calculated according to the distance between each member satellite:

$$m_i=\frac{1}{2}\min (d_{ij})(1\le j\le N,i\ne j)$$

(6)

2

Calculate the volume of each member satellite’s monopolized sphere:

$$V_i=\frac{4}{3}\pi m_i^3$$

(7)

3

Calculate the radius of circumscribed sphere of all satellites’ monopolized spheres:

$$M=\frac{1}{2}\max (d_{ij}+m_i+m_j)(1\le i\le N,1\le j\le N,i\ne j)$$

(8)

4

Calculate the volume of the circumscribed sphere:

$$V=\frac{4}{3}\pi M^3$$

(9)

5

Finally, the distribution uniformity of the cluster is calculated by Equation (5).

2.2. Formation Control Algorithm

2.2.1. Dynamic Model of Satellite Cluster

Consider a satellite cluster composed of $n$ satellites, the motion equation of member satellite $i$ in the satellite cluster can be expressed as:

$$\{\begin{array}{l}{\dot{\mathit{x}}}_i={\mathit{v}}_i,\hfill \\ m_i{\dot{\mathit{v}}}_i={\mathit{F}}_i={\mathit{F}}_{ic}+{\mathit{F}}_{if}-k_v{\mathit{v}}_i,\hfill \end{array}$$

(10)

where ${\mathit{x}}_i$ and ${\mathit{v}}_i$ are the position and velocity of satellite $i$ in the earth-centered inertial frame, $m_i$ is the mass of satellite $i$, and ${\mathit{F}}_i$ is the resultant force on satellite $i$. ${\mathit{F}}_i$ is composed of collision avoidance control force ${\mathit{F}}_{ic}$, formation control force ${\mathit{F}}_{if}$, and velocity feedback $-k_v{\mathit{v}}_i$, where $k_v$ is a positive control gain to control the amplitude of the dissipation in units of kg/s.

The APF method is widely used in the formation control by building artificial potential functions to create virtual attraction or repulsion potential fields, similar to the electric fields in electromagnetics. The superposition of each potential field forms the system’s potential field. Satellite under the system potential field will move towards the global minimum point along the negative gradient of the system’s potential field. Two potential fields are considered in this research. One is collision avoidance potential field aim to create ${\mathit{F}}_{ic}$, and the other is bifurcating potential field used to create ${\mathit{F}}_{if}$.

Due to the relatively compact structure of the satellite cluster, the relative motion between the satellites can be described by the Clohessy–Wiltshire (CW) equation. The CW equation, also known as Hill equation, was first proposed by George William Hill [27] when studying the motion of the moon. Hill equation studies the motion law of the lunar relative to the Earth. The CW equation is a linearized equation proposed by Clohessy et al. [28] to describe the close relative motion of two satellites, which is expressed in the target orbital frame (defined in Figure 5) as:

$$\{\begin{array}{l}\ddot{x}-2\omega \dot{z}=u_x,\hfill \\ \ddot{y}+{\omega }^{2}y=u_y,\hfill \\ \ddot{z}-3{\omega }^{2}z+2\omega \dot{x}=u_z,\hfill \end{array}$$

(11)

where $x,y,z$ are the coordinates of the satellite in the target orbital frame, $\omega$ is the mean angular velocity of the target orbital frame, and $u_x,u_y,u_z$ are the accelerations applied by the satellite thrusters in the $x,y,z$ directions.

2.2.2. Collision Avoidance Potential Function

The collision avoidance potential field has two effects. One is to disperse satellites during the formation configuration process, and the other is to avoid collisions when the inter-satellite distance is less than the safe value. The magnitude of the repulsive force is related to the inter-satellite distance. The repulsive force is zero when the inter-satellite distance is within the safe range. When the inter-satellite distance is outside the safe range, the repulsive force increases as the inter-satellite distance decreases.

The inter-satellite collision avoidance potential function is established based on the inter-satellite distance as follows:

$$U_{ij}=\{\begin{array}{cc}-k_c\left({\mathit{x}}_{ij}^T{\mathit{x}}_{ij}-s\Vert {\mathit{x}}_{ij}\Vert \right)& if\Vert {\mathit{x}}_{ij}\Vert \le s,\\ 0& if\Vert {\mathit{x}}_{ij}\Vert >s,\end{array}$$

(12)

where $U_{ij}$ is the avoidance potential between satellite $i$ and satellite $j$, $k_c$ is a positive control gain to adjust the amplitude of the potential function in units of $\mathrm{kg}/{\mathrm{s}}^2$, $s$ is a constant describing the desired safe distance between satellite $i$ and satellite $j$, and ${\mathit{x}}_{ij}$ is the position of satellite $i$ relative to satellite $j$ in units of km.

The corresponding repulsive force of satellite $i$ by satellite $j$ is:

$${\mathit{F}}_{ij}=-{\nabla }_{{\mathit{x}}_{ij}}{U}_{ij}=\{\begin{array}{cc}{k}_c\left(1-\frac{s}{\Vert {\mathit{x}}_{ij}\Vert }\right){\mathit{x}}_{ij}& if\Vert {\mathit{x}}_{ij}\Vert \le s.\\ 0& if\Vert {\mathit{x}}_{ij}\Vert >s.\end{array}$$

(13)

The collision avoidance control force of satellite $i$ is equal to the sum of the repulsive force of satellite $i$ by all other member satellites in the cluster:

$${\mathit{F}}_{ic}={\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^n{\mathit{F}}_{ij}}.$$

(14)

The repulsion force ${\mathit{F}}_{ic}$, the formation control force ${\mathit{F}}_{if}$, and velocity feedback $-k_v{\mathit{v}}_i$ together guide satellites toward the goal state. The repulsive force is used to ensure that satellites do not collide with each other when they are steered towards the goal state. Once all satellites are driven to the goal equilibrium state, the repulsive force also ensures that they are equally spaced for symmetric formations.

Let $k_c=1 \mathrm{kg}/{\mathrm{s}}^2,s=0.05 \mathrm{km}$, the collision avoidance potential function and the collision avoidance control force are shown in Figure 6.

2.2.3. Bifurcating Potential Function

Bifurcation theory is a mathematical study of qualitative or topological structure changes of a given family of curves, often used in the mathematical study of dynamical systems. A continuous-time dynamic system can be described as:

$$\dot{x}=f\left(x,\mu \right),$$

(15)

where $x$ is the state variable, $x=x(t)$, $\dot{x}=dx/dt$, and $\mu$ is the bifurcation parameter. The independent variable $t$ represents the time, and it is usually omitted as in Equation (15) for simplicity.

Define a bifurcation according to the requirements of satellite cluster formation reconfiguration, and its typical form is:

$$\dot{x}=f\left(x,\mu \right)=x\left(\mu -\left|x\right|\right).$$

(16)

The equilibrium points of the system are:

$$x_e=\{\begin{array}{cc}0& if \mu <0.\\ 0,\pm \mu & if \mu >0.\end{array}$$

(17)

Linearize the system described by Equation (16) at the equilibrium points:

$$\dot{x}=Ax.$$

(18)

$$A={\frac{\partial f}{\partial x}|}_{x=x_e}={\left(\mu -2\left|x\right|\right)|}_{x=x_e}.$$

(19)

The eigenvalues of $A$ are:

$$\lambda =\{\begin{array}{ll}\mu \hfill & if \mu <0,x_e=0.\hfill \\ \mu \hfill & if \mu >0,x_e=0.\hfill \\ -\mu \hfill & if \mu >0,x_e=\pm \mu .\hfill \end{array}$$

(20)

The stability of each equilibrium point can be judged based on Lyapunov theory. When $\mu <0$, the equilibrium point $x_e=0$ is stable. When $\mu >0$, the equilibrium point $x_e=0$ is unstable, while $x_e=\pm \mu$ are stable. Based on this bifurcation, the bifurcating potential field is constructed to obtain a stable cluster configuration. The bifurcating potential function is expressed as:

$$U=\frac{{\left|x\right|}^3}{3}-\mu \frac{x^2}{2}.$$

(21)

Figure 7 shows the shape of the bifurcating potential when $\mu =\pm 5$. The minimum point of bifurcating potential function corresponds to the stable equilibrium point of bifurcation dynamics.

We propose a two-dimensional bifurcating potential function for the in-plane formation control of satellite clusters:

$$U_f=k_f\left(\frac{{\left|\sqrt{b^2{x}^2+a^2{z}^2}-abr\right|}^3}{3}-\mu \frac{{\left(\sqrt{b^2{x}^2+a^2{z}^2}-abr\right)}^2}{2}+\sigma y^2\right),$$

(22)

where $k_f$ is a positive control gain to adjust the amplitude of the potential function in units of $\mathrm{kg}/\mathrm{km}{\mathrm{s}}^2$, $a,b$ are dimensionless configuration adjustment parameters, $r$ is another configuration adjustment parameters in units of $\mathrm{km}$, $\mu$ is the bifurcation parameter in units of $\mathrm{km}$, $\sigma$ is the out-of-plane potential parameter in units of $\mathrm{km}$, and $x,y,z$ are the coordinates of satellite $i$ relative to the desired formation center in units of $\mathrm{km}$, denoted as $x_i=[\begin{array}{ccc}x& y& z\end{array}]$. $\sigma y^2$ is the out-of-plane potential used to control the satellite cluster operates in the orbital plane.

The corresponding formation control force of satellite $i$ under this two-dimensional bifurcating potential function is:

$${\mathit{F}}_{if}=-{\nabla }_{{\mathit{x}}_i}{U}_{if}.$$

(23)

According to the two-dimensional bifurcating potential function, a variety of satellite cluster formation configurations can be formed in the orbital plane. Each configuration has clear correspondence with the bifurcating parameters. The specific corresponding relationships are shown in

Table 1. In-plane formation configurations and bifurcating parameters.

In-Plane Formation Configuration	$\mathit{a}$	$\mathit{b}$	$\mathit{r} (\mathbf{km})$	$\mathit{\mu } (\mathbf{km})$
circle (radius r)	$1$	$1$	$r$	$<0$
$\mathrm{concentric} \mathrm{double} \mathrm{circle} (\mathrm{radius} r\pm \mu$)	$1$	$1$	$r$	$>0$
disk	$1$	$1$	$0$	$<0$
ellipse semi-major/minor axis ar/br)	$a$	$b$	$1$	$<0$
concentric double ellipse (semi-major/minor axis $ar\pm \mu /br\pm \mu$)	$a$	$b$	$1$	$>0$
ellipse disk	$a$	$b$	$0$	$<0$
$\mathrm{line} (\mathrm{formation} \mathrm{along} z=0$)	$a$	$0$	$r$	$<0$
$\mathrm{double} \mathrm{line} (\mathrm{formation} \mathrm{along} z=\pm \mu /a$)	$a$	$0$	$r$	$>0$

.

The bifurcating potential field shows different equilibrium states depending on its varying bifurcating parameters. Let $k_f=1 \mathrm{kg}/\mathrm{km}{\mathrm{s}}^2,a=b=1$, three different potential field surfaces are shown in Figure 8 when $r=0.6 \mathrm{km},\mu =-0.1 \mathrm{km}$, $r=0.6 \mathrm{km},\mu =0.4 \mathrm{km}$, and $r=0 \mathrm{km},\mu =-0.1 \mathrm{km}$ respectively. Satellites under the bifurcating potential field will move towards the global minimum point along the negative gradient of the bifurcating potential field. Under the circle potential field shown in Figure 8a, the satellites finally converge to the circle with radius $r$ centered on the reference point, forming a circle configuration. Under the concentric double circle potential field shown in Figure 8b, the satellites finally converge to the concentric double circle with radius $r\pm \mu$ centered on the reference point, forming the concentric double circle configuration. Under the disk potential field shown in Figure 8c, the satellites finally converge near the reference point and form a disk configuration under the collision avoidance potential field.

Let $k_f=1 \mathrm{kg}/\mathrm{km}{\mathrm{s}}^2,a=1,b=0.5$; three different potential field surfaces are shown in Figure 9 when $r=1 \mathrm{km},\mu =-1 \mathrm{km}$, $r=1 \mathrm{km},\mu =0.4 \mathrm{km}$, and $r=0 \mathrm{km},\mu =-1 \mathrm{km}$ respectively. Satellites under these potential field will move towards the global minimum point, forming the ellipse configuration with semi-major/minor axis $ar/br$, the concentric double ellipse configuration with semi-major/minor axis $ar\pm \mu /br\pm \mu$, and the ellipse disk configuration centered on the reference point, respectively.

Let $k_f=1 \mathrm{kg}/{\mathrm{kms}}^2,a=1,b=0$, two different potential field surfaces are shown in Figure 10 when $\mu =-1 \mathrm{km}$ and $\mu =0.5 \mathrm{km}$, respectively. In this case, $r$ loses the role of formation configuration adjustment. Satellites under these potential field will move towards the global minimum point, forming the line configuration along $z=0$ and the double line configuration along $z=\pm \mu /a$, respectively.

We propose a three-dimensional bifurcating potential function for the stereo formation control of satellite clusters:

$$U_{if}=k_f\left(\frac{{\left|d-abcr\right|}^3}{3}-\mu \frac{{\left(d-abcr\right)}^2}{2}\right)$$

(24)

$$d=\sqrt{b^2{c}^2{x}^2+a^2{c}^2{y}^2+a^2{b}^2{z}^2}$$

(25)

where $k_f$ is a positive control gain to adjust the amplitude of the potential function in units of $\mathrm{kg}/\mathrm{km}{\mathrm{s}}^2$, $a,b,c$ are dimensionless configuration adjustment parameters, $r$ is another configuration adjustment parameters in units of $\mathrm{km}$, $\mu$ is the bifurcation parameter in units of $\mathrm{km}$, and $x,y,z$ are the coordinates of satellite $i$ relative to the desired formation center in units of $\mathrm{km}$, denoted as ${\mathit{x}}_i=[\begin{array}{ccc}x& y& z\end{array}]$.

The corresponding formation control force of satellite $i$ under this three-dimensional bifurcating potential function is:

$${\mathit{F}}_{if}=-{\nabla }_{{\mathit{x}}_i}{U}_{if}.$$

(26)

According to the three-dimensional bifurcating potential function, a variety of satellite cluster formation configurations can be formed, and each configuration has clear correspondence with the bifurcating parameters. The specific corresponding relationships are shown in

Table 2. Stereo formation configurations and bifurcating parameters.

Stereo Formation Configuration	$\mathit{a}$	$\mathit{b}$	$\mathit{c}$	$\mathit{r} (\mathrm{km})$	$\mathit{\mu } (\mathrm{km})$
spherical surface (radius r)	$1$	$1$	$1$	$r$	$<0$
$\mathrm{concentric} \mathrm{double} \mathrm{spherical} \mathrm{surface} (\mathrm{radius} r\pm \mu$)	$1$	$1$	$1$	$r$	$>0$
spherical space	$1$	$1$	$1$	$0$	$<0$
ellipsoidal surface (semi-axis $ar/br/cr$)	$a$	$b$	$c$	$1$	$<0$
concentric double ellipsoidal surface (semi-axis $ar\pm \mu /br\pm \mu /cr\pm \mu$)	$a$	$b$	$c$	$1$	$>0$
ellipsoidal space	$a$	$b$	$c$	$0$	$<0$

.

3. Results and Discussion

Several cases of satellite cluster formation configuration and reconfiguration are simulated to test the bifurcating control law. The satellite cluster is composed of 50 satellites with a weight of 100 kg. In this paper, the control law was designed based on the CW equation because the satellites move in close proximity (tens to hundreds of meters). In the numerical simulation, the Runge–Kutta Fehlberg 78 method is used to propagate the satellite orbit, and the influence of J2 term is considered.

3.1. Formation Configuration

3.1.1. Line and Double Line Formation

A satellite cluster composed of 50 satellites with a weight of 100 kg initially operates in spherical airspace with a radius of 0.1 km. The cluster center runs on a nearly circular orbit with an altitude of 600 km, an inclination of 98°, and an eccentricity of 0.0002. At the initial time, the satellite cluster receives the control command. It is required to establish an in-plane line formation along $z=0 \mathrm{km}$. When $a\ne 0,b=0,\mu <0$, the satellite cluster can establish a line formation along $z=0 \mathrm{km}$. The simulation control parameters were set as follows:

$$k_c=0.25 \mathrm{kg}/{\mathrm{s}}^2,k_f=10 \mathrm{kg}/{\mathrm{kms}}^2,k_v=2000 \mathrm{kg}/\mathrm{s},a=1,b=0,r=1 \mathrm{km},s=0.2 \mathrm{km},\mu =-1 \mathrm{km}.$$

When $a\ne 0,b=0,\mu >0$, the satellite cluster can establish a double line formation along $z=\pm \mu /a$. Change the control command to establish an in-plane double line formation along $z=\pm 0.5 \mathrm{km}$. The simulation control parameters were set as follows:

$$k_c=0.25 \mathrm{kg}/{\mathrm{s}}^2,k_f=10 \mathrm{kg}/{\mathrm{kms}}^2,k_v=2000 \mathrm{kg}/\mathrm{s},a=1,b=0,r=1 \mathrm{km},s=0.2 \mathrm{km},\mu =0.5 \mathrm{km}.$$

Figure 11 shows the motion trajectory of the satellite cluster to establish the line formation along $z=0 \mathrm{km}$ and the double line formation along $z=\pm 0.5 \mathrm{km}$ in the target orbital frame. The satellites gather in spherical airspace with a radius of 0.1 km at the initial time. Under the formation control law, each satellite can spontaneously choose its motion path according to its relative state in the potential field.

3.1.2. Circle Formation

A satellite cluster composed of 50 satellites with a weight of 100 kg initially operates in spherical airspace with a radius of 0.1 km. The cluster center runs on a nearly circular orbit with an altitude of 600 km, an inclination of 98°, and an eccentricity of 0.0002. At the initial time, the satellite cluster receives the control command. It is required to establish an in-plane circle formation with a radius of 0.5 km, taking the cluster center as the formation center. The simulation control parameters were set as follows:

$$k_c=0.25 \mathrm{kg}/{\mathrm{s}}^2,k_f=3 \mathrm{kg}/\mathrm{km}{\mathrm{s}}^2,k_v=2000 \mathrm{kg}/\mathrm{s},a=b=1,r=0.5 \mathrm{km},s=r,\mu =-5 \mathrm{km}.$$

Figure 12 shows the motion trajectory of the satellite cluster to establish an in-plane circle formation with a radius of $r$ in the target orbital frame. The satellites gather in spherical airspace with a radius of 0.1 km at the initial time. Under the formation control law, each satellite can spontaneously choose its motion path according to its relative state in the potential field. In 0–1000 s, the in-plane circle formation is primarily established by the bifurcating potential function, and in 1000–2000 s, the formation configuration is evenly distributed by the collision avoidance potential function.

Figure 13 shows the three-axis velocity increments of each satellite in the process of establishing an in-plane circle formation. The velocity increment of each satellite is constrained within 1 m/s, which is in line with the actual application of micro-nano satellites.

Figure 14 shows that when the target formation is established and evenly distributed, the distribution uniformity of the cluster, $L$, tends to the steady-state value of 0.01.

According to the definition of $L$, for the evenly distributed in-plane circle formation with a radius of $r$, the monopolized sphere radius of each satellite is equal, which is approximate as follows:

$$m_i=\frac{2\pi r}{2N}=\frac{\pi r}{N},$$

(27)

where $N$ is the cluster size.

The distribution uniformity of the cluster is:

$$L=\frac{{\displaystyle \sum _1^n{V}_i}}{V}=\frac{N\times V_i}{V}=\frac{N\times 4/3\pi m_i^3}{4/3\pi {\left(r+2m_i\right)}^3}=N{\left(\frac{\pi }{N+2\pi }\right)}^3.$$

(28)

When $N=50$, the distribution uniformity of the cluster $L\approx 0.01$. The simulation results show that the control law designed in this research can control the satellite cluster to establish the desired evenly distributed in-plane circle formation.

Adjust the value of $\mu$, but keep it negative while leaving the rest of the parameters unchanged. The value of $\mu$ does not change the position of the minimum potential field but will affect the gradient of the potential field. The satellite cluster can still spontaneously establish a circle formation, but the time and energy consumption are different. The simulation control parameters were set as follows:

$$k_c=0.25 \mathrm{kg}/{\mathrm{s}}^2,k_f=3 \mathrm{kg}/\mathrm{km}{\mathrm{s}}^2,k_v=2000 \mathrm{kg}/\mathrm{s},a=b=1,r=0.5 \mathrm{km},s=r,\mu =-20 \mathrm{km}.$$

Figure 15 shows the motion trajectory of the satellite cluster to establish an in-plane circle formation in the target orbital frame when $\mu =-20 \mathrm{km}$. Compared with $\mu =-5$, the time required to establish the in-plane circle formation is reduced by about half. From the relative motion trajectory of 0–500 s, it can also be seen that the control force is too large, resulting in unnecessary fuel consumption.

By adjusting the value of $r$, the radius of circle formation can be easily changed. Figure 16 shows the motion trajectory to establish an in-plane circle formation with a radius of 5 km, and the simulation control parameters were set as follows:

$$k_c=0.25 \mathrm{kg}/{\mathrm{s}}^2,k_f=3 \mathrm{kg}/\mathrm{km}{\mathrm{s}}^2,k_v=2000 \mathrm{kg}/\mathrm{s},a=b=1,r=5 \mathrm{km},s=r,\mu =-5 \mathrm{km}.$$

3.1.3. Concentric Double Circle Formation

When $\mu >0$, the satellite cluster can establish a concentric double circle formation with the radius of $r\pm \mu$. At the initial time, satellites randomly operate in certain airspace and receive the control command to establish an in-plane concentric double circle formation with the radius of 1.5 km and 0.5 km, taking a target spacecraft as the formation center. The target spacecraft runs on a nearly circular orbit with an altitude of 600 km, an inclination of 98°, and an eccentricity of 0.0002. The simulation control parameters were set as follows:

$$k_c=0.25 \mathrm{kg}/{\mathrm{s}}^2,k_f=3 \mathrm{kg}/\mathrm{km}{\mathrm{s}}^2,k_v=2000 \mathrm{kg}/\mathrm{s},a=b=1,r=1 \mathrm{km},\mu =0.5 \mathrm{km}.$$

Figure 17 shows the motion trajectory of the satellite cluster to establish in-plane concentric double circle formation with the radius of 1.5 km and 0.5 km in the target orbital frame.

3.1.4. Disk Formation

In particular, when $a=b=1,r=0 \mathrm{km},\mu <0$, the satellite cluster can establish a disk formation. The radius of the disk formation is proportional to the inter-satellite repulsive force. Figure 18 shows the motion trajectory of satellite cluster to establish an in-plane disk formation in the target orbital frame, and the simulation control parameters were set as follows:

$$k_c=0.25 \mathrm{kg}/{\mathrm{s}}^2,k_f=3 \mathrm{kg}/\mathrm{km}{\mathrm{s}}^2,k_v=2000 \mathrm{kg}/\mathrm{s},a=b=1,r=0 \mathrm{km},s=1 \mathrm{km},\mu =-1 \mathrm{km}.$$

3.1.5. A Variety of Formations

The formation control algorithm based on bifurcating potential field was verified by simulation of several cases of satellite cluster formation. A variety of formation topologies were achieved by simply changing the bifurcation parameter and configuration adjustment parameters, as shown in Figure 19.

Furthermore, these formations can be analyzed by the five-element characterization method. When the cluster size $N=50$ and the outer boundary of the group $d_{\max }=1 \mathrm{km}$, the distribution uniformity of the cluster corresponding to various evenly distributed formations is shown in

Table 3. The distribution uniformity of the cluster corresponds to different formations.

Evenly Distributed Formation (km)	$\mathit{L}$
Ellipse (semi-major/minor axis 1/0.6)	0.005165486521583
Circle (radius 1)	0.010232661152041
Concentric double circle (radius 0.7 ± 0.3)	0.024920907802092
Ellipse disk (semi-major/minor axis 1/0.5)	0.050473947816051
Disk (radius 1)	0.071896892978728
Ellipsoidal surface (semi-axis 1/0.8/0.6)	0.202119553694005
Spherical surface (radius 1)	0.341012298531677
Ellipsoidal space (semi-axis 1/0.8/0.6)	0.482978703021702
Spherical space (radius 1)	0.765401006065638

.

The distribution uniformity of the cluster can reflect the running state of the satellite cluster to a certain extent. These data suggest that:

• The distribution uniformity of in-plane formation is lower than that of stereo formation;
• Among in-plane formations, the distribution uniformity of linear formation is lower than that of disk formation;
• The distribution uniformity of monolayer formation is lower than that of bilayer formation.

3.2. Formation Reconfiguration

A satellite cluster composed of 50 satellites with a weight of 100 kg randomly operates in certain airspace. At the initial time, the satellite cluster receives the control command to establish an in-plane disk formation with a radius of 0.5 km, taking a target spacecraft as the formation center. The target spacecraft runs on a nearly circular orbit with an altitude of 600 km, an inclination of 98°, and an eccentricity of 0.0002. With the transformation of the space mission, the satellite cluster formation is successively reconfigured to the circle formation with a radius of 0.5 km, the line formation along $z=0 \mathrm{km}$, the double line formation along $z=\pm 0.5 \mathrm{km}$, the ellipse formation with a semi-major axis of 1 km and a semi-minor axis of 0.6 km, and the concentric double circle formation with the radius of 1.5 km and 0.5 km. The simulation control parameters are shown in

Table 4. Formation reconfiguration control parameters.

Formation Reconfiguration (km)	$\mathit{a}$	$\mathit{b}$	$\mathit{r} (\mathrm{km})$	$\mathit{\mu } (\mathrm{km})$
Disk (radius 0.5)	1	1	0	−1
Circle (radius 0.5)	1	1	0.5	−5
$\mathrm{Line} (\mathrm{along} z=0$)	1	0	1	−1
$\mathrm{Double} \mathrm{line} (\mathrm{along} z=\pm 0.5$)	1	0	1	0.5
Ellipse (semi-major/minor axis 1/0.6)	1	0.6	1	−5
Concentric double circle (radius 1 ± 0.5)	1	1	1	0.5

.

Figure 20 shows the relative motion trajectory of satellite cluster formation reconfiguration based on the bifurcating potential field. The satellite cluster formation can be easily reconfigured by adjusting the bifurcation and configuration parameters. This method does not need to set the desired target state of each satellite in the cluster. After receiving the command to switch formation configuration, each satellite can spontaneously approach the target formation according to its current relative state. For example, in the process of satellite cluster switching from the ellipse formation to the concentric double circle formation, the satellites closer to the inner circle are finally evenly distributed on the inner circle, while the satellites closer to the outer circle are evenly distributed on the outer circle.

The five-element characterization method can represent the characteristics of the orbiting cluster at any time. In the process of switching from disk formation with a radius of 0.5 km to circle formation with a radius of 0.5 km, the cluster’s characteristic parameters at several selected moments are shown in

Table 5. The cluster’s characteristic parameters.

Parameters	0 s	1000 s	1050 s	1100 s	1500 s	2000 s
$N$	50	50	50	50	50	50
$\left({\mathit{r}}_o,{\mathit{v}}_o\right)/$ $(\mathrm{km},\mathrm{km}/\mathrm{s})$	(−443.54295	(5850.6303	(6037.1365	(6205.9070	(6869.2879	(5917.3955
	899.34186	1530.4312	1517.0148	1499.1419	1204.6086	533.21440
	−6898.0861	−3471.9062	−3143.3833	−2805.6007	106.34651	3654.0846
	7.4271828	3.9039044	3.5545128	3.1946381	0.069789521	−3.7823796
	1.3587345	−0.22348015	−0.31304389	−0.40170255	−1.0478601	−1.5714632
	−0.29933858)	6.4713810)	6.6663189)	6.8416710)	7.4884619)	6.3527428)
$d_{\max }/\mathrm{km}$	0.9234	0.9214	0.9661	1.004	1.003	1.003
$d_{\min }/\mathrm{km}$	0.09163	0.09440	0.05093	0.003256	0.02745	0.04681
$L$	0.06967	0.07092	0.02099	0.004614	0.006609	0.008613

.

Figure 21 shows the change of cluster’s distribution uniformity during the process of switching from disk formation to circle formation. At 0–1000 s, the satellite cluster maintains the disk formation configuration, and the distribution uniformity of clusters remains at 0.7. At 1000 s, the cluster receives the control command to transform into a circle formation. When the target formation is established and evenly distributed, the distribution uniformity of the cluster tends to the steady-state value of 0.01.

4. Conclusions

Aiming at the mathematical characterization of large-scale satellite clusters, we provide the five-element characterization method to represent the cluster characteristics and internal correlation characteristics of orbiting satellite clusters. We also propose a satellite cluster formation control method based on bifurcating potential field, which solves the problems of multiple formation configuration control and realizes the dynamic migration and rapid reconfiguration of satellite cluster formation. According to the bifurcating potential function, a variety of satellite cluster formation configurations can be formed, and each configuration has clear correspondence with the bifurcating parameters $a,b,c,r,\mu$. The feasibility and effectiveness of the proposed formation control algorithm were verified by simulation with a cluster composed of 50 satellites. The results show that multiple formation topologies can be realized by simply changing the bifurcating parameters. The proposed five-element characterization method can describe the operation state of orbiting satellite clusters.

At present, the control parameters $k_c,k_f,k_v$ are assigned according to the resultant force limit, but they are not optimal. Future work will focus on the automatic optimal selection of $k_c,k_f, k_v$ by considering the total fuel consumption, the equilibrium of fuel consumption or actual mission constraints. More quantitative analysis of fuel consumption is also needed for real satellite cluster missions.