Adaptive Coding in Control of Two Satellites Relative Motion Over the Packet Erasure Communication Channel

The paper deals with the navigation data exchange between two satellites moving in a swarm. It is focused on reduction of the inter-satellite demanded communication channel capacity taking into account dynamics of the satellites relative motion and possible erasing in the channel the navigation data. The feedback control law is designed ensuring 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. Dependence of the regulation time on the data transmission rate is numerically evaluated. 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 based on the extensive simulations. These results provide dependence of the required load of the communication channel on the desired quality of the regulation process. It is shown that both data transmission error and regulation time depend inversely proportional on the communication rate, and that for the significantly high data transmission rate erasure of data in the channel with probability up to 0.3 does not make an effect on the regulation time.


Introduction
In the 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 the 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 one, whilst the relative orbits of the followers are considered to be the Projected Circular Orbit, which is the relative orbit between the master and slave satellites, being a circular one when it is projected onto the local horizontal plane [19]. PCO predicted distances [20] between master and slave satellites are constant. The HCW equation describes the relative motion between the satellites in the local reference frame LVLH, with the origin on the master satellite. Monakhova and Ivanov [21] considered the problem of constructing a swarm from a large number of nanosatellites immediately after their separation from the launcher. To solve the problem of eliminating the drift between satellites, in [21] the decentralized control based on the force of aerodynamic drag is employed.
focus of the paper is in the navigation data exchange between the satellites to be used to keep the satellites 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 Sec. 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, various coding/decoding schemes are described. Section 3 is devoted to dynamics of two satellites relative motion in a near-circular orbit. The main result is concentrated in Sec. 4. This section is started with designing the control law, which ensures asymptotic regulation of the satellites relative motion (subsection 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, cf. [36]. 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 transmission of position between the satellites in the formation, employing the kinematics process description, is introduced in subsection 4.2. It is worth to mention that the application of this procedure makes it possible to avoid measuring of 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 subsection 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 dependance on the communication channel bitrate and erasure probability obtained based on the extensive simulations are are presented in subsection 4.4. Concluding remarks and the future work intentions given in Sec. 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 otherwise ideal. The cases of packet erasure channel and 'blinking' channel are widely appear in various real-world systems, see, e.g. [40][41][42][43][44][45][46][47][48][49]. Therefore, for 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. Rigorous examination of the influence of time sampling and the level quantization is a complex nonlinear analysis problem, that usually has not an exact analytical solution. In the early studies, the level quantization in digital control systems was usually considered as 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 the LTI plants. However, if the quantization level is relatively high (for example, in the case of binary quantization), this can lead to 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 effect 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 the 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 ∈ R m are the state, output and control vectors, respectively; A, B, C are the matrices of the corresponding dimensions; k ∈ Z + denotes the step number (the discrete time). It is assumed that pair (A, B) is reachable, (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 inequality [31]:R where η j are the eigenvalues of matrix A , j = 1, . . . , n. The right-hand side of (2), denoted asR NE = ∑ η j ≥ρ log 2 η j ρ (the Nair-Evans-, or NE-number), gives a tight admissible bound when the ρ-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 the several reasons the data bitrate is usually much greater than R NE , but NE-number can serve as a measure showing the maximum available possibilities of estimation and control over the existing communication channel. Also a promising approach is to employ the event-triggered control instead of 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 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, 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-79], where 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; ϕ(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 estimatê y(t) quantized with given M[k] as follows: where the estimateŷ(t) is generated by the following observeṙx 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; continuous-time observation errorσ (t) is found as an extension of σ [k] over the sampling interval. In the case of the zero-order extrapolation,

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 Hill-Clohessy-Witshire (HCW) model, see [21,[80][81][82][83][84] is used.
In the present paper, the LVLH reference system is employed, where OZ axis is directed from the center of the Earth, OY axis is directed normal to the orbital plane, OX axis complements the triple to the right coordinate frame, see Fig. 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 hold, 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 OX axis, which is controlled by the satellite turning through attack angle α. Based on 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 satisfyes the expression ω = µ/a 3 where µ = GM, G is the gravitational constant, M is the mass of the central body (for the Earth, µ = 398603×10 9 m 3 s −2 ), a is the semi-major axis of the satellite orbit. HCW equations (9), (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; 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 a 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 the oscillator. To resolve this contradiction, Leonard [17] has 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 ; 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 i j ≤ 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 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 forṁ where matrices A, B, K are of form (13).
The design problem consists in 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 the 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
Application of the 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]): , respectively, l 1 , l 2 are the observer design parameters.

Erasure channel description
By 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. Also 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 f i = − ρV 2 2 C α (∆S sin α i − S 0 )m −1 ms −2 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; m is the mass of the satellite. Therefore the Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 November 2020 differential drag can be found as f i j = − ρV 2 2 C α (∆S(sin α j − sin α i ), which dives 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: 24 kg, so that µ = 3.986×10 14 m 3 s −2 ; 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 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 = 54000 s=15 hours.
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 naturally 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 Figs  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 hours (Figs. 2, 3) and T * = 4.61 hours (Figs. 4, 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.

Satellite Position Transmission Over the Digital Communication Channel with Finite Capacity
Secondly, consider the case when the satellites exchange each other their LVLH coordinates over the binary communication channels, described in Secs. 2.2, 4.2. It is assumed that each satellite transmits its coordinates x i , z i , i = 1, 2 to the accompanying one, where the estimatesx i [k],ẑ i [k] andx i [k],ẑ i [k] are obtained. The following quantization and data transmission algorithms are used: binary quantization (4), adaptive zooming (5), 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ŷ [ depending on the transmitted coordinate (x or z) and the satellite number i = 1, 2. Initial conditions y[0], V [0] are set to zero for all four channels. Therefore, two data transmission channels are supposed to be used on each satellite, so the data transmission rate for each binary communication channel R = 1/T 0 should be doubled to find the overall channel load.
The simulation results for T 0 = 0.667 s (R = 1.5 bit/s for each coordinate x i and z i , transmitted from i-th satellite) and initial conditions

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 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   data transmission rate R (for each one channel) at various erasure probabilities p ∈ {0, 0.1, 0.2, 0.3} is shown in Fig. 10. The curves in Fig. 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 make an effect on the regulation time. This time is defined by the system dynamics regardless the communication channel capacity. To illustrate the system performance, the particular processes for T 0 = 0.667 s, x 12 (0) = 200 m,ẋ 12 (0) = 0.025 m/s, z 12 (0) = −50 m,ż 12 (0) = −0.025 m/s and probability p = 0.2 of erasing data in the communication channel are plotted in Figs. 11, 12. Regulation time T * = 7.07 hours is found by the simulation. The time histories and trajectories for the ideal and the 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 Figs. 13, 14. Adaptively tuning quantizer range M[k] in accordance with algorithm (5) is depicted in Fig. 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 is designed ensuring regulation of the relative satellites motion in a swarm. Unlike the previous papers the limited data transmission rate over the communication channels is taken into account. The adaptive coding/decoding procedure for transmission of position between the satellites in the formation, employing the kinematics process description is proposed and studied for the cases of ideal and erasure channel. Dependence of the closed-loop system performance and accuracy of the data transmission algorithm on the data transmission rate is numerically evaluated. It is shown that both data transmission error and regulation time depend inversely proportional on the communication rate. In the future research it is planned to take into account disturbances and noise in the communication channel.

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

Abbreviations
The following abbreviations are used in this manuscript: