Modeling and Control of Satellite Formations: A Survey

This survey deals with the problem of the group motion of spacecraft, which is rapidly developing and relevant for many applications, in terms of developing the onboard control algorithms to ensure the fulfillment of a given mission. The paper provides a comprehensive overview of spacecraft formation flight control. The bibliography is divided into three main sections: the multipleinput–multiple-output approach, in which the formation is treated as a single entity with multiple inputs and multiple outputs; the leader–follower formation, in which individual spacecraft controllers are linked hierarchically; and a virtual structure formation, in which spacecraft are treated as rigid bodies embedded in a common virtual rigid body. This survey expands a 2004 survey and updates it with recent results.


Introduction
Satellite formations find many applications in the various fields, including Earth surveillance, geodesy, astronomy, space exploration, etc. Following [1], the satellite formation flight (SFF) can be defined as the movement of a group of more than one spacecraft whose dynamic states are coupled by a common control law (see also [2]). At the same time, at least one satellite (formation agent) must track the desired state relative to another agent, and the tracking law must minimally depend on the state of this other agent. For example, even if specific relative positions are actively maintained, the Global Positioning System (GPS) satellites form a formation because only the position and speed of an individual satellite need to be used to correct their orbital motion.
Let us present some examples. The distributed space missions with application to Earth observation with special emphasis on radar payloads are presented in [3]. Schilling et al. [4] described the distributed satellite system Telematics Int. Mission (TIM), which is a formation of nine cooperating small satellites for photogrammetric Earth observation. TIM's application scenario focuses on the characterization of ash clouds from volcano eruptions, providing useful information to plan detour maneuvers for airplanes to avoid damage to engines. The CloudCT mission [5] uses a formation of 10 nano-satellites to detect the 3D properties of clouds. CloudCT integrates interdisciplinary synergies from nanosatellite system engineering, cloud modeling, and tomographic imaging to enable a sensor network approach to innovative Earth observation. As Freimann et al. [6] mentioned, nano-satellites offer an affordable way to realize distributed systems for various scientific and commercial applications in the area of Earth observation and global communication services. Freimann et al. [6] presented a contact plan design approach that solves the Medium Access Control (MAC) problem by predicting additive multi-node interference and scheduling interference-free links based on these predictions.
The UWE-program aimed for in-orbit demonstrations of key technologies to enable formations of cooperating distributed spacecraft at the pico-satellite level is presented in [7]. In this program framework, the CubeSat UWE-3 is used for the experiments of evaluation of real-time attitude determination and control. The spacecraft formation consists of a group of satellites performing a single mission, which is usually difficult to accomplish with a single device [8][9][10]. For example, Bik et al. [9] described a joint project between ESA and NASA called the Laser Interferometer Space Antenna (LISA), which will become the first space observatory for gravitational waves. Its main purpose is to detect and observe gravitational waves from massive black holes and binary galaxies in the frequency range 10 −4 -10 −1 Hz. LISA consists of three identical spacecraft located at a distance of 5 × 10 6 km from each other in the form of an equilateral triangle. LISA tracks the stretching of space-time caused by gravitational waves by measuring the change in distance between two formation agents using laser interferometry. Three agents will provide independent information about the polarization of the gravitational wave. One shoulder is determined by the distance between the test masses at the ends of each branch of the triangle, which are covered by satellites so that they are not affected by non-gravitational forces. The satellites counteract non-gravitational forces using micro-Newtonian thrusters to follow the trial mass path through what is called drag-free control [11]. Schilling [12] considered technologies for implementing sensor networks in orbit for joint measurements by small satellites weighing several kilograms. It is shown that through sensor data fusion and post-processing, innovative distributed methods of Earth observation can be obtained. In [12], exemplary nano-satellite missions in the field of formations are also outlined. The first demonstration of collision avoidance and orbit control for UWE-4 pico-satellites is presented in Kramer and Schilling [13]. A possible collision with a piece of space debris was avoided by UWE-4 by use of its orbit control capabilities for an evasive maneuver in early July 2020. This maneuver increased the miss distance from 822 m to more than 6000 m. The results, presented in [13] demonstrate perspectives for small spacecraft to comply with regulations on de-orbiting and capabilities for collision avoidance maneuvers. Ref. Mathavaraj and Padhi [14] used the optimal and adaptive control techniques to develop the satellite formation flying-high precision guidance.
The mentioned examples are related to the SFF missions and the formation flying guidance (FFG), which is considered as the generation of the reference trajectories used as an input for a formation agent's relative state tracking control law. The examples mentioned are related to the SFF missions and the formation flying guidance (FFG), which is considered as the generation of any reference trajectories used as an input for a formation agents' relative state tracking control law [1]. A detailed and informative review of the results on the control of spacecraft formations can also be found in the monograph [15], which presents fundamental problems and approaches related to group flights of spacecraft, provides information on the dynamics of satellite motion, considers methods of static optimization, feedback control, and filtration, the nonlinear models of the relative spacecraft motion, including the rendezvous of satellites, orbital perturbations and methods for their effect mitigating, as well as the methods for groups of spacecraft control algorithms development.
To implement the given navigation laws by individual agents of the group, control algorithms are developed-both local control algorithms for each satellite separately, and algorithms based on the receipt by individual agents of data on the state (e.g., positions and velocities) of other satellites.
This paper provides a comprehensive overview of spacecraft formation flight control (SFFC), covering design methods and stability results analysis for these coupled state control laws. Since the SFFC architecture defines a general approach to the development of a particular SFFC algorithm, it is reasonable to organize publications according to the formation architecture. Therefore, following the lines of the 2004 review on the subject [1], in the present paper, the SFFC literature is divided into three main types: Multiple-Input-Multiple-Output (MIMO), in which the formation is treated as a single entity with multiple inputs and multiple outputs; the leader-follower (LF) formation, in which individual spacecraft controllers are linked hierarchically; and a virtual framework as the virtual structure (VS) formation, in which spacecraft are treated as rigid bodies embedded in a common virtual rigid body, and the master/slave configuration. This survey is aimed to expand [1] and present modern results of recent years. When picking up papers for reviewing, the authors have largely been guided by the influence of the papers on other publications in the field. For example, works [16][17][18] by the time of the review submission received 237, 193, 156 citations, respectively, in the database Scopus. In addition, the selection of material in this work was undoubtedly influenced by the scientific interests of the authors, who apologize in advance to colleagues whose contributions were outside the main focus of this review.
The rest of the paper is organized as follows. Section 2 deals with the MIMO formation architecture. The circle formation is described in Section 3, whilst Section 4 is devoted to leader-follower architecture. Concluding remarks of Section 5 finalize the paper. In general, the references inside each Section are given in chronological order.

Multiple-Input-Multiple-Output Formation Architecture
Vassar and Sherwoodt [18] developed a feedback controller to maintain a movement of a swarm of satellites using the calculated optimal trajectory. Chemical thrusters are used for the slave satellite as the actuators. The situation is considered when two satellites fly at a distance of 700 m from the geophysical orbit. The results of [18] show that movement in the swarm can be maintained at a distance of 21 m of satellites from each other.
Ulybyshev [19] studied keeping the satellite formation in near-circular orbits. The tangential maneuvers of a group of satellites are determined by continuous relative motion equations. During the entire movement, the satellites move relative to each other, and these movements are presented in the form of a graph. In [19], (2N − 1) linearized equations, where N is the total number of satellites are used. A linear-quadratic controller has been developed to control satellite formation. An analytical solution of linear-quadratic equations for the formation of two satellites is found. The simulation results for low Earth orbits (LEO) are presented where atmospheric drag is taken into account.
Inalhan et al. [20] focused on the planning of missions that provide the necessary conditions for the initialization of a periodic formation of a satellite constellation that provides a good distributed image of the Earth, for example, using synthetic aperture radar. These passive apertures were previously designed using the closed form solution provided by the Hill equations (see [21,22]), also known as the Clohessy-Wiltshire equations [23], which are linearized with respect to the near-circular reference orbit. A study was carried out to develop apertures that are insensitive to differential noise J 2 (the second zonal harmonic J 2 characterizes the polar compression). Based on nonlinear dynamic models (see [16]), Inalhan et al. [20] used the results by [24][25][26] on the derivation and solution of homogeneous equations of relative motion for several spacecraft. They present an initialization procedure for a large constellation of satellites with an eccentric reference orbit. The main result is derived from homogeneous solutions of the linearized relative equations of satellite motion. These solutions are used to find necessary conditions on the initial states for generations of T-periodic solutions so that the satellites return to their initial relative states at the end of each orbit. This periodicity condition and the resulting initialization procedure, originally defined (in compact form) at the perigee of the reference orbit, are generalized to include initialization at any point relative to the reference orbit. In particular, in [20], an algorithm is presented that minimizes the fuel consumption associated with the initialization of the satellite state to values that are consistent with periodic relative motion. Periodicity conditions and homogeneous solutions can also be used to estimate relative motion errors and approximate costs of fuel associated with neglect of the reference orbit eccentricity. The simulation results are presented, showing that ignoring the eccentricity leads to an error close to the error from disturbances, caused by differential accelerations of gravity.
The movement of a satellite formation relative to an eccentric reference orbit is discussed in [27], where optimal fuel/time ratio algorithms for low-level maintenance of station movement in the presence of interference are given. The algorithm by Tillerson and How [27] is optimized by applying the linear programming method with respect to a time-varying dynamics model linearized with respect to an eccentric reference orbit. Numerous results of nonlinear system simulation are given, demonstrating the efficiency of this general approach to control. In [27], it is shown that even in the presence of differential interference J 2 , the proposed approach to control the movement of the formation is very efficient, which requires fuel consumption (expressed in terms of the change in velocity during the orbit) ∆V ∈ [5,15] mm/s/orbit, depending on the scenario. Simulations showed that Lawden's equations [24] are needed to determine the desired state for periodic relative motion, but the Hill equations are sufficient for posing the linear programming control problem. This result is important because the use of the time-independent Hill equations [22] significantly reduces the computational costs required to solve a linear programming problem.
Two nonlinear feedback control laws for recovering the desired J 2 -invariant relative orbit are presented in [16]. According to [28], there are two necessary conditions for two adjacent orbits to be invariant to the perturbation of the gravitational harmonic J 2 in the sense that both orbits will demonstrate the same secular angular drift velocities. In particular, these conditions ensure that the velocities of the ascending nodeΩ and the mean values of the latitudeθ M =Ṁ +ω, where M is the mean anomaly and ω is the argument of perigee, which are equal to each other. Since it is convenient to describe the relative orbit of the follower ("daughter") satellite relative to the lead ("parent") satellite in terms of the mean differences of the orbital elements and since the conditions for the relative invariance of J 2 relative to the orbit are expressed in terms of the mean orbit elements, the first control law returns errors in terms of mean orbit elements. Using the designations of the orbital elements with the index 1 for the lead satellite, the differences in the momenta of the J 2 -invariant relative orbit elements of the slave satellite are written in [16]: where a is the semi-major axis of the orbit, i is the inclination, the eccentricity index η is defined as η = √ 1 − e 2 , e is the eccentricity. After choosing a certain middle element of the difference δa, δe or δi, the remaining two differences of the momentum elements are given by two constraints in the Equations (1) and (2). Since only the pulse elements a, e and i affect the secular drift caused by J 2 , the mean angles M, ω and Ω can be chosen arbitrarily. To find the corresponding vectors of inertial position and velocity, the elements of the mean orbit are first translated into the corresponding oscillating elements of the orbit using the theory of the artificial Brouwer satellite [29]. Since the relative orbit is described in terms of the relative differences in the mean orbit elements when establishing J 2 -invariant relative orbits, Schaub et al. [16] constructed the feedback law in terms of the mean orbit elements instead of the more traditional feedback approach on position and velocity vector errors. Not all orbital position errors are the same. The ascending node error should be controlled at a different position in orbit than the pitch error. Gaussian variational equations of motion (see [30]) provide a convenient set of equations relating the influence of the control vector of accelerations u = [ u r u θ u h ] T with time derivatives of the osculating element of the orbit: d d where u r is directed radially away from the Earth, u h is directed along the orbital angular momentum vector, u θ is orthogonal to the previous directions. Parameter f is the true anomaly, r stands for the orbital radius, p = a 1 − e 2 denotes the focal parameter (semilatus rectum), latitude θ = ω + f . The average angular velocity n is determined by the expression These equations in variations were obtained for the Keplerian motion. In matrix form, they take the formė osc = 0 0 0 0 0 n T + B(e osc ) u, where e = a e i Ω M T is the vector of the osculating orbit elements, and B(e osc ) is the control influence matrix of size 6 × 3 corresponding to the right-hand sides of Equations (4)- (9). The vector of the mean orbit e = a e i Ω M T is introduced and the analytical transformation e = ξ(e osc ) from osculating elements of the orbit e osc to the mean orbit vector e. Schaub et al. [16] used the first order truncation of Brouwer analytical solution [29]. Taking into account the influence of the disturbance J 2 , the Gaussian variational equations for the averaged motion are written in the formė where (6 × 1)-matrix A(e) of the controlled plant corresponds to the right-hand sides of (4)- (9). Considering the Brouwer transformation between the oscillating and mean elements of the orbit, one can see that matrix ∂ξ ∂e osc is approximately 6 × 6-identity matrix with off-diagonal elements of the order of J 2 or less. Therefore, in [16], the equation for the average velocity of an orbital element is approximated in the forṁ The plant model matrix [A] accurately describes the behavior of the elements of the mean orbit. Perturbation J 2 has no secular impact on elements a, e, i. The control influence matrix [B], obtained in the Gaussian variational equations, as seen from (4)-(9), makes it possible to calculate the change in the osculating orbital elements due to control vector u. It is assumed that changes in these vibrational elements of the orbit, as follows from (12), are directly reflected in the corresponding changes in the elements of the mean orbit. Working with the elements of the mean orbit has the advantage that short-term period fluctuations are not perceived as tracking errors; only long-term errors are compensated.
Using the quadratic Lyapunov function V(δe) = 1 2 δe T δe, where δe = e 2 − e d is the tracking error of the slave satellite (e 2 ) orbit of the lead satellite (e d ) in terms of mean orbit elements, in [16] the following control law for the slave satellite was obtained where the time-variable matrix of parameters [P] is chosen when synthesizing the algorithm. Ref. Schaub et al. [16] recommended one way of doing this. The second control law, proposed in [16], uses the position and velocity errors in Cartesian coordinates r = (x, y, z). Based on the mismatch δr = r 2 − r 2 d between the current r 2 and the desired r 2 d the position of the slave satellite, based on the Lyapunov method, the control law is obtained where function f(r) is the right-hand side of the equation of free motion of a satellite in orbit,r = f(r).
In [16], comparative results of a numerical study of the proposed control laws are presented, and further in [31], the main method for constructing orbits with different Cartesian coordinates is presented. The numerical study illustrates the accuracy with which the coordinate transformation occurs after linearization. A hybrid law with continuous feedback is proposed, in which the desired orbit geometry is explicitly given in the form of a control error, and the actual orbit is given in local Cartesian coordinates. Numerical simulations illustrate the efficiency and limitations of such feedback control laws. Using a linearized mapping between relative orbital coordinates results in a slight degradation in quality. However, this method gives the advantage of working in the main space of elements in the case when the relative errors in the construction of the orbit are determined.
One of the fundamental differences between satellite movement and conventional rendezvous operations, as noted in [32], is the need to increase the duration of the movement. Therefore, if the formation is to be maintained, then long-term disturbances, in particular those caused by the Earth's compression, must be corrected. Williams and Wang [32] explored the use of motionless aids to counteract the differential precession of the orbital plane, which is the primary perturbation affecting the formation caused by the Earth's compression. The approach used is to use the pressure of solar radiation acting on a relatively small surface called the solar wing, which is attached to the satellite. The resulting torque causes the orbit to precess; if the wing is the correct size, this motion will, on average, suppress the motion caused by flattening, thereby maintaining the formation without using propellant. Ref. Williams and Wang [32] considered satellites that are nominally Earth-oriented; most small satellites that are stabilized by the gravitational field gradient or have sensors designed to observe the Earth fall into this category. It is assumed that a relatively small reflective surface is attached to the spacecraft at an angle α (usually 45 • ) to the velocity vector. The reflector surface illuminated by the Sun will alternate as the satellite moves around its orbit, creating an out-of-plane component of the alternating solar radiation force. Consequently, the torque that is created in orbit by this force will tend to increase rather than decrease during the entire revolution, leading to a total torque along the direction to the Sun. If the wing area is appropriately sized, this torque can be used to zero out the differential nodal regression velocity between the two orbits, thus maintaining the formation without a propulsion system. Williams and Wang [32] demonstrated that the long-term orbital effects of the solar wing's steering input (wing orientation angle) are highly non-linear and show a strong relationship between orbital inclination and ascending node longitude, and if the wing is properly sized, node drift due to solar torque, on average, is compensated by the differential speed of modality, so the formation is preserved without the use of propellant. Numerical results are presented to illustrate the approach.
Hussein et al. [33] dealt with establishing relationships with the motion of a satellite group in a physical two-dimensional space with visualization problems. A class of movements is presented that satisfies the visualization tasks, and simple solutions to the problem of motion design and formation control for visualization applications are found, allowing to achieve a compromise between image quality and fuel consumption for the proposed class of maneuvers.
As part of a study by the European Space Agency (ESA), Gill and Runge [34] discussed the potential advantages, disadvantages, and challenges of moving satellites in a formation in close proximity to implement synthetic aperture radar (SAR) for trajectory interferometry. Gill and Runge [34] noted that the technology of using SAR is well developed. Radar signals reflected from the Earth's surface, water, and ice have provided high-resolution images for areas such as forestry, geology, ocean wave monitoring, and ice mapping (see [35]). In addition, progress in the processing of SAR data has been made through the use of phase differences between SAR images of the same scene, a technique referred to as interferometric SAR (InSAR). Distinguishing SAR images from two laterally-spaced antennas generally yields terrain elevation measurements. In contrast, along-track interferometry (ATI) is not based on a specific interferometric baseline, but on the time interval between two moving SAR devices separated in the direction along the track, and displaying the same terrain, which provides measurements of the speed of objects on the surface. Space ICT was successfully demonstrated using the Shuttle Radar Topography Mission (SRTM) in 2000 [36,37]. A serious limitation of single-satellite IWT is due to the physical dimensions of the spacecraft. It can be overcome on the basis of the distributed sensor concept, when two SAR antennas are located on different platforms, thus forming variable baseline. This can be conducted with a small satellite that is passive in terms of radar, but active in terms of orbital maneuvers that are conducted relative to the larger active SAR satellite. Alternatively, two similar SAR satellites that are capable of moving in a controlled formation can be used.
In IVT interferometry, two antennas are used, separated in the direction of motion along the spacecraft trajectory by a certain amount ∆L. Considering the common orientation of the antenna, both antennas thus display the same scene with a time lag of ∆t = ∆L/v sc , where v sc denotes the spacecraft speed. When one antenna works for both transmit and receive, and the other only for receive, then the two subsequent sets of SAR phase measurements differ by the value where v r is the projection of the target speed onto the direction of the line of sight, λ is the radar radiation wavelength. Thus, the measurement accuracy can be increased either by increasing the base ∆L, or by decreasing the wavelength λ. It should be noted that large values of the base cause the effect of phase wrapping and associated ambiguities, if the velocity range of interest is mapped to a phase interval of more than 2π. Based on the constraints imposed by the decorrelation of the SAR phase in oceanographic applications for extended baselines, in [34] the upper limit of the relative baseline along the path of the spacecraft formation is 50 ms −1 GHz / f v sc , which is about 300 m in the L-band and 40 m in the X-band. (Recall that L-band: 1.0-2.0 GHz, X-band: 8.0-12.0 GHz.) To provide a margin of accuracy with respect to the decorrelation threshold, the appropriate separation of spacecraft in the formation along the trajectory would be 225 ± 75 m in the L-band and 30 ± 10 m in the X-band. There are two main options for controlling a formation: ground-in-a-loop control and autonomous formation control. Obviously, in the second case, information is required on board the spacecraft about their relative position and speed. To the relative position of the two spacecraft prediction, it is necessary to integrate the equations of their motion taking into account all relevant perturbations, including the effects of a complex form of the Earth's gravitational field, the influence of the Sun, and the Moon, atmospheric drag and solar radiation pressure [38]. For the motion of two spacecraft in the limiting case of near-circular orbits and close formations (with a distance of less than 1 km), closed forms of the Clohessy-Wiltshire equations of relative motion [23] can be obtained.
A number of works are investigating the possibility of using the natural electrostatic forces of spacecraft in High Earth Orbit (HEO), or in deep space [39]. Investigation of [40], showed that the absolute charge of spacecraft can reach the level of kilovolts in the Geostationary Earth Orbit, GEO. If spacecraft fly at a distance of several tens of meters from each other, this causes disturbances ranging from micro-to millinewtons. In an orbit, these disturbances can lead to hundreds of meters of deflection. The concept of Coulomb thrust proposes to use active control of the spacecraft charge to ensure its desired values and use this disturbing force to directly control the relative motion KLA [41,42]. The electrostatic charge of the spacecraft is partially hidden from another nearby ship due to interaction with ions and electrons of cosmic plasma. The strength of this shielding is measured by the Debye shielding length [43]. The cold and dense plasma of the environment in low-orbit orbits (LEO) has a Debye length of about a centimeter, which makes it impossible for Coulomb interaction at a distance of tens of meters. In a geostationary orbit, the Debye length is in the range of 100-1000 m, which makes it possible to use the Coulomb thrust between satellites. In deep space, at a distance of 1 AU from the Sun, the Debye length is 30-50 m [40]. Coulomb thrust usually requires only a watt level of electricity, while consuming essentially no fuel. The specific impulse I sp of the Coulomb thrust is in the range 10 9 -10 12 width [40,44]. Spacecraft charge control has been demonstrated on the SCATHA [45], ATS [46], Cluster [47,48] missions. Coulomb structures have electrostatic forces that fully compensate for differential gravitational acceleration. As a result, a cluster of spacecraft is formed, the positions of individual vehicles which seem to be frozen relative to the rotating main locally-vertically-locally-horizontal system [39,40,[49][50][51]. However, all static equilibrium solutions of relative position in orbit or deep space are unstable, and active recharging feedback is required to stabilize them.
Due to the distributed nature of spacecraft formations, methods of decentralized filtration [52] are of great interest. In [53][54][55], the estimation is performed using the parallel operation of full-order observers with local dimensions. For example, [54] considers the dynamics of N agents given by the linear difference state equations where x[k] ∈ R n is the system state vector, u[k] ∈ R mN is vector control, y i [k] are the measured outputs of agents, n i [k]-measurement errors, i = 1, . . . , N. It is assumed that each vehicle has the same number m of executive devices, and also that each agent with the number i has a local controller that generates the corresponding control signal is used such that the matrix of equations of state of the entire closed-loop system A clp = A + BK is asymptotically stable (its spectrum lies inside a circle of unit radius). This leads to N dynamic controllers according to the state estimates are generated by the observerŝ where L is the observer's feedback matrix to be chosen in the synthesis. Each control signal u i [k] is generated based on a local estimate of the state of the entire systemx i [k]. This approach uses an estimate of the actuator input applied by each of the other vehicles and will generally not be correct for non-local components. The incompleteness of information about the input data of the object leads to the fact that the dynamics of the estimation error become bounded. Further, [54] uses the well-known method of communication specification based on graph theory [56,57]. It was found that communication between vehicle observers can be used to eliminate or mitigate the consequences of the interconnected dynamics of estimation errors. The communication structure and analysis presented in [54] provide a means of designing a communication network that can determine the dynamics of interconnected errors. It is shown that in this structure the minimum number of communication channels required to eliminate these dynamics grows only linearly with the size of the formation. Another result is that in terms of network robustness and resiliency, "receivers are much more valuable than transmitters". A similar problem is considered in [58], where the exchange of navigation data between two satellites is considered, aimed at reducing the load of the inter-satellite communication channel, taking into account the dynamics of the relative motion of the satellites and possible erasures in the navigation data of the channel. For this, an adaptive binary coding/decoding procedure is used to transmit navigation data between satellites and a closed-loop control law that regulates the relative motion of the satellites.
Using the calculus of variations, Lawden [24] deduced the necessary conditions for the optimal fuel consumption problem of multi-impulse rendezvous, known as the primer vector theory (see also [59]). In [60], a numerical method using Lawden's necessary conditions is proposed. Lawden's theory is applied in [61] to the task of rendezvous of the Apollo spacecraft. Research in this direction was developed in [62][63][64][65][66]. Prussing et al. [63] used basis vector theory to find optimal time-fixed multi-impulse coplanar and bounded class non-coplanar rendezvous solutions. In [63], the spacecraft motion equations in the gravitational field are used through the orbital radius vector r [67]: where Γ is the thrust force modulus expressed in acceleration units (0 Γ Γ max ), u is the unit vector in the direction of thrust action, J is the characteristic velocity to be minimized. Equations (18)- (20) describe the behavior of the state vector x T = r T v T J under the action of free-fall acceleration g(r) and control variables Γ and u. For a high-thrust engine, one can make an approximation of impulsive thrust, assuming that the magnitude of thrust is not limited (Γ max = ∞). In this case, the engine is either turned off (Γ = 0), or provides an impulse thrust of infinitely short duration. The thrust impulse solution requires the determination of the times, locations, and directions of thrust impulses that meet the specified boundary conditions for orbit transfer, interception, or rendezvous. Determining the solution with the minimum amount of fuel requires solving the optimal control problem in the time interval t ∈ [t 0 , t F ], which minimizes the final value of J and satisfies the equations of motion and the orbital boundary conditions of the problem. The Hamiltonian following from (18)-(20) must be maximized taking into account the equations for the conjugate variableṡ where G(r) is the symmetric gravity gradient matrix. Since the characteristic velocity J is not limited, λ J (t) = −1 is satisfied. The boundary values of the remaining conjugate variables depend on the terminal constraints r(t F ) and v(t F ) on state variables. The Hamiltonian H is maximized by choosing the direction of the thrust force that maximizes the scalar product λ T v u, that is, combining the thrust vector with the direction conjugate to the speed, which is the basis vector (see [24]), denoted below by p. Then, from (21)-(24) one can deduce thaẗ where p denotes the amplitude of the basis vector. From (26) for the Hamiltonian it follows that the thrust switch function is p − 1. In the case of continuous thrust, the Hamiltonian is maximized by choosing Γ = 0 when p < 1, and Γ = Γ max when p > 1. In the impulse case Γ = 0, when p < 1 with impulses arising at those times when p(t) touches p = 1 from below [24] Further, in [63] multi-impulse solutions are obtained for both coplanar and a limited class of non-coplanar optimal fixed in time rendezvous "circular orbit/circular orbit". For sufficiently long transition times, solutions become well-known time-open solutions, such as the Hohmann transfer for the coplanar case. (In orbital mechanics, the Hohmann transfer orbit is an elliptical orbit used to transfer between two circular orbits of different radii around the central body in the same plane. An orbital maneuver uses two engine impulses to execute it: one to put the spacecraft into a transfer orbit, and the second to lift it from it.) The results of [63] can be used to find a trade-off between time and fuel consumption for missions that have time constraints, such as rescue missions or evasive maneuvers in outer space.
The studies were continued in [64,66], where impulsive and time-fixed solutions with minimum fuel consumption for the problem of rendezvous and interception in orbit with internal trajectory constraints were obtained. Transitions between coplanar circular orbits in a gravitational field along a circular trajectory representing the minimum or maximum radius of the admissible orbit are considered. Basis vector theory has been extended by introducing trajectory constraints. In [64,66], Equations (18) and (19) are used, supplemented by the objective function to be minimized J = T 0 Γ(t)dt corresponding fuel consumption. For the n-pulse solution, the thrust force (per unit mass) is described by the expression Γ = ∑ n k=1 ∆v k δ(t − t k ), where 0 t 1 . . . t n T and ∆v k δ(t − t k ) represent impulses at times t k , expressed in terms of the δ denotes the Dirac function. Therefore, the objective function takes the form J = ∑ n k=1 ∆υ k , where ∆υ k are the discontinuities of the amplitudes of the velocity vector. In [64,66], a restriction of the inner minimum or maximum radius of the inner trajectory is introduced in the form r(t) r min or r(t) r max and is implemented by forcibly limiting the radius of the periapsis or apocenter. Initial orbit and restrictions on the place of meeting, or interception, are defined in terms of the constraint function Ψ r(0), v(0), r(T), v(T) = 0. Otherwise, without introducing restrictions, the problem statement is the same as in [63]. In [64,66] the optimal number of pulses is determined, as well as their time and moments of appearance. The existence of bounding boundary arcs is investigated, as well as the optimality of the class of singular arc solutions. In [66], a bifurcation phenomenon is discovered that can occur in an optimal solution for a constrained encounter problem. It consists of the fact that for symmetric boundary conditions, symmetric optimal solutions are expected and they found for sufficiently short transit times. However, when the transition time increases to some critical value, some of the extreme solutions become asymmetric and ambiguous, and the only symmetric extreme solution splits into two asymmetric ones.
For the problem of controlling a formation of two satellites (according to the masterslave scheme), [68] (see also the review [69]) considers various methods for applying continuous control in sliding modes, including second-order sliding mode algorithms [70], for example, super-twisting algorithm, continuous sliding mode control with a sliding mode disturbance observer (SMDO) [71], algorithms with integral sliding surface (ISS) [72]. Pulse Width Modulation (PWM) is used to convert continuous control signals into pulse trains that can be implemented using satellite jet engines. A numerical example is considered in which the position of the slave satellite relative to the master is described on the appropriate time scale by the following equations where x, y, z are the coordinates of the slave satellite relative to the master u x , u y , u z are the control forces acting on the slave satellite, d x , d y , d z includes various perturbations acting on a two-satellite system, caused, for example, by the non-centrality of the gravitational field, solar pressure, lunar-solar tides, the rotational-vibrational motion of the Earth in space, oscillations of the Earth's pole, uneven rotation of the Earth, and also the deviation linearization error position. The desired relative position of the satellites is given by the variables x c , y c , and z c . Appropriate sliding variables are formed, including an integral term to compensate for constant perturbations (see [71] (pp. 147-153)): (28) with some parameters of the desired motion in the sliding mode For numerical analysis, the values c i,0 = 1.8, c i,−1 = 1.0 are taken. The following formation control algorithms are considered: (1) Sliding mode controller and sliding mode perturbation observer Here,v i eq is the "equivalent control", the result of passing signals where ρ i > 0, through low-pass filters. (According to the authors of the review,v i eq is not exactly equivalent control, as defined in [71], since, unlike the indicated works, in (29) the filtering occurs in a closed rather than open control loop.) As noted in [68], thev i eq signals, after a finite transient time, are an accurate estimate of the matched reduced perturbations ψ 0 i (·) (linear combinations of object output and e i tracking error (see [68,72] for details)).
(2) The following second order sliding modes super-twisting algorithm, ensuring convergence of variables on the sliding surface and their derivatives to zero in a finite time (3) The sliding modes continuous algorithm with an integral sliding surface, based on the Lyapunov method. This algorithm implements continuous control in a sliding mode, ensuring robust stabilization of σ i with respect to matched perturbations. It uses two types of moving variables: σ i and additional integral moving variables η i [72]. It can be written as whereρ i > 0,k 0i > 0, 0.5 < a i < 1.0, are the constants, chosen by a designer;ũ = E(x)u, where E(x) denotes the nonsingular transformation matrix to a new basis in which the original system is reduced to m independent subsystems. Whenũ is found, the control u is obtained by the inverse transformation u = E −1 (x)ũ. For the considered problem, (27), algorithm (31) and (32) read as where a i = 0.9,k 0i = 0.5,ρ i = 10, i ∈ {x, y, z}. The control signals u i in (29)-(33) are continuous, as discontinuous high frequency components are filtered or integrated. The comparative results of applying the proposed control algorithms are presented in [68]. Kim et al. [10] explored the three types of reconfiguration maneuvers described in [73] such as resizing, reassigning, and reorienting. Resizing is the change in distance between spacecraft in the formation. In [10], a hybrid optimization strategy with a genetic algorithm and basis vector theory is used for fuel-optimal reconfiguration of a satellite group. In hybrid optimization, the genetic algorithm performs a global search to create two-impulse trajectories at each transition time, and vector analysis of basis vectors serves for local optimization with two-impulse trajectories as initial approximations. The final optimal trajectory has multiple pulses and also lower fuel consumption than two-pulse trajectories. The mission planner must investigate the trade-off between fuel saved by multiple impulses and the increased mission complexity due to additional impulses. In many missions, the trajectory obtained by the genetic algorithm will have an advantage for the transition due to its simplicity and almost optimal fuel consumption.
In a series of papers, Udwadia and Kalaba [74,75], Udwadia [76], Udwadia and Kalaba [77], Udwadia [78,79,80,81], Udwadia and Kalaba [82] the approach is proposed based on the Gauss principle to obtaining models of dynamical systems with constraints. In particular, Udwadia and Kalaba [74] noted: "The principles of analytical mechanics laid down by D'Alembert, Lagrange (1787), and Gauss (1829) are all-encompassing, and therefore it naturally follows that there can be no new fundamental principles of the theory of motion and equilibrium of discrete [finite-dimensional] dynamical systems. Despite this, additional perspectives can be helpful in understanding the laws of Nature from new perspectives, in particular, if they can help in solving problems of particular importance and in providing a deeper understanding of how Nature works. Although the general problem of limited motion was formulated at least as early as the time of Lagrange, the definition of explicit equations of motion for discrete dynamical systems with constraints, even within the limited view of Lagrange mechanics, was a serious obstacle. The Lagrange multiplier method is based on problem-oriented approaches to determining the multipliers; they are often very difficult to find and, therefore, to obtain explicit equations of motion (both analytically and numerically) for systems that have a large number of degrees of freedom and many non-integrable constraints". In [74], a finite-dimensional dynamical system whose configuration is described by n generalized coordinates q = [q 1 , q 2 , . . . , q n ] T is considered. Its motion can be described, using the Newtonian or Lagrangian formalism, by equations where the n × n matrix M is symmetric and positive definite. Generalized accelerations a of the system without restrictions, therefore, are given by the expression Now, let the system be subject to m consistent constraints (not necessarily linearly independent) A(q,q, t)q = b(q,q, t), where A is m× matrix, commonly called the "constraint matrix", b is the m -dimensional vector. Differentiation of the usual constraint equations used in Lagrangian mechanics, which often have a Pfaffian form, leads to equations of the form (36). These constraint equations therefore include, among others, conventional holonomic, nonholonomic, scleronomic, rheonomic types of restrictions and their combinations. Thus, Equation (35) describe an unbounded system, and the Equation (36) describe the connections superimposed on this system, which covers all Lagrangian mechanics. Constraints (36) are more general than those in the usual framework of the Lagrangian mechanics. The presence of constraints (36) sets additional "generalized constraint forces" on the system so that the explicit equations of motion for the constrained system take the form where the additional term Q c on the right-hand side arises due to the constraint imposed according to (36). The main result for a constrained system is formulated in [74] in three equivalent forms.
1. Explicit equations of motion that describe the evolution of a system with constraints: or where the matrix K(q,q, t) = M 2. An additional term on the right-hand side of the Equation (37), representing the generalized force of constraints, is explicitly given by the formula 3. Equation (38) after multiplication by M −1 can be rewritten as or where vector ∆a =q − a represents the deviation (at time t) of the bounded generalized accelerationq from the corresponding unbounded acceleration a; the error vector e = b − Aa represents the degree to which accelerations at time t corresponding to unbounded motion do not satisfy the constraint equations (36); matrix (Matrix K 1 is called the weighted pseudoinverse Moore-Penrose matrix, [74].) The last form of [74] results in the following fundamental principle of Lagrangian mechanics: "The motion of a discrete dynamical system subject to constraints evolves at every moment of time. time in such a way that the deviation of its accelerations from those that it would have at this moment in the absence of restrictions is directly proportional to the degree to which accelerations corresponding to its unbounded motion at time t do not satisfy restrictions; the proportionality matrix is the weighted generalized inverse Moore-Penrose matrix for the weighted constraint matrix A, and the dissatisfaction measure constraints are provided by the vector e." The Udwadia-Kalaba equation, in contrast to the Lagrange equation, is equally easily applicable to both holonomic and nonholonomic constraints.
In space systems with a large number of satellites, orbit coordination between satellites is required throughout the life of a mission. The main focus has usually been on the limited relative motion between satellites in a constellation, Zhang and Gurfil [83] also considered the additional degree of freedom to manipulate an arbitrary number of orbital elements, which is represented as coordinating the overall orbital transfer and constellation in space. The basic concept uses consensus theory to characterize the properties of a control object, as in a multi-agent system. For this purpose, Zhang and Gurfil [83] assumed that communication in the network satellite system is represented in the form of an undirected graph, and then the dynamics of the control system are realized in an affine-controlled form, as described in the variational Gauss equations. For the general problem of orbital transmission, a controller based on the edge error has been developed and proven to be asymptotically stable. Several strategies for managing the grouping are discussed, namely by changing variables or a two-phase control process based on a dynamic structure. Numerical simulations are performed to validate the analysis and demonstrate the results.
Monakhova and Ivanov [84] addressed the problem of building a swarm of nanosatellites immediately after they are separated from the launch vehicle. The decentralized control is proposed using the aerodynamic drag force to eliminate the relative drift between satellites in a swarm. The effect of dividing a swarm into several independent groups is studied, which is considered a violation of the integrity of the swarm and is undesirable. An aerodynamic force application is also studied in [85], where for the formation of two satellites the control algorithms based on the modal approach (pole placement technique), on passification and variable structure control, linear and transient-speed optimal for partial stabilization algorithms are developed and examined. The stability of control systems to the plant model parameters is studied and the possibility of the occurrence of stable oscillations as a consequence of the control action (aerodynamic resistance) limitation is shown.
The problem of deploying a swarm of nanosatellites immediately after their separation from the launcher was studied in [86]. Since some error in ejection velocity during the launch is inevitable, nanosatellites have different orbital periods, and as a result, gradually disperse in orbit. The decentralized differential control using the force of resistance is applied to the swarm formation problem. It is assumed that each satellite has information about the position of others in a particular communication area. The proposed control algorithm eliminates the relative drift between adjacent satellites. Ref. Ivanov et al. [86] studied the separation effect that occurs when a swarm is divided into several independent groups. This effect depends on the size of the communication area, the number of communication satellites, and the initial conditions. The boundary values of these parameters are studied for the example of twenty 3U cubesatellites. The effect of J 2 harmonics and uncertainty on the atmosphere's density is also examined in [86].
Knowledge of position and velocity is a practical requirement for all spacecraft in low Earth orbit. Some satellite systems require autonomous orbit determination on board. The methods for estimating the absolute positions of two satellites using relative distance and azimuth are developed by Zhou and Li [87]. The spherical coordinates measurements are converted to Cartesian ones with the converted offset and covariance calculated. The transformed measurements use the Unbiased Transformed Kalman Measurement (UTKM) [88,89] algorithm to estimate the position and velocity of the satellites. A square root filter is introduced to avoid losing the positive definiteness of the state covariance matrix. The simulation results show that the algorithm proposed in [87] converges faster and has a higher estimation accuracy.
Shouman et al. [90] introduced a new control action for the formation of satellites moving in LEO, depending on the difference between the force of aerodynamic drag and the force generated by the propeller. A parameterized algorithm for regulating the output signal for the formation of flight missions has been developed, based on the equations of relative dynamics by Schweighart-Sedwick [91,92]. It is implemented in order to accurately track different trajectories of the reference relative motion and to exclude the influence of disturbances from the second harmonic of the Earth's gravitational potential.
Invention [93] relates to satellite group motion control, in which the average angular velocity of all artificial earth satellites (AES) in the group is maintained equal to the average angular velocity of the passive satellite per turn. The latter is located in the central orbit of the group. Active satellites maintain their orbital position relative to passive satellites through periodic reactive correction. The technical result of the invention is to provide a given configuration of the AES system, observed from certain places on the Earth's surface.

Circle Formation Architecture
An approach to communicating satellites in a ring is also discussed by Tragesser [94]. To maintain the shape of the formation and the integrity of the bonds, the system rotates along an axis perpendicular to the plane of the ring. In [94], the influence of the Earth on this system is considered, therefore the canonical Likins-Pringle Relation (LPR) [95,96] is used to construct the orbit for maintaining the rotation of the group around the axis in the strict direction to the Earth. The equations of motion are constructed for an annular formation with an arbitrary number of satellites. Simulation is performed for the formation of three satellites and the assessment of the system stability is provided.
The first virtual structures with feedback and charge stabilization were proposed within the framework of the Coulomb cable concept with two devices [97][98][99]. While the physical tether must always be taut, the Coulomb tether can exert forces of attraction and repulsion between two ships. However, the concept of the Coulomb cable is viable only for relatively short distances up to 100 m, while the size of a typical space with a leash length of kilometers is considered within the framework of other concepts. These Coulomb mission concepts consider static scenarios in which ships are in nominally fixed locations. The first passively stable virtual Coulomb structure is a rotating system of two devices [51]. Without the screening effect of plasma charges, the force of attraction between two oppositely charged bodies is mathematically equivalent to the gravitational force of two bodies. The resulting trajectories of the two vehicles are also conical solutions, which are orbitally stable. Taking into account the finite Debye lengths and charge confinement, the relative trajectories are no longer closed conical solutions. In [51], it is shown that if ordinary circular relative motion has a diameter less than the Debye length, then the resulting non-linear motion remains stable.
In [100], the generality of the problems of three bodies rotation under the action of gravitational and electrostatic forces is investigated for finding invariant forms. A schematic illustration of the structure discussed in [100] is shown in Figure 1. Hussein and Schaub [39] examined the general relative equilibria of such a system, and also discusses the stability of an open-loop system. It is assumed that the charged spacecraft moves in circular heliocentric orbits, far from the gravitational potential fields of planets or other celestial, bodies, therefore Hussein and Schaub [39] did not take into account the relative gravitational forces and focuses on the motion of freely moving bodies under the influence of electrostatic forces. In [39], general conditions are obtained, the solutions of which are relative equilibria for a rotating group of three spacecraft with Coulomb cables. Using the methods of linear control theory, stabilizing feedback is obtained in the form of PD and PID controllers for a nonlinear system, which guarantees its convergence to the neighborhood of the desired relative equilibrium. It is noted that for asymptotic stabilization of the desired relative equilibrium, independent of the value of the angular momentum h, it will be necessary to change the nominal equilibrium value of the charge in order to reflect the actual value of the angular momentum of the system. It is also indicated that the use of general methods of nonlinear stabilization can provide larger regions of stability than the methods based on linearization considered in [39]. The decentralized state estimation algorithms and their application to the spacecraft group are discussed in [101], where three main architectures are considered: Centralized-only one point is involved, which performs primary calculations for a group of devices based on data collected at remote sites.
Decentralized-each apparatus in the group is involved, making an equal contribution to the assessment process. Each unit performs its own calculations based on the data obtained either as measurements or on the basis of the data transmitted from the entire group. In a decentralized architecture, there is no central entity to unify results.
Hierarchical-includes hybrids of the above two architectures. The large spacecraft group can be divided into smaller ones, each with its own architecture and grading system.
The simulation results given in [101] show that the decentralized reduced-order filter leads to simultaneous close to optimal estimates, as well as to the balance of communication and computational resources between the spacecraft. In addition, Ferguson and How [101] presented a hierarchical architecture for embedding decentralized assessment systems when scaling a problem for a large number of spacecraft in a group, as well as a necessary condition for a communication topology that guarantees the stability of the operation of simultaneously parallel computers and controllers.
The purpose of the research presented in [102] is an analysis of the dynamics of orientation and control of a connected satellite formation moving where the connected nodes are modeled as extended rigid bodies. Ref. [102] uses a three-in-a-row structure, and the general formula for the equations of motion of the system is obtained using the Lagrange method. The motions of the associated satellite system are analyzed in a threedimensional coordinate system in outer space. The System Dependent Riccati Equation Regulator (SDRE) is used to eliminate orientation errors. The stability region for an SDREdriven coupled satellite system is estimated numerically to show the global asymptotic robustness of the steering method. Centralized and decentralized approaches are applied to a dynamic system to compare targeting performance. It was found that the SDRE works well with both centralized and decentralized approaches to control the orientation of linked satellites in the formation of motion. The design of a formation flying orbit control system for LEO satellites based on the SDRE control approach may be also found in [103].
In [104], the problem of controlling the motion of an electromagnetic formation using the finite time method is investigated. The principle of electromagnetic formation flight (EMFF) uses magnetic fields electrically generated by all satellites, which then allows you to control the relative degrees of freedom. In EMFF, each satellite is equipped with three orthogonal magnetic coils and three orthogonal flywheels [105].
Zeng and Hu [104] presented an electromagnetic force model and analyzes the effect of the Earth's magnetic field on EMFF satellites. Then, the equations of relative motion and the method for describing the general formation are established. The robust sliding mode controller is designed to track the trajectory in the presence of model uncertainties and external noise. The proposed controller, which combines the advantages of linear and terminal sliding mode control, can guarantee the convergence of tracking errors in a finite time, and not in an asymptotic sense. By constructing a specific Lyapunov function, the closed-loop system is proved to be globally stable and convergent. Numerical simulations of the maintenance and formation reconfiguration are then presented to show the efficiency of the designed controller.
Optimal reconfigurations of the Coulomb formation of two spacecraft are found in [106] using nonlinear optimal control methods. The aim of reconfigurations is to maneuver two spacecraft between two charged equilibrium configurations. Four optimality criteria are considered: minimum time, minimum separation acceleration, minimum use of fuel for an electric motor, and minimum power consumption. Reconfiguration between equilibria is first considered by changing the desired dilution distance. In the radial configuration of the relative equilibrium, only the Coulomb force is required to control the motion in the plane and to control the satellites from their initial to final radial positions. In this reconfiguration maneuver, the torque of the gravity gradient is used to stabilize motion in the plane. For equilibrium positions along the path and normal orbits, the reconfiguration maneuver requires hybrid control. Here, the Coulomb force is changed to control the separation distance, and the inertial micromotors are activated to control the lateral directions. Second, a reconfiguration involving hybrid steering is used to maneuver ships from any initial equilibrium position to a final one. The goal is to determine optimal maneuvers that maximize the use of the Coulomb propulsion unit while minimizing the use of the electric propulsion unit. The formulation of the two-point boundary value optimization problem is numerically solved using pseudospectral methods. The Pontryagin minimum principle checks the optimality of solutions without feedback.
In [107], a quaternion-based approach is considered for solving the time-limited problem of synchronizing and stabilizing satellites while moving in a formation. Sufficient conditions are given for the finite temporal convergence and stability of the distributed consensus problem. The nonlinear control law based on the finite time control method is designed in such a way that the exact position of the spacecraft will be consistent and will converge to the leader's position, while the angular velocity will converge to zero in a finite time. A term in the form of an integral of a power function has been added to the Lyapunov function. In order to reduce the load on the network, a modified control law is proposed, where an estimate of the sliding mode with time constraints is introduced. This aims to ensure that only one satellite communicates with the leader. The modeling results are presented, showing the efficiency of the proposed method and its potential advantages.
Nair et al. [108] dealt with the development of a formation control strategy for a circular satellite constellation formation. In [108], the artificial potential field method is used to plan the trajectory, and the sliding mode control method is used to create a robust controller. The fuzzy inference engine is used to reduce the chattering phenomena inherent in conventional sliding mode systems. An algorithm for adaptive tuning of tuning a fuzzy parameter is derived on the basis of the Lyapunov stability theory. The proposed method based on a fuzzy system with a sliding mode is intended to compensate for the model uncertainties in practical applications. The simulation results for a group of five satellites forming a circular formation confirm the stability and robustness of the proposed scheme.
The purpose of [109] is to investigate the control of an electromagnetic tethered satellite system using Model Predictive Control (MPC). An electromagnetic tethered satellite system is powered by electromagnetic coils to generate control forces. The dynamic system model is described with high and low levels of accuracy, which are used to design the control system. Several horizons of the predictive model are used to bring the formation to the desired state. The presented control law not only satisfies the input and output constraints but also has the corresponding characteristics of optimality. The main advantage of using multiple horizons is multiple control based on predictive models, having a lower computational load compared to classical MPC. The results of numerical simulation and their comparison with sliding mode control are presented to demonstrate the effectiveness of the proposed control method and its advantages over both classical MPC control methods and sliding mode. The results obtained show a dramatic reduction in computational time and energy consumption compared to the classical MPC control methods and the sliding model, respectively.

Leader-Follower Formation Architecture
Fulfillment of the laws of linear control for the satellites formation in the presence of gravitational disturbances is considered by Sparks [110]. The use of a linear quadratic controller that minimizes the error between the actual and the required relative motion of the satellite is evaluated for the ability to maintain a specific geometry of the formation in the presence of such a gravitational disturbance as the uneven gravity of the Earth. The required formation geometry is based on solving linear equations of relative motion without disturbances. In particular, a group of satellites is selected whose projected motion onto the tangential plane of the Earth is a circle with a radius of one kilometer. It is shown that linear control laws support the formation within the error in the presence of gravitational disturbances. In addition, the simulation provides an estimate that takes into account the required fuel consumption to maneuver such a formation, thereby providing studies based on reducing costs while moving along the desired trajectories and according to the selected control strategies.
In [111], the formation control problem is considered in the following general setting. A complex system consisting of many subsystems (mobile agents) is described by the following equationsẋ where k is the number of subsystems; x i ∈ R n i is the state vector of the i-th subsystem; r i express interrelations between subsystems and are defined as functions of (x j , u j ) for i = j; u i ∈ R m i denotes the control input action of the i-th subsystem, y i = h i (x i ) ∈ R p is the subsystem's output. For example, for a formation of land vehicles, h i corresponds to the coordinates of each vehicle. It is considered that p is the same for all subsystems. Formation is defined in the coordinate system moving along with the desired trajectory y d (s) ∈ R p , where s is a certain parameter, and the orthonormal vector F (s) = e 1 (s), . . . , e p (s) ∈ R p , specifying a moving RMS, start which lies at the point y d (s). The formation F is thus specified by k by the points P 1 , . . . , P k ∈ F (s), F = {P 1 , . . . , P k }, where P i = ∑ p j=1 α ij e j , and in the general case it changes over time.
The key parameter is defined as the action reference s, on the basis of which the desired trajectory y d (s) is formed. To calculate the reference action, the reference projections are used, which transform the sensor data into a reference action s, on the basis of which the time is then calculated, which is used to coordinate lower-level feedbacks. The synthesis of the control law consists of several steps. On the first of them, given y d (s) and {P 1 , . . . , P k }, the desired trajectory y d i (s) = y d (s) + ∑ p j=1 α ij e j (s) for each agent (subsystem) in the formation is found. The speed of the formation along the trajectory y d (s) is specified by means of a strictly increasing function s = v(t). The formation control law u = u(x) must ensure the fulfillment of the control goal as lim t→∞ y i (t) − y d i v(t) = 0. With this control, there is an initial state x 0 = [x 01 , . . . , x 0k ] T such that the corresponding trajectory x di (s) coincides with the given one, that is, h i x di (s) = y di (t). In the second step, the control law u i = u i (x, t) is synthesized separately for each subsystem using existing tracking methods ( signal tracking) or movement along a given trajectory ( path following). Different methods of controller synthesis can be used for different subsystems. To coordinate the subsystems, the third step uses the projection mapping. A reference projection is defined as a transformation s = γ(x) such that γ x d (s) = s, that is, if the state x d (s) is on the desired trajectory, then γ must give the corresponding value of s also on the desired trajectory x d (s). For example, for any state x 0 , vector x d (s 0 ) can be an orthogonal projection of x 0 onto x d (s 0 ). Orthographic design is not the only way to define γ, and changing it can fundamentally change the way of subsystems mutual coordination. The last stage of the controller design is the synthesis of a non-time-based feedback law that is used to control a multi-vehicle system. The control law is obtained by simple substitution and has the form A number of examples for selecting reference projections to ensure coordination of the formation are presented in [111]. In particular, the satellite clusters, forming a highresolution synthetic aperture imaging formation are described. A flat formation consists of a group of satellites occupying the same orbital plane OXY in the inertial terrestrial coordinate system OXYZ and separated by the mean anomaly. The desired trajectory of movement of each satellite is considered to be a circular orbit of radius r 0 as x d = r 0 cos(ωt), y d = r 0 sin(ωt). (45) As is known, ω = µ r 3 0 1/2 , µ = 3.986 × 10 14 m 3 / s 2 .
Following [112], the synthesis of the control law in [111] is performed by the feedback linearization method. For a near-circular orbit, the reference action s is defined as s = θ, where θ is the angle between the satellite radius vector and X-axis on the OXY plane. For the desired trajectory, θ = ωt. The relative position of satellites in the formation is determined through γ i . For example, one can accept where θ 0 is the desired angle between two satellites in the formation. According to (46), satellites #1, #3 are the leaders of the formation. Kang et al. [111] presented the theorems on the stability of the formation and the simulation results of the formation reconfiguration during movement. Milam et al. [113] discussed the positioning and reorientation control of a formation of fully steerable low-thrust microsatellites. A general control methodology based on optimization is proposed for solving the problems of constructing a limited trajectory for positioning and reorientation. Taking advantage of the all-wheel drive microsatellites, onboard calculations can be performed on them. The software package Nonlinear Trajectory Generation (NTG) was used for typical space missions.
A systematic approach for studying the formations of multi-agent systems is presented in [114]. In particular, undirected formations for centralized formations and directed formations for decentralized formations, are considered. Paper [114] focuses on the feasibility problem: given the kinematics of multiple agents together with inter-agent constraints, one should find out if there are agent trajectories satisfying constraints. If the answer is yes, the task is to find a "smaller" control system that supports the movement of the formation along its trajectory, which, from the side of higher control levels, allows one to consider the formation as a whole. In [114], n heterogeneous agents with states x i ∈ M i , i = 1, . . . , n are considered. The agent's kinematics are described by drift-free controlled distributions [115,116] on manifolds M i as where U i is the control space, the vector fields X i form the distribution basis. Controlled distributions are general enough to describe an underactuatio. The formation of a group of agents is determined by its graph, which completely describes the kinematics of individual agents and global inter-agency constraints. The graph of the formation F = (V, E, C), consists of: a finite set V of vortices, the cardinality of which is number of agents, where each vertex v i : M i × U i → TM i is a distribution ∆ i describing kinematics of each individual agent according to (47); binary relation E ⊂ V × V, representing communication between agents; families of constraints C indexed by the set E, C = {c e } e∈E . Two different types of formation graphs are considered: undirected formations, where (V, E) is an undirected graph, and directed formations, where (V, E) is a directed graph. In [114], conditions are obtained for determining formation feasibility for two types of formations, and also a management abstraction for a group is found, which allows one to describe a formation as a single object. When directional formation is not feasible it is possible to extract a feasible formation by eliminating the degrees of freedom that cannot be worked out by the follower agents.
In [117], nonlinear equations of orbital dynamics of relative motion are introduced for the problem of motion in a formation with separation into constant distances. This general equation of orbital dynamics allows for elliptical, non-coplanar maneuvers at large distances between spacecraft, as well as typical near-circular, coplanar maneuvers at short distances. In addition, changes in the equations of position and velocity have been introduced for a path in the plane of a formation with a large, constant separation angle between satellites. A nonlinear control method with the Riccati equation depending on the state is used to solve the problem of controlling the movement of the formation. This new control method for a non-linear system provides a clear design trade-off between control action and state error, similar to the classical linear quadratic control method. Numerical modeling shows the effectiveness of the new control method with the Riccati equation depending on the state with the developed equations of relative motion.
Palmer [118] presented a general analytical formulation of the problem for optimal trajectories of satellites flying in a formation, based on the circular version of Hill's problem. Optimization is performed to minimize the displacement energy received from the motor. The optimization problem is the choice of the trajectory along which the satellite should move during the maneuver, while the required time and boundary conditions are fixed. Low thrust engines are suitable for maneuvering in formation. It is assumed that the engine is running throughout the entire maneuver, and the optimal program calculates the magnitude and direction of the thrust as a function of time. The possibility of using your own dynamics to obtain additional fuel economy by changing the boundary conditions is considered.
A method for constructing a satellite trajectory capable of creating a formation of a given configuration from identical spacecraft is given in [119]. The method uses a behaviorbased approach to implement autonomous and distributed control of relative configuration using limited sensor information. For each satellite, the desired speed is defined as the sum of the various components of the general high-level behavior patterns. Further, the model of behavior is determined by an inverse dynamic calculation called the formation of equilibrium. Introduced several feedbacks for speed control. The results of the study show that in the developed method it is advisable to use control in a sliding mode.
Kumar et al. [120] discussed achieving a given formation using thrust only in the direction along the trajectory. A system consisting of the leading satellite and the following satellites is considered, and linear control law is developed to obtain a formation of a given type. The orbital motion of the leader satellite is determined by the radial distance r c from the center of the Earth and the true anomaly θ. The movement of the slave satellite is described with respect to the movement of the leader satellite using the relative r.m.s. S-xyz, anchored in the center of the leader satellite. The x-axis is directed along the local vertical, the z-axis is taken along the normal to the orbital plane, and the y-axis represents the third axis of this selected right RMS. The geometry of the orbital motion of the lead and slave satellites according to [120] is shown in the Figure 2. The orbital equations of motion for the lead satellite and the equations of motion of the slave satellite relative to the lead satellite are of the form [21,23]r where r = (r c + x) 2 + y 2 + z 2 1/2 ; f y is the acceleration of the control force in the longitudinal direction; f dj , j ∈ {x, y, z} are the accelerations of the disturbing force along the directions x, y and z, respectively. To synthesize the control law, the following linearized model of the (50)-(52) system is used, obtained under the assumption that the leading satellite is moving in a circular orbitẍ whereθ = const is taken. Figure 2. Geometry of the orbital motion of the lead and trailing satellites according to [120].
In [120] the following linear proportional-differential (PD) control law is proposed where x d , y d are the desired (specified) values of the coordinates x, y, k p1 , k d1 , k p2 , k d2 are the regulator coefficients. The performance of the proposed regulator was investigated by the Routh-Hurwitz criterion and tested using the simulation of nonlinear equations of motion of the system taking into account several factors, including changes in the initial conditions, the size of the formation, the presence of disturbing forces and eccentricity of the orbit of the lead satellite. It was found that a regulator using only thrust along the trajectory can provide a bounded relative position error, without damping, with a maximum control acceleration (in m/s 2 ) of 0.0718θ 2 times the error along the trajectory and, therefore, has a limited ability to neutralize disturbances. However, by modeling, it was found that the proposed controller successfully provided limited relative position errors below ±10 m.
Sengupta et al. [121] presented expressions describing the average relative motion of two satellites in adjacent orbits around a planet flattened at the poles. The theory assumes small relative distances between satellites, which is equally true for all elliptical orbits, as well as for the special case of a circular orbit due to the use of nonsingular orbital elements. An analytical filter has been introduced, which compensates for short-period changes in the relative state without using built-in digital filters, which are usually used to filter noise in the design of control systems. The application of this filter to maintain the formation on a given relative trajectory is considered.
In [122], orbital elements are defined as a set of parameters that characterize the relative motion of a satellite formation based on the principle of geostationary spacecraft control. From these parameters, the orbital motion and maneuvering equations are derived. In addition, in [122], acceleration control laws are proposed, including out-of-plane and in-plane control. To demonstrate them as orbital maneuvers, the following are considered: the creation of a formation of two satellites, its reconfiguration, maneuvering at long distances, and the preservation of the formation. The results of modeling the dynamics of the movement of the formation under the influence of the gravitational field and eccentricity when applying the proposed law of control with feedback are presented, which show the effectiveness of the proposed method and the possibility of correct control of the slave satellite from the lead satellite.
In [123], simple and general analytical solutions are considered for the optimal reconfiguration of a satellite constellation controlled by various linear dynamic equations. The calculus of variations is used for the analytical search for optimal trajectories and control. The proposed method makes it possible to predict in advance the exact form of optimal solutions without the need to solve the problem. It is shown that the optimal thrust vector as a function of the fundamental matrix of the given equations of state. The analytical solutions presented in this paper can be applied to most problems of satellite reconfiguration governed by linear dynamic equations. Numerical modeling confirms the brevity and accuracy of the solutions obtained.
During the LISA missions [9], it may be necessary to reconfigure the formation to maintain or improve its functionality [10]. This reconfiguration can be described by relative dynamics, as presented, for example, in [23,124]. Using relative dynamics in projection onto a circular orbit, the authors of [17] found analytical solutions for optimal two-pulse trajectories of reconfiguration of a formation of two spacecraft in orbit. The desired formation is characterized by nonsingular different orbital elements. To achieve the desired orbital-element differences, [17] proposes an analytical two-pulse solution found using the Gaussian variational equations. The solutions obtained can be easily implemented using onboard computing devices. In order to maximize mission life during reconfiguration, many papers consider the problem of optimization in terms of fuel consumption. The works by Kurzhansky and his colleagues significantly developed a solution to the problem of control synthesis for impulsive systems, based on generalizations of the Hamilton-Jacobi-Bellman variational inequalities, which made it possible, within the framework of a single formalization, to study control problems for hybrid systems containing jump-like rearrangements of states [125,126].
Optimization methods based on the calculus of variations, including the theory of basis vectors, encounter difficulties in finding the global optimal solution, since these methods can only search in a convex neighborhood around the initial approximation. Therefore, the global optimality of the solution is not guaranteed, especially with substantially non-convex objective functions. To mitigate these difficulties, several methods have been proposed, one of which is the "hybrid optimization"-a combination of the calculus of variations and global extremum search [127,128]. Paper [10] is apparently the first attempt to use hybrid optimization with a genetic algorithm to find reconfiguration trajectories, although the genetic algorithm has often been used for other trajectory optimization problems. In [10], the relative motion in the LVLH reference frame is described using the HCW equations of motion.
The satellite constellation consists of a head satellite and surrounding slave satellites. The orbit of the lead satellite is considered a reference orbit, and the relative orbits of the follower are considered to be PCO. PCO predicted distances between the master and slave satellites [129] are constant. The HCW equation describes the relative motion between the master and slave satellites in the local reference frame LVLH with the origin on the master satellite. In this system, the OX axis is defined in the radial direction of the main orbit, and the OY axis is defined as the longitudinal direction of the main orbit. The OZ axis is defined to complete the right coordinate system. The relative position vector r and the relative velocity vector v of the slave satellite are defined in the LVLH reference system. State vector x of relative dynamics is defined as Then, [10] displays the linearized equation of relative motion for two bodies. The distances between the master and slave satellites are ignored during linearization since they are significantly less than the distance between the main satellite and the center of gravity (Earth, when moving in the LEO). Thus, the linear equations of state are obtaineḋ where O 3 , I 3 are the zero and unit matrices of order 3, n = √ GM/a 3 is the average speed of the lead spacecraft, G is the gravitational constant, M is the mass of the central body (for the Earth, GM = 398, 603 · 10 9 m 3 · c −2 ), a is the semi-major axis of the orbit of the lead satellite. The homogeneous Equations (57) are then analytically integrated using the t is the fundamental matrix of the system (57) , t 0 stnads for the current time, x(0) denotes the initial state vector.
(There are some inaccuracies in this expression for matrix Φ in [10].) Kim et al. [10] explored the three types of reconfiguration maneuvers described in [73] namely resizing, reassigning, and reorienting. Resizing is the change in distance between spacecraft in the formation. In [10], a hybrid optimization strategy with a genetic algorithm and basis vector theory is used for fuel-optimal reconfiguration of a satellite group. In hybrid optimization, the genetic algorithm performs a global search to create two-impulse trajectories at each transition time, and vector analysis of basis vectors serves for local optimization with two-impulse trajectories as initial approximations. The final optimal trajectory has multiple pulses and also lower fuel consumption than two-pulse trajectories. The mission planner must investigate the trade-off between fuel saved by multiple impulses and the increased mission complexity due to additional impulses. In many missions, the trajectory obtained by the genetic algorithm will have an advantage for the transition due to its simplicity and almost optimal fuel consumption.
The Udwadia-Calaba method is used in [130], where under the assumption of unlimited and time-varying throttle starts throughout the entire maneuver, precisely and explicitly without any restrictions on the distance between satellites, the formation retention problem is solved. Cho and Yu [130] discussed a system of two satellites orbiting a central body (around the Earth). One of them is called the "master", and the other is the "deputy" satellite, which must maintain a given location relative to the head. The equation of motion of two bodies for an inertial frame of reference is used, which is supplemented by the requirement to hold the station as constraints for obtaining an equation of motion that fully takes into account all nonlinearities. From this, the governing force is derived as an explicit function of state and time required to maintain the formation. The results obtained are also applicable to the case when the distance between the satellites is so great that the linearized relative equations of motion are not applicable. The dynamics described in the inertial frame of reference is then transformed for convenience into a local moving frame of reference through the transformation matrix. As a practical limitation, in [130] a relative configuration that is circular when projected onto the local horizontal plane is chosen. It is commonly referred to as the PCO [129,131]. Unrestricted motion is relative motion without a control force acting on the satellite system. The master and slave satellites move only under the influence of gravity, and therefore "unlimited" and uncontrolled motions have the same meaning. In the case of movement with restrictions, control forces must be applied to the system. Therefore, "limited" and "controlled" movements also have the same meaning. When a satellite with a point mass orbits the Earth, its acceleration a is written in the Earth's inertial coordinate system as follows [132]: where G is the universal gravitational constant, M is the mass of the central body (the mass of the Earth). The Equation (58) is written in the inertial terrestrial coordinate system originating from the center of the Earth, the X axis points to the vernal equinox, the Z axis passes through the North Pole, and the Y axis completes the right n.v. (see Figure 3).  [130]. (Earth s.c. (X-Y-Z) and s.c. Hill (x-y-z). The distance between the leading and trailing satellites on the yz-plane of the local rotating coordinate system x-y-z must be equal to the specified ρ).
Assuming that the master satellite is in a circular orbit with a fixed inclination without any disturbance, the goal is to maintain a constant distance between it and the follower, for example, when projecting onto the local horizontal plane. In addition, to prevent discrepancy in X coordinate, a constraint is also applied between the X and Z coordinates. Hill's LVLH coordinate system [21,23] is used to describe these limitations. The origin of this system is located on the leading satellite, X-axis lies in the radial direction, Y-axis is in the path direction, and Z-axis lies along the orbital angular momentum vector, thus completing the right-handed coordinate system see Figure 3. In LVLH coordinate system, the constraint equations can be written as where ρ is the constant distance between satellites. Since the Equations (59) and (60) describe holonomic constraints, they are twice differentiated in time, which gives Further, in [130], transformations of the rotation of variables from the Earth's r.s. to the Hill coordinate system are performed with the rotation matrix R, depending on the longitude of the ascending node of the lead satellite Ω, the angle between the equatorial plane and the orbital plane of the lead satellite i, latitude of the lead satellite ω = ω 0 + n(t − t 0 ), n = GM r 3 0 , and the radius r 0 of the circular orbit of the lead satellite. This circular is based on linearized equation solutions, and therefore non-linear behavior of relative motion will violate this configuration. To avoid this, additional control force is used. The approach of [130] provides a general methodology that can be easily applied to any type of relative configuration. To solve the problem in [130], further transformations of constraints are performed (61) and (62) to the form (36) and then use (42). For this purpose, equations (59) and (60) are first reduced to the form Finally, the system equations are reduced to the form [81]: whence the expression for the control force F c is derived: The presentation is illustrated by numerical examples. It is worth mentioning that in [130] there are no explicit restrictions on either the control action or the fuel consumption.
Simple and accurate formation control schemes are presented by Cho and Udwadia [133], where a method of analytic dynamics by Udwadia [76], Udwadia and Kalaba [77], Udwadia [78,79] is employed. In [133], the explicit control inputs for accurately maintaining a given formation configuration, using continuous thrust systems, are found. A much simpler and more explicit expression is obtained to accurately satisfy the formation conservation constraints for Keplerian reference orbits. Cho and Udwadia [133] also includes explicit control results when the follower is put into orbit with incorrect initial conditions, as is usually the case in practice. Extensive computational simulations demonstrate that the approach of [133] is easy to implement and numerically precise.
The problems of attitude synchronization and tracking in the spacecraft formation in the presence of model uncertainties and external interference are considered by Wu et al. [134]. A decentralized control law for the adaptive slip mode using the topology of undirected communication between spacecraft is proposed, which is analyzed on the basis of the theory of algebraic graphs. A sliding manifold is obtained for several spacecraft, on which each spacecraft approaches the desired time-varying orientation and angular velocity, while maintaining orientation synchronization with other spacecraft in the formation. Next, a control law is developed to ensure convergence to a sliding manifold. The stability of the resulting closed-loop system is proved by applying the Barbalat lemma. The simulation results demonstrate the effectiveness of the proposed synchronization and tracking methodology for orientation problems.
The problem of designing a distributed controller for a continuous time system consisting of a number of identical dynamically coupled subsystems is discussed in [135]. The main mathematical conclusion is based on the Kronecker product and a special similarity transformation built from a matrix of interconnection patterns that breaks the system into modal subsystems. This, along with the result of formulating the extended linear matrix inequality (LMI) for continuous-time systems, allows one to obtain explicit expressions for calculating the parameters of distributed controllers for both state feedback and dynamic output feedback. The main contribution of [135] is the solution of the distributed control problem for systems with continuous time under the conditions of H ∞ /α-stability and H 2 /αstability characteristics by solving the set of LMIs with non-common Lyapunov variables. The efficiency of this method is demonstrated by the example of a satellite formation.
Autonomous navigation of micro-UAVs is usually based on the integration of inexpensive Global Navigation Satellite System (GNSS) receivers and inertial and magnetic sensors based on microelectromechanical systems (MEMS) for motion stabilization and control. The resulting navigation performance in terms of position and orientation accuracy may be insufficient for other mission tasks, such as those related to precise sensor pointing. In this context, [136] introduces a collaborative UAV navigation algorithm that allows the main vehicle equipped with inertial and magnetic sensors, a global positioning system (GPS) receiver, and a vision system to improve its navigation performance (in real time or in real-time). The focus is on the outdoor environment, and the key concept is to use differential GPS (DGPS) among vehicles and vision-based tracking (DGPS/Vision) to create a virtual auxiliary navigation sensor, which information is then integrated into a sensor fusion algorithm based on an extended Kalman filter (EKF). The developed concept and processing architecture are described with an emphasis on the DGPS/Vision orientation determination algorithm. Efficiency assessment is carried out on the basis of numerical modeling and flight tests. In the latter of these, the navigation estimates obtained from the DGPS/Vision approach are compared with the estimates provided by the onboard autopilot system of the tuned quad. The analysis shows the potential of the developed approach, mainly due to the ability to use magnetic and inertial-independent accurate orientation information.
The distributed robust control method is developed relies, on a spaceborne distributed telescope, is presented in [137], where a combination of robust H ∞ control and distributed control using the consensus approach is proposed, focused on the application to satellite formation motion.
The method of simultaneous control of the relative motion and orientation of the solar radiation pressure was proposed in [138]. The aim of the control is to stabilize a given closed relative orbit. The basic idea is to use special materials for solar sails that can change their optical properties. A solar sail consists of several rectangular cells. Each of them can change its optical properties. The required steering force is generated by varying the average reflectivity of the solar sail, and the steering torque is achieved using the appropriate black and white mesh pattern.
Mashtakov et al. [139] discussed the use of sunlight pressure to control a satellite formation. The control is based on the rotation of the sail normal, which provides for a change in the optical properties of its surface.
Monakhova et al. [140] dealt with the problem of satellite formation immediately after launch. When detached from the launcher, there is a certain error in the magnitude of the speed jump, which leads to slightly different periods of rotation of the satellites. As a result, relative paths become open. Differential control using a resistance force is considered. The satellite orientation is controlled by a magnetic actuator.
The problem of creating a robust formation controller for a group of satellites whose dynamics are associated with nonlinearities and uncertainties is considered by Liu et al. [141]. For a constellation of satellites, a formation adjuster is proposed that includes a position adjuster for the formation of the desired precise formation and a position adjuster for aligning the positions of the satellites. It is shown that errors in trajectory tracking and orientation of the entire closed-loop control system can converge to a given neighborhood of the origin in a finite time. Simulation results are presented to demonstrate the benefits of the proposed formation control scheme.
In [142], trajectory analysis, mission design, and the control law for several microsatellites for joint traversal of the spacecraft are considered. The problem of the cooperative movement of a group of satellites along an elliptical trajectory, i.e. cooperative circumnavigation, (CCN) is first defined to bring the group of microsatellites to a predetermined plane ellipse relative to the host spacecraft while maintaining the geometric configuration of the formation. Then, a CCN control law is proposed, in which artificial potential functions and the Laplace affine matrix are applied to fulfill the requirements of the CCN mission. Simulation results are presented to show the feasibility of this approach.
In [143], a joint navigation approach for unmanned aerial vehicles (UAVs) is presented, which makes it possible to reliably and accurately determine the orientation of the main aircraft flying in formation with other UAVs (alternates, slaves). The proposed method is based on a closely related extended Kalman filter using the spatial diversity of measurements from global navigation satellite systems (GNSS) and vision systems integrated with data from inertial and magnetic sensors. The main focus is on the external environment, and the innovative idea is to extend the attitude estimation approaches based on multiple GNSS antennas to a multipurpose system where differential GNSS and UAV visual tracking of UAVs are used to create a virtual auxiliary navigation sensor. The processing concept and architecture are described with an emphasis on an ECP measurement update phase that is applicable to any number of collaborating proxies and to different GNSS processing architectures. The effectiveness of the proposed method is evaluated by experimental tests using two multi-rotors and two fixed ground antennas, one of which is used as a ground reference point for the analysis of pointing accuracy. The results show the feasibility of the developed approach in terms of accuracy and providing drift-free estimates in real-time or in post-processing scenarios.
A group of four satellites forming a tetrahedron is discussed in [144]. The task is to maintain relative orbits so that the tetrahedron retains its shape and size over time according to the introduced criterion. The problem of constructing relative reference orbits of satellites is solved, which, within the framework of the linear HCW model, ensure the preservation of the tetrahedron parameters in the absence of control. In the framework of the Schweigart-Sedwick linear model [91,92], which takes into account the influence of the second harmonic of the geopotential, as well as in the framework of a nonlinear motion model with the second harmonic, algorithms are proposed for constructing uniaxial motion control of satellites oriented along the Earth's magnetic field, which ensure the maintenance of relative orbits.
In [145], the task of controlling a group of four satellites in a low near-circular Earth orbit is considered. The satellites form a tetrahedron of a certain shape and size. The initial data were selected in such a way that the change in the shape and size of this tetrahedron during the movement of satellites in orbit was minimal. The choice of the initial data was carried out analytically and then numerically refined. Michael Ross et al. [146] presented an approach to solving a class of problems arising in the development of satellite swarms.
The key issue addressed in this article is "simultaneous" fusion and orbit control to achieve swarm configuration. Although any design criterion can be used, a fuel consumption approach is used, since the fuel cost (loss) function is extremely high for a spacecraft (SC). It is shown how certain elements of the theory of optimal periodic control provide a very natural formulation of this problem. Using versatile dynamic optimization software it is shown how one can easily design satellite formations without using any analytical results. If a natural zero fuel solution does not exist, a byproduct of the proposed approach automatically determines the minimum amount of fuel and the associated controls required to maintain the formation.

Conclusions
This paper provides an overview of the results on the control of a formation flight of spacecraft, expanding and supplementing the review [1], related to 2004, with modern publications. Three types of formation architectures are considered: MIMO, in which the formation is treated as a single entity with multiple inputs and multiple outputs; LF, in which individual spacecraft controllers have linked hierarchically, and a virtual structure in which each spacecraft in formation is considered as a rigid body, embedded in a common virtual structure. The MIMO architecture differs from other ones in its optimality and stability. However, the use of all system states leads to the need to introduce high requirements on the information exchange, so such algorithms are usually unstable to local failures. The leader-follower architecture uses information about leaders only, which simplifies formation navigation. The problem of reliability is solved by increasing the number of leaders. However, this approach is inferior in optimality, and for ensuring high system stability, information requirements can reach those of the MIMO case. The cyclic architecture is a middle ground between LF and MIMO architectures. The information requirements are as high as those of MIMO case, and connections between some devices are allowed. However, round-robin algorithms can be completely decentralized. There is no coordinating agent and no instability due to single failures. In the process of research in this area, a number of directions have appeared in which, in particular, an increase in the stability and reliability of systems, reducing requirements for information restrictions, and increasing the spacecraft autonomy in formations. In the paper, examples of modern group spacecraft missions are also outlined.

Funding:
The work was carried out in BSTU "VOENMEH" with financial support from the Ministry of Science and Higher Education of the Russian Federation (Government contract agreement No 075-03-2020-045/2 of 9 June 2020).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: