Control of Two Satellites Relative Motion over the Packet Erasure Communication Channel with Limited Transmission Rate Based on Adaptive Coder

: The paper deals with the navigation data exchange between two satellites moving in a swarm. It is focused on the reduction of the inter-satellite demanded communication channel capacity taking into account the dynamics of the satellites relative motion and possible erasures in the channel navigation data. The feedback control law is designed ensuring the regulation of the relative satellites motion. The adaptive binary coding/decoding procedure for the satellites navigation data transmission over the limited capacity communication channel is proposed and studied for the cases of ideal and erasure channels. Results of the numerical study of the closed-loop system performance and accuracy of the data transmission algorithm on the communication channel bitrate and erasure probability are obtained by extensive simulations. It is shown that both data transmission error and regulation time depend approximately inversely proportionally on the communication rate. In addition the erasure of data in the channel with probability up to 0.3 does not inﬂuence the regulation time for sufﬁciently high data transmission rate.


Introduction
In recent years, there has been a growing interest in using the differential force (i.e., the difference between the aerodynamic drag forces applied to the satellites) to eliminate the relative drift between the satellites in a swarm moving in a group (without the mandatory requirement to maintain relative position, cf. [1][2][3]). Various control algorithms using differential aerodynamic drag have been proposed in numerous publications, see [4][5][6][7][8][9][10][11][12][13][14][15][16]. One of the fundamental publications is the work by Leonard [17] where, based on the assumption of the possibility of changing the effective cross-section of satellites, a method of switching control has been developed. The differential force is created by changing the attack angles of the plates, located on the satellites, due to the rotation of the satellites with respect to the incident airflow. Kim et al. [18] deal with a satellite constellation, consisting of a leader satellite and surrounding slave ones. The orbit of the leader is considered as a reference, whilst the relative orbits of the followers are considered to be the Projected Circular Orbit (PCO), which is the relative orbit between the master and slave satellites, the data transmission rate for the given power consumption can make it possible to increase the power of transmitted signals, expanding the area of inter-satellite interaction.
The present paper is focused on the reduction of the inter-satellite demanded communication channel capacity taking into account dynamics of the satellite's relative motion and possibility of erasing the navigation data in the channel.
The limitations of control under constraints imposed by the limitations of the communication channel capacity, have been deeply studied within the control theoretic literature, see [31][32][33][34][35] and the references therein. A fundamental result establishing the smallest value for which the stabilization (estimation) problem for linear time invariant (LTI) systems is obtained by Nair and Evans [31] and presented in the form of the seminal Data Rate Theorem. However, the study of the application control problems in aerospace under communication constraints is still very limited.
In the present paper, a system of two coupled satellites is considered. As in [21], the satellites are assumed to be launched at the starting time in accordance to the specified separation conditions. It is assumed that the satellites move in a low circular near-Earth orbit and are controlled using the aerodynamic drag force, which is achieved by rotating the satellite relative to the incoming flow using a flywheel attitude control system. The main focus of the paper is in the navigation data exchange between the satellites to be used to keep the satellite motion in a swarm. To this end, the adaptive coding procedure is proposed and is studied for the cases of ideal and erasure communication channel. The regulation time is taken as the performance criterion, and its dependence on the data transmission rate is numerically studied.
The reminder of the paper is organized as follows. The existing results on control and estimation under information constraints are briefly recalled in Section 2, where the minimum necessary data rate for the estimation and control of Linear Time Invariant (LTI) systems in the form of the data rate theorem is specified and various coding/decoding schemes are described. Section 3 is devoted to the dynamics of two satellites' relative motion in a near-circular orbit. The main result is concentrated in Section 4. This section starts with designing the control law, which ensures the asymptotic regulation of the satellites relative motion (Section 4.1). The design is based on the linearized dynamics model without taking into account the control signal saturation. The classical modal control approach based on the pole-placement technique is employed [36]. The behavior of the system with the saturation in control is studied in the subsequent sections by the simulations. The next stages of the present study are dedicated to the evaluation of the proposed scheme of the inter-satellites data transmission over a digital communication channel, aimed for reducing the necessary channel capacity. To this end, the adaptive coding procedure for the position transmission between the satellites in the formation, employing the kinematics process description, is introduced in Section 4.2. It is worth mentioning that the application of this procedure makes it possible to avoid measuring the time derivatives of the satellites relative position (these derivatives are needed for control) due to the state observer, which is embedded into the adaptive coder/decoder pair. Then in Section 4.3, the model of the erasure communication channel, adopted in the present research, is described. Results of the numerical study of the closed-loop system performance and accuracy of the data transmission algorithm dependence on the communication channel bitrate and erasure probability obtained based on the extensive simulations are presented in Section 4.4. Concluding remarks and the future work intentions given in Section 5 finalize the paper.

Problem Description
Let us consider control and observation (estimation) systems containing a digital communication channel. For these systems, plant output measured by the sensor at discrete instants t k = kT 0 , where T 0 denotes the sampling interval, k = 0, 1, . . . , is converted by the coder into the characters of the coding alphabet S. The sequence of characters is transmitted over a digital communication channel to the decoder. The decoder transforms messages from the transmitted form in a form adequate to subsequent calculations and transformations by the controller. In [37][38][39] the communication channel was considered of limited capacity, but was otherwise ideal. The cases of packet erasure channel and 'blinking' channel widely appear in various real-world systems, see, e.g., [40][41][42][43][44][45][46][47][48][49]. Therefore, for a more realistic analysis, the properties of the communication channel, such as distortion, erasure, and data loss, should be taken into account. In the present study, the effect of data erasure is considered.
Signal quantization introduces essentially non-linear properties into the system, characterized by the presence of the dead zone, discontinuities, and saturation (associated with bit grid overflow). Additionally, the signal sampling on time involves the hybrid (continuous-discrete) system description. A rigorous examination of the influence of time sampling and the level quantization is a complex nonlinear analysis problem, that usually does not have an exact analytical solution. In the early studies, the level quantization in digital control systems was usually considered a source of independent additive random noise affecting the system. This assumption makes it possible to significantly simplify the study of level quantized systems, especially for LTI plants. However, if the quantization level is relatively high (for example, in the case of binary quantization), this can lead to the emergence of self-oscillations and even the system divergence, see [50][51][52][53]. Therefore, to analyze the system, its nonlinear model is required. Besides, the possibility of the bit grid overflow can also affect the quantizer, as a result of which the saturation is introduced to the control loop [54,55].

Minimum Necessary Data Rate for Estimation and Control
The limitation of the data transmission rate over the communication channel can be expressed in informational terms. Assume that coding alphabet S consists of µ elements. Then at each step k = 0, 1, . . . over the channel an amount ofR = log 2 µ bit can be transmitted per each step. Let the data transmission be carried out at discrete instants t k = kT 0 , where T 0 is the sampling interval. Then the data transmission rate in bits per second is as R = T −1 0 log 2 µ bit/s. In this regard, one speaks of "information constraints" in control and estimation problems.
The problem of determining the minimum bandwidth of the communication channel, at which it is possible to provide the required estimation accuracy, is posed and partially solved by Nair and Evans [56]. The sufficient condition for the value obtained in Nair and Evans [56] was developed in subsequent works in the form of the Data Rate Theorem which is a fundamental result establishing the smallest value for which the stabilization (estimation) problem for linear systems is solvable in principle. Nair and Evans [31] studied the exponential stabilizability of LTI plants in the sense of achieving an exponential moment stability. For a deterministic initial state case the result of [31] can be roughly presented in the following form [57].
Let the LTI discrete-time plant be described by the difference equation: where x[k] ∈ R n , y[k] ∈ R l , and u[k] ∈ R m are the state, output, and control vectors, respectively; A, B, and C are the matrices of the corresponding dimensions; and k ∈ Z + denotes the step number (the discrete time). It is assumed that pair (A, B) is reachable and (C, A) is observable. Let the sensor be connected to the controller over a digital communication channel, and no more thanR bits of data can be transmitted at each step k. Then the necessary and sufficient conditions for ρ-exponential (with the prespecified stability bound ρ > 0) stabilization are given by the inequality [31]: where η j are the eigenvalues of matrix A , j = 1, . . . , n. The right-hand side of (2), denoted as , gives a tight admissible bound when ρ-exponential stabilization can be achieved. For real-time systems with the constant sampling interval T 0 , NE-number R NE in bits per second has a form: Extensions of this result for stabilizing nonlinear systems in the vicinity of the origin and observing nonlinear systems through finite capacity communication channels, including large networks, were obtained in the series of the subsequent papers, see [58][59][60][61][62] to mention a few.
It is worth mentioning that on practice, the data transmission rate can not be taken as small as the NE-number (3) gives due to several reasons for the data bitrate being usually much greater than R NE , but the NE-number can serve as a measure showing the maximum available possibilities of estimation and control over the existing communication channel. Furthermore, a promising approach is to employ the event-triggered control instead of the control with constant sampling time, see e.g., [63][64][65].

Coding/Decoding Schemes
Under the assumption that the sampling time T 0 can be chosen arbitrarily, optimality of the binary coding in the sense of the required transmission rate (in bit-per-second) has been proven in [39], see also [66]. Therefore in the present study the binary quantizer is used as a core element for the coding procedure.

Static Binary Quantizer
Let σ[k] be a scalar information signal to be transmitted over the digital communication channel in discrete instants t k = kT 0 , where k = 0, 1, · · · ∈ Z + is a sequence of natural numbers, and T 0 is the sampling interval. Let us introduce the following static quantizer: where sign(·) is the signum function: Parameter M is referred to as a quantizer range. The output signal of the quantizer is represented as one-bit information symbol from the coding alphabet S = {−1, 1}, and is transmitted to the decoder. Note that for the binary coder, the transmission rate is as R = 1/T 0 bit per second. It is assumed that the equi-memory condition is fulfilled, i.e., the coder and decoder make decisions based on the same information [67,68]. The binary output codeword s ∈ S is transmitted to the decoder.

Zooming Strategies
In time-varying quantizers [33,39,[69][70][71][72] range M is updated with time. Using such a zooming strategy improves the steady-state accuracy of the transmission procedure and at the same time prevents the encoder saturation at the process beginning. The values of M[k] can be precomputed (the time-based zooming) [39,73,74], or current quantized measurements can be used at each step for updating M[k] (the event-based zooming). For an audio channel, Moreno-Alvarado et al. [75] developed the coding schemes with the capacity to simultaneously encrypt and compress audio signals, which makes possible increasing necessity for transmitting sensitive audio information over insecure communication channels.
The event-based zooming can be realized in the form of the adaptive zooming [76][77][78][79], where the quantizer's range is adjusted automatically depending on the current variations of the transmitted signal.
For the binary quantizer the following adaptive zooming algorithm was proposed and experimentally studied in [78]: where

Coders with Memory
The coding/decoding procedure can include the embedded observer, which adds a memory to the coder. The following model of drive process is used: where x(t) ∈ R n is the process state space vector; y(t) is the scalar measured signal; A ∈ R n × n , B ∈ R n × 1 are given real matrices; and ϕ(t) is the external input signal which is assumed to be the same both on the transmitter and the receiver sides (so that the equi-memory condition be fulfilled).
The quantized observation errorσ[k] is defined as a deviation between measured signal y(t) and its estimateŷ(t) quantized with given M[k] as follows: where the estimateŷ(t) is generated by the following observer: wherex(t) ∈ R n is the state estimation vector;ŷ(t) is the estimate of the drive process; L is (n × 1)-matrix (the column-vector) of the observer parameters; and the continuous-time observation errorσ(t) is found as an extension ofσ[k] over the sampling interval. In the case of the zero-order extrapolation,σ(t) has a form

Relative Motion Dynamics of Two Satellites in a Near-Circular Orbit
For describing the dynamics of the relative motion of two satellites, moving in near-circular orbits at the Earth's central gravitational field, equations in the Local Vertical Local Horizontal (LVLH) reference system in relative coordinates according to the HCW model is used, see [21,[80][81][82][83][84].
In the present paper, the LVLH reference system is employed, where the OZ axis is directed from the center of the Earth, the OY axis is directed normal to the orbital plane, and the OX axis complements the triple to the right coordinate frame, see Figure 1. Taking into account that the motion along the normal to the orbital plane (along OY axis) is isolated, and proceeding from the universal gravitation law and Kepler's third law under the assumption that for distance ρ of the satellite to the center of the Earth inequalities ρ x, z are held, one obtains the following HCW equations of satellite motion in the LVLH local coordinate system: where F x , F z are the components of the vector of non-gravitational forces applied to the satellite, expressed in acceleration units. These forces include the aerodynamic drag force f x along the OX axis, which is controlled by the satellite turning through the attack angle α. Based on the physical considerations it is valid that f ). Therefore it is natural considering angle α within 0 ≤ α ≤ π/2. In (10), ω denotes the averaged angular velocity of the spacecraft in orbit and satisfies the expression ω = µ/a 3 where µ = GM, G is the gravitational constant, M is the mass of the central body (for the Earth, µ = 398,603 ×10 9 m 3 s −2 ), and a is the semi-major axis of the satellite orbit. HCW Equations (9) and (10) are derived under the assumptions: √ x 2 + z 2 is small compared to ρ; the aerodynamic force is small, the acceleration caused by it does not exceed 1.7 × 10 −6 m/s 2 ; the eccentricity of the orbit is small, the orbit is very close to the circular one; and the orbital rate ω is approximately constant. For the resulting system (9), (10) in [17], a non-degenerate coordinate transformation is performed to the real Jordan form [36], as a result of which the model (9), (10) is represented in the form of a double integrator and a harmonic oscillator, connected in parallel. For the double integrator, a speed-optimal control is created in the form of relay feedback with a quadratic switching function. A similar, but technically more complex approach is used for the harmonic oscillator. Since the control is scalar, it is not possible to apply both control laws simultaneously and independently to both the double integrator and oscillator. To resolve this contradiction, Leonard [17] developed and studied for different flight scenarios an algorithm for switching from one control law to another, based on the specifics of the rendezvous problem requirements. The linearized equations of the satellites relative motion in the OXZ plane (9), (10) in the state-space form read as (cf. [5,21]): where χ = x 12ẋ12 z 12ż12 T ∈ R 4 is the system state vector, where x 12 denotes the difference in coordinates of the second and first satellites along the OX axis, x 12 = x 2 − x 1 ; z 12 = z 2 − z 1 ; and u denotes the control action (difference of aerodynamic forces acting on satellites, in units of acceleration). It is limited in modulus by the value u max . Equation (12) closes the feedback loop. In it, y is the output of the linear controller in the feedback, ϕ(·) is a nonlinear function describing the limitation of the control action, which is assumed to be the saturation function, ϕ(y) = sat u max (y). Since for each satellite the aerodynamic drag force is negative (that is, it acts against the direction of motion along the longitudinal axis and lies in the interval [−u max , 0]), then in order to provide the desired differential control −u max ≤ u ij ≤ u max , the actual steering (braking) action should only be applied to the first (along X-axis) satellite of the pair, see [21]. Matrices of dynamics model (11) have the form: Coefficients k i , i = 1, . . . , 4 are chosen at the stage of the control law synthesis.

Control Law Design
Suppose that each satellite has the onboard navigation equipment, e.g., the GLONASS/GPS receivers [85], for determining its position and is able to share this information via the digital inter-satellite communication channel. Therefore it is assumed that the onboard control system of each satellite is supplied by the values of all the relative coordinates x 12 = x 2 − x 1 , z 12 = z 2 − z 1 and their time derivatives. Control signals for each satellite are generated so as to provide the difference u = u 2 − u 1 , generated in accordance with the chosen control law u = U(x 12 , z 12 ,ẋ 12 ,ż 12 ). In the present study, for control law design the pole-placement technique is employed. The design procedure is made for a LTI system model under the assumption that the restriction on the control signal is not "active", i.e., u does not go beyond the boundaries of the linear region, that is |u| ≤ u max . State-space equations of the closed-loop system (disregarding disturbances) have the form:χ where matrices A, B, K are of form (13). The design problem consists of choosing the coefficients of the controller k i so as to provide the required spectrum {λ ABC } of the matrix (A − BK) of the closed-loop system (14). The fourth order Butterworth polynomial: is used, where parameter Ω is the geometric mean root of the characteristic polynomial. This parameter determines the desired transient time of the closed-loop system. Note, that the controller design is of a secondary meaning for this study that is focused on the data exchange between the satellites, and the other control schemes can be used for improving the formation dynamics. For example, in the case of essential parametric uncertainty, the Implicit Reference Model adaptive control [86,87] or the sliding mode control of [88][89][90] can be also employed.

Adaptive Coding for Transmission of Position Between Satellites in Formation
The application of adaptive coding for navigation data transmission between satellites in the formation is often based on the kinematic representation of the vehicle motion by the second-order model for each channel under the assumption that the satellite speed is constant, which leads to the following representation of the data source generator (cf. [91][92][93]): where y[k], V[k] are the satellite position and speed with respect to the given direction, variable ξ[k] denotes the unmodeled variations of the vehicle speed and is considered as an unknown disturbance.
To estimate the change of the coded signal, the following embedded observer is introduced to the encoder and decoder ( [78,93]):

Erasure Channel Description
By an analogy with [47,93,94], let us assume that output measurement is encoded by an encoder and transmitted to a decoder through packet erasure channel with erasure probability p. In addition, suppose that there exists a feedback from decoder to the encoder for acknowledgment whether the packet was erased or not. Therefore the encoder knows what information has been delivered to the decoder (i.e., the aforementioned equi-memory condition is fulfilled). Let the acknowledgment signal at time k which is sent by the decoder and received by the encoder be represented by ζ[k] ∈ {0, 1} as follows: Random variables ζ[k], k = 0, 1, . . . are assumed to be independent and identically distributed with common distribution: P r (ζ[k] = 0) = 1 − p and P r (ζ[k] = 1) = p.

Satellite Model Parameters
For the numerical investigations, two satellites are represented by 3U-cubesates in the form of rectangular parallelepipeds with dimensions 10 cm × 10 cm × 30 cm. The aerodynamic drag force acting to i-th satellite in the acceleration units is given by the relation where V denotes the satellite running speed with respect to the Earth's atmosphere; ρ = ρ(h) is the air density on the height h of the satellite orbit; α i ∈ [0, π/2] stands for the angle of attack; C α is the drag derivative coefficient with respect to attack angle; ∆S denotes the difference between the satellite cross-sectional areas for the cases of α = 0 and α = π/2 rad; S 0 is the cross-sectional area as α = 0; and m is the mass of the satellite. Therefore the differential drag can be found as which divides the following maximal control action u max = ρV 2 2 C α ∆Sm −1 ms −2 . The satellite running speed V is related to the orbital angular velocity ω as follows: ×10 24 kg, so that µ = 3.986 ×10 14 m 3 s −2 and r 0 = R earth + h, R earth = 6.371 × 10 6 m.

Simulations for Ideal Communication Channel
The simulations were performed in the MATLAB/Simulink software environment. For the numerical solution of the differential equations the MATLAB variable-step routine ode45 with the relative tolerance 10 −3 was used. For hybrid (continuous-discrete systems), including the coding/decoding procedure model, the blocks from the standard Simulink/Discrete blockset were included. The simulation time was confined to t fin = 54,000 s = 15 h.
Firstly, consider the "ideal" case where the position information is transferred over the channel without sampling and level quantization. For the performance criterion, let as pick up instant T * when the relative trajectory on the (X, Z) plane reaches the circle with given radius Q and does not leave it in the future, i.e., T * = max t x 2 12 + z 2 12 > Q (T * is further called the regulation time). It is natural to assume that within this circle the regulation rule switches to the rule ensuring collision avoidance, which is beyond the scope of this paper.
The simulation results for various initial conditions are depicted in The simulation results show that, despite of the control signal saturation at the beginning of the process, the targeting manifold is attained for both sets of the initial conditions with regulation time T * = 3.27 h (Figures 2 and 3) and T * = 4.61 h (Figures 4 and 5), respectively. Moreover, based on the harmonic balance method arguments [52,53] the conclusion can be made that the closed-loop system is globally asymptotically stable.   (5), and embedded state estimation (17). For the numerical study, the coding-decoding procedure parameters are picked up as follows: M 0 = 1, m c = 0, ρ = e −0.01T 0 , gains l 1 , l 2 of observer (17) are found by the pole-placement technique for the discrete-time systemŷ [

Case of Erasure Communication Channel
Thirdly, let us study how the erasure of data in the communication channel affects the data transmission accuracy and the regulation time for the relative position of the satellites.
For evaluating the data transmission accuracy, the transients of the transmission procedure were excluded by consideration only the last 0.3t fin interval of the simulation time. For this interval, the standard deviations σ e x 1 and σ eẋ 1 of the data transmission errors e x 1 [k] = x 1 (t k ) −x 1 (t k ), and eẋ 1 [k] =ẋ 1 (t k ) −ẋ 1 (t k ) (respectively) were calculated. Logarithmically scaled functions σ e x 1 , σ eẋ 1 versus transmission rate R for various values of erasure probability p are plotted in Figures 8 and 9. The plots show that the data transmission errors decrease monotonically, looking like inversely proportional functions on communication rate R and are practically negligibly small as R > 1 bit/s for all considered values of p.
A summary graph of the dependence of regulation time T * on the data transmission rate R (for each one channel) at various erasure probabilities p ∈ {0, 0.1, 0.2, 0.3} is shown in Figure 10. The curves in Figure 10 make an impression about the required load of the communication channel and the quality of the stabilization process with its use. It is seen from the plot that for significantly high data transmission rate (exceeding 2 bit/s in our example) erasure of data in the channel with probability up to 0.3 does not have an effect on the regulation time. This time is defined by the system dynamics regardless of the communication channel capacity. To illustrate the system performance, the particular processes for T 0 = 0.667 s, x 12 (0) = 200 m, x 12 (0) = 0.025 m/s, z 12 (0) = −50 m, andż 12 (0) = −0.025 m/s and probability p = 0.2 of erasing data in the communication channel are plotted in Figures 11 and 12. Regulation time T * = 7.07 h is found by the simulation. The time histories and trajectories for the ideal and limited capacity erasure communication channels significantly differ from the case of the ideal channel, and from an application viewpoint, the process quality for these conditions is a boundary one. .
Illustrations of the adaptive coding procedure are given in Figures 13 and 14. Adaptively tuning quantizer range M[k] in accordance with algorithm (5) is depicted in Figure 13. The plot shows how the range is automatically increased at the initial stage of the process and decays then, which leads to reducing the data transmission error. Figure 14 illustrates an influence of data erasure on the codewords, transmitted over the channel. Signal s[k] from the coder is received in the form of s ζ [k] on the decoder side. The difference between these signals causes the additional data transmission errors.

Conclusions
In the paper the feedback control law was designed ensuring the regulation of the relative satellites motion in a swarm. Unlike previous papers, the limited data transmission rate over the communication channels was taken into account. The adaptive coding/decoding procedure for the transmission of position between the satellites in the formation, employing the kinematics process description was studied for the cases of the ideal and erasure channel. Note that an adaptive coder used in the paper was not new since it was employed earlier for other problems in [77,78]. However such a coder was not applied to the control of satellite swarms previously and this is an additional novelty in the paper.
The dependence of the closed-loop system performance and accuracy of the data transmission algorithm on the data transmission rate was numerically evaluated. It was shown that both data transmission error and regulation time depend approximately inversely proportionally on the communication rate.
In future research, disturbances and noise in the communication channel will be taken into account.

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

Abbreviations
The following abbreviations are used in this manuscript: