Assessment of Satellite Contacts Using Predictive Algorithms for Autonomous Satellite Networks

Upcoming Low Earth Orbit Satellite Networks will provide low-latency and high downlink capacity necessary for future broadband communications and Earth Observation missions. This architecture was proposed at the beginning of the 90’s, although it has just recently re-gained popularity thanks to the so-called Mega-Constellations. This network is composed of satellites that have Inter-Satellite Links (ISL) to communicate between them. Due to the satellite motion, an ISL is a temporal contact between two satellites characterized by a lifetime in which the communication remains feasible. The determination of a route between distant satellites is a challenging problem in this context. However, the satellite follows a well-known deterministic orbit trajectory, being feasible the prediction of its position by propagating a trajectory model over time. The Contact Graph Routing protocol uses this feature to determine the evolution of the routes by pre-computing on-ground a planning of the satellite contacts. This centralized ground-dependent solution cannot be directly applied in the Internet of Satellites paradigm, which proposes the autonomous deployment of heterogeneous satellite networks without pre-assuming any specific satellite system architecture. Following this concept, the present work proposes a distributed algorithm by which a satellite predicts neighbor contacts, and generates a global contact plan without trajectory propagation. To achieve this solution, an ISL has been modeled as a “close approach” between two satellites, which is characterized by their relative motion. The present work details the predictive algorithm, and evaluates its performance in two scenarios with a hybrid satellite constellation and a mega-constellation.


I. INTRODUCTION
The Horizon 2020 Operational Network of Individual Observation Nodes (ONION) project [1] identified the current Earth Observation (EO) market needs and requirements. In particular, new satellite systems are requested to monitor weather disasters, and sea-ice level in different regions [2]. These new solutions have to provide high spatial resolution The associate editor coordinating the review of this manuscript and approving it for publication was Datong Liu . measurements as well as fast access to this data, as authors in [3] state. In the last years, Distributed Satellite System (DSS) architectures [4] have become a promising solution to satisfy these stringent requirements, by distributing a common goal among different satellites. Low Earth Orbit (LEO) Satellite Networks are DSS in which the satellites have means to communicate among them, using point-to-point links called Inter-Satellite Links (ISL). Due to satellite dynamics, an ISL is characterized by a lifetime in which the communication remains active. This temporal satellite contact makes the route definition (i.e. the routing) in the network a considerable challenge. Authors in [5] have addressed it by designing a satellite constellation as a static mesh topology, mitigating the satellite mobility impact. This approach directly fixes the satellite constellation architecture, limiting the number of possible applications.
Alternatively, other researchers have proposed heterogeneous systems to satisfy the requirements, such as the Federated Satellite System (FSS) concept [6]. This one has emerged as a promising solution in which satellites decide to sporadically establish a win-win collaboration to share unused resources, such as memory storage or downlink opportunities. This interaction between satellites, called federation, is generated opportunistically depending on the resource availability. The federation is established using an ISL, which is limited by its lifetime and effective communications range. For this reason, the Internet of Satellites (IoSat) paradigm [7] extends the concept to a multi-hop scenario in which a set of satellites decide to deploy an ad-hoc network, while they are performing their own mission. This network, called an Inter-Satellite Network (ISN), is opportunistically deployed over heterogeneous satellites with different orbit trajectories, resource capabilities, and mission goals. Therefore, an ISN is composed of satellites that do not conform a predefined constellation or architecture, instead, it represents a hybrid satellite system. For this reason, routing protocols conceived for LEO Satellite Networks cannot be directly applied in this context. Previous work in [8] presented a compatibility assessment of Mobile Ad-hoc Networks (MANET) [9] protocols to deploy ISNs. However, the performance of these solutions is impacted by the network disruption. This behavior is commonly found in satellite networks, and it represents the fragmentation of the network in different sets of satellites, which are connected between them. Due to satellite motion, this fragmentation evolves over time, making a satellite not always reachable.
This behavior has been investigated in Delay/Disruptive Tolerant Networks [10], proposing a solution that takes advantage of the satellite orbit trajectory determinism. This is the case of the Contact Graph Routing (CGR) protocol [11], [12] which determines a route between a source and a destination over time using a contact plan. This plan is a scheduled sequence of contacts between all the satellites that compose the system. The CGR protocol was designed to establish communications in the NASA Interplanetary Overlay Network (ION) [13], but it has gained interest for LEO systems. The presented protocol standard draft does not clarify however the procedure to generate the contact plan.
Traditionally, a centralized solution in which the ground facilities propagate orbit trajectories and compute the corresponding plan is used. The concept of the Contact Plan Computation Element (CPCE) to support the generation of the contact plan and the route definiton is presented in [14]. This CPCE pre-computes the contact plan and the routes onground to periodically upload them into the satellites using the ground station network. An example of a protocol that uses this centralized architecture is the Predictable Link-State Routing (PLSR) protocol [15]. Although this centralized approach can generate accurate contact plans for complex systems, the upload process of the plan needs still to be further investigated. In particular, if a new satellite is included in the system, the recomputation and upload of the entire contact plan must be conducted, making the solution difficult to be scalable. Moreover, this centralized facility entails an important operations cost, and makes the satellite system grounddependent. Due to these features, the integration of this architecture in the IoSat paradigm, where heterogeneous systems with satellites from different operators interact, becomes not suitable. Additionally, this centralized solution would require high resources for massive satellite constellations with thousands of satellites, such as StarLink [16] or Telesat [17].
To overcome these limitations, new investigations aim at providing autonomy to satellite systems by enabling the capability to identify route changes. Authors in [18] propose a routing protocol that selects routes according to its stability, characterized by its Route Expiration Time (RET). This RET is determined by the minimum lifetime of the links that compose a route. This minimum lifetime of the links is also called Link Expiration Time (LET). Featuring the protocol with predictive capabilities reduces the overhead as well as mitigates the disruption. Although a model to compute the LET is presented, its accuracy and the estimation procedure is not detailed. Additionally, the LET is determined by propagating the orbit trajectory of the satellites over time, which means more processing capacity for the satellites.
Additionally, other prediction mechanisms have been proposed from the DTN research community. The Opportunistic CGR (OCGR) [19] proposes that each node estimates contacts and their duration using a history log. This log accumulates information of contacts that have already achieved, which is updated according to the confidence parameter. This parameter represents how reliable is the new predicted contact with a satellite, according to its history log. This mechanism is characterized by having an important learning curve, which as compared to others that use determinism to estimate contacts (like CGR) degrades its performance. The amount of data that can be exchanged during a contact, well-known as contact capacity, is another parameter used by the CGR to determine the routes. Authors in [20] propose a mechanism to estimate this capacity between satellites and ground stations, modelling the features of the transceiver and the communications loss.
As opposed to the previous proposals, this work presents an autonomous and de-centralized protocol by which a satellite identifies other neighbor satellites, and predicts when a contact happens assuming an omnidirectional antenna pattern. The proposed solution uses a predictor which estimates the contacts benefiting from the determinism of the satellite trajectory. Therefore, it uses the satellite orbital elements to formulate simple model, without propagating the entire trajectory. To define the prediction algorithm, a contact between two satellites has been modeled as a satellite proximity or closeness. This scenario has largely been studied for satellite or debris collisions [21], [22], defining the Probability of Close Approach (PCA) of two satellites. As this work demonstrates, a predictor based on the PCA is impacted by the initial satellite conditions, scheduling wrong contacts. For that reason, an additional predictor is defined using the relative orbit motion theory [23], which has normally been used in flight formation and spacecraft rendezvous. With this relative orbit motion, the position of a deputy satellite with respect to a chief satellite can be determined and modelled. To the best of our knowledge, these models have not been used to estimate the time interval of a satellite contact. Note that the contact estimation is achieved considering a communications range, which is determined by the link budget of the satellites. This work also evaluates the performance of the proposed protocol in two different scenarios. The former with an hybrid satellite system composed of tropical and polar constellations that demonstrates the correctness of the contact plan construction procedure. The last scenario evaluates the time of the satellites to receive the contact plan from all the other satellites in a massive satellite constellation.
The main contribution of the presented work is the extension of [24] by providing 1) a predictor model based on close approach theory, 2) an alternative predictor definition based on relative orbit dynamics, 3) an accuracy assessment of the different predictors, 4) a protocol with which the satellites can autonomously generate a contact plan, 5) a feasibility analysis executing the protocol in a hybrid satellite constellation, and 6) an evaluation of the time to construct the contact plan in a Mega-constellation.
The remaining of the article is structured as follows. First, Section II presents a global view of the predictive algorithm. Characteristics of the predictor based on the PCA are detailed in Section III. Section IV formulates the predictor based on relative orbit motion. The mechanism to generate the contact plan and the corresponding performance analysis is presented in Section V. Finally, Section VI presents the conclusions of the work.

II. PREDICTIVE ALGORITHM OVERVIEW
In a satellite network, the nodes are constantly moving, which modifies the topology over time. However, during a lapse of time the topology remains stable generating a sample called snapshot. The topology evolution is represented thus as a sequence of snapshots, not evenly spaced in time. Due to the periodicity of satellite motion, this sequence becomes periodic repeating the same snapshots defined in m time segments. Note that this periodicity is achieved in short integration periods, because secluar effects (e.g. J2). Although these effects, the sequence model is valid due to their slow impact on the orbit. Each snapshot can be modeled as a graph G, in which n nodes (i.e. satellites) denoted by V = {v 1 ,v 2 , · · · ,v n } are connected by a set of edges E. An edge (v j ,v k ) ∈ E represents an active ISL from v j to v k nodes. The graph of a snapshot is thus determined by the nodes, and the different edges: G = (V , E). In this dynamic environment, the definition of a route between distant nodes is a considerable challenge. Authors in [25] propose a solution that computes feasible routes in each snapshot graph independently.
Furthermore, a route can also exist between snapshots, being possible to define routes in the time-domain. For instance, Fig. 1 presents a snapshot sequence (s 0 , s 1 , and s 2 ) over which node 1 can establish a route to node 5 through node 2. The space-time graph model [26] presents a succession of graphs in the snapshot sequence {G(t i )|i = 1, · · · , m}. Two types of edges cohabit in this model. The former is the temporal edge (v j ,v j ) t i that connects the same node v j between snapshots at t i and t i+1 , storing any possible data. The other edge is the spatial one (v j ,v k ) t i (where j = k) that connects different nodes in the time segment [t i , t i+1 ). Using both models it is possible to define the routes in space and time domains.
The CGR protocol, as part of the NASA ION [13], implements an algorithm based on these graph models. In particular, it defines routes over graphs using a contact plan. A contact between satellites is represented as a time interval [t j , t k ) in which an ISL remains feasible. The schedule of all the contacts between satellites conforms the contact plan, refered in this work as the global contact plan. This protocol assumes that each satellite holds this plan, which normally is pre-computed and distributed through a centralized ground entity (i.e. the operations center). Although this strategy is suitable for planned and scheduled networks, it cannot be applied to heterogeneous and autonomous scenarios, such as the IoSat paradigm [7]. To overcome this limitation, this work presents a predictive algorithm that is executed in each satellite to self-learn and self-construct the contact plan. As it is presented in the following sections, this de-centralized approach is suitable for heterogeneous satellite systems (e.g. hybrid satellite constellations) or massive satellite constellations. In order to understand the proposed algorithm, different definitions need to be presented first.
Definition 1 (Contact Sequence): Considering a node v j has a sub-set of nodes that are its neighbours N v j ⊂ V , the contact sequence of node v j with a neighbor v k ∈ N v j is the schedule of contacts with this neighbor Definition 2 (Local Contact Plan): Considering a node v j , the union of its contact sequences represents the local contact plan (i.e. only with its neighbors N v j ).
Definition 3 (Global Contact Plan): Considering the nodes that compose the entire network v j ∈ V , the union of their local contact plans C v j represents the global contact plan C.
The proposed algorithm is based on periodically predicting the different contact sequences to generate a local contact plan which is shared through the network to construct a global contact plan. Fig. 2 presents a block diagram of the modules that compose the algorithm executed in each satellite. The former component is the predictor which iteratively computes the local contact plan of the satellite from time t m until time t m+1 (where T = t m+1 −t m is the prediction period) using the orbital elements vector e = [a, e, i, ω, , M ], where a is the semimajor axis, e is the eccentricity, i is the inclination, ω is the argument of periapsis, the argument of the ascending node, and M the mean anomaly of the satellite position. 1 The goal of this predictor is to estimate future contacts without propagating an orbit trajectory model, just with a simple and direct formulation. With the different estimated contacts, the predictor constructs the local contact plan. After computing the local contact plan, the global contact plan is updated in each satellite. As this global plan is uncompleted, the satellite exchanges it with its neighbors (and vice-versa) 1 The position of a satellite can be defined by the mean anomaly angle M , the eccentric anomaly angle E, the true anomaly angle ν v, and the distance from the orbit center r. Further information of those anomaly angles is presented in Section III, and in [27]  when a contact starts. At the end, the satellite is able to have a complete contact plan, which can be used by a routing protocol (e.g. the CGR protocol).
In this solution, the definition of the predictor becomes crucial to ensure an accurate contact identification. For this reason, two different predictor models have been evaluated in Sections III and IV respectively. This predictor must be accurate enough to correctly schedule the estimated contacts. For instance, Fig. 3 illustrates two predicted contacts which both are partially correct with respect to the real contacts. In particular, the striped areas are the time slots in which the prediction and the reality match. This predictor has unsuccessfully estimated the beginning of the former estimation, and the contact duration of the other one. This situation could exist depending on the predictor definition. If the predictor is based on constantly and periodically propagate an orbit trajectory (traditional method), the error between predicted and real contacts would not exist (consuming more resources). However, if the predictor uses a relative trajectory model due to its definition, errors in the predictions can appear, being the predictor not accurate enough. Therefore, two features characterize the accuracy of a contact predictor: the punctuality, and the duration matching.
Definition 4 (Punctuality of an Estimated Contact): It is the probability to estimate a contact inside the real time slot. Considering J contacts, the end timet e,j and the start timê t s,j of the estimated contact j ∈ J (wheret s,j <t e,j ), and the end time t e,j and start time t s,j of the real contact j ∈ J (where t s,j < t e,j ), the punctuality of the jth estimated contact P j is defined as follows.
ift s,j < t s,j and t s,j <t e,j < t e,j t e,j −t s,ĵ t e,j −t s,j if t s,j <t s,j < t e,j andt e,j > t e,j 0 ift s,j ≥ t e,j ort e,j ≤ t s,j 1 otherwise (4) Definition 5 (Predictor Punctuality): Considering J contacts, the punctuality P is the average of the punctuality of all the estimated contacts P j (∀j ∈ J ). Its definition is VOLUME 8, 2020 presented in (5).
Definition 6 (Duration Matching of an Estimated Contact): It is the probability to correctly estimate the contact duration. Considering J contacts, the start timet s,j and the end timet e,j of the estimated contact j ∈ J (wheret s,j <t e,j ), the start time t s,j and the end time t e,j of the real contact j ∈ J (where t s,j < t e,j ), the duration matching of the jth estimated contact L j is defined as follows. (7).

Definition 7 (Predictor Duration Matching): Considering J contacts, the duration matching L is the average of the duration matching of all the estimated contacts L j (∀j ∈ J ). Its definition is presented in
Definition 8 (Predictor Accuracy): The accuracy of a predictor Q is the combination of the punctuality P, and the duration matching L. Q = P · L (8) Using these definitions, the following sections present the details of the predictors and their accuracy. The predictors have been defined for satellite constellations targeting EO missions. This kind of architecture is characterized by circular orbits at the similar altitudes.

III. PREDICTION WITH CLOSE APPROACH PROBABILITY
A satellite contact can be modeled as a temporal closeness of two satellites. In the course of the past years, space object proximity has concerned the community due to the possibility to experience spacecraft or space debris collisons. Consequently, large efforts intended to statistically appraise these proximities. One of the outcomes of this research has been the definition of the Probability of Close Approach (PCA) [21] between a deputy satellite with respect to a chief one. The PCA estimates the possibility of two satellites to be close enough to consider it a nearness, i.e. within a distance threshold. The corresponding formulation of the PCA is detailed in [21], provinding the following compact definition: where v c corresponds to the true anomaly of the chief satellite, M d the mean anomaly of the deputy satellite, and M c the mean anomaly of chief satellite. To better understand the PCA definition, let's inspect each term that composes its definition. The derivative term of the mean anomaly dM c (v c ) dv c characterizes the deviation of the mean anomaly M with respect to the true anomaly v due to the orbit eccentricity e: The other term M d (v c ) corresponds to the mean anomaly regions in which the distance between the deputy satellite and the chief one is less than a distance threshold, so-called close approach regions. This term can be represented by the union of each jth close approach region in the deputy satellite orbit, defined by the true anomalies v d/j1 and v d/j2 : where J is the amount of close approach regions. The method of how to compute these boundaries and regions is highly detailed in [21]. EO and broadband communications missions are normally composed of satellites that follow circular orbits (i.e. e ≈ 0). In these scenarios β ≈ 1, and the derivative term of the mean anomaly dM (v) dv does not drive the PCA behavior (i.e. dM (v) dv ≈ 1). Therefore, the PCA is defined as follows.
This continuous formulation cannot be implemented in a computational algorithm, being necessary to represent it in the discrete domain. Considering the integral step v which fragments the space  probability in each region p: where M dk (v c ) corresponds to the mean anomaly difference in the kth region of the chief true anomaly space. Equation (13) highlights the importance of the integral step v, which drives the resulting probability of the kth region (i.e. p k ). Therefore, this work explores for different steps the behavior of this probability between two satellites following orbits with 90 • of inclination (i), and periapsis argument (w) difference (Fig. 4). Potential closeness true anomaly regions are identified in all the three cases (i.e. p k = 0). However, while the resolution improves decreasing v, the probability inside the potential region does not remain homogeneous. This makes some regions more favorable to have a satellite proximity, presented in Figure 4 with blue areas. It is thus clearly possible to identify and define these potential regions using this probability. For instance, the close approach regions, in terms of true anomaly, in Fig. 4c  Following this conduct, a selection algorithm has been conceived to provide these true anomaly regions in which a contact is probable. Fig. 5 presents an example of how this algorithm works. It starts by computing the probability with v = 2π, which preliminary evaluates if two satellites could have a closeness (having a non-zero probability). If this condition is not achieved, the algorithm stops indicating that no contact is possible between the satellites. Otherwise, the integration step is decreased, and the probability is computed in each of the new regions. If the resulting probability is zero, the region is discarded. Contrarily, if the region can have a contact, it is explored in the next iteration. This procedure is repeated until reaching the minimum integral step value (configurable). At the end, the true anomaly boundaries are determined with a resolution of the last integral step used.
After retrieving the true anomaly boundaries, the start and ending times of contacts are computed with an orbit trajectory  model. This model presents the formulation to characterize the position of a satellite at each instant time. In our case, the Keplerian orbit [27] has been used as the trajectory model, because of its low-computation requirements and enough precision for mission analyses. This orbit model presents the satellite motion considering only the attraction of the Earth body, neglecting any orbital perturbation. Using Kepler's formulation, the true anomaly v, eccentric anomaly E, and mean anomaly M are defined according to the time evolution.
where M 0 represents the initial mean anomaly, e the orbit eccentricity, and n the mean angular motion.
To evaluate the orbital element difference impact on the predictions, Fig. 6 presents the estimation of the contact time interval applying the selection algorithm between two satellites with just 30 • of inclination difference. The blue line corresponds to the evolution of the distance between the satellites, and the orange one is the distance threshold d th (in this case, 2500 km). The propagation of the Kepler orbit (presented previously) in the Earth-Centered Intertial (ECI) frame [27] has been used to compute just the reference distance between the satellites (blue line). 2 When the distance is less (or equal) to the threshold value, a contact between the satellites is achieved. Using the presented predictive algorithm, the contacts are estimated and scheduled with blue dots and black dots, indicating the contact start and end respectively. The results indicate that, in this case, this solution estimates the contacts with an accuracy of 98.9%.
Extending the previous results, Table 1 presents the performance of the proposed predictor in scenarios with different orbits. The search algorithm has been executed until achieving a resolution of π 10 4 in the integral step v. This predictor is capable to correctly estimate the contact duration of all the contacts (i.e. L ≈ 100%), due to the high resolution of the selection algorithm. However, this high resolution provokes more iterations, increasing the processing time. Additionally, the punctuality of the predictor P fluctuates more, from 10.8% to 97.7%. This predictor can thus correctly estimate the duration of each contact, using more computation, but it cannot always correctly schedule them. Moreover, in cases 16, and 17 the predictor schedules perfectly the contacts, because it is able to detect that the satellites are in constant line-of-sight. Cases 5, 10, 12, 15, and 18 are scenarios in which satellites would not have any contact. The predictor cannot deal with this situation, scheduling unrealistic contacts that would never happend. For that reason, the accuracy Q in these scenarios is zero. This behavior appears because of the PCA definition itself, which does not consider the initial distance between the satellites. This initial state is driven by the argument of periapsis w, the argument of the ascending node , and the initial true anomaly M . This algorithm has also been evaluated between three NASA EO satellites: Calipso, Terra, and Aqua. As in previous cases, when the orbit elements difference increases the punctuality of the predictor is degraded. However, in case 19 the accuracy is well enough to be applied in this scenario.
Concluding this section, the predictor based on the PCA cannot always correctly estimate satellite contacts.
Additionally, due to the search algorithm, it presents an iterative execution which depends on the integration step, that could be computational expensive. For these features, this predictor is not suitable to generate a reliable contact plan. An alternative predictor model has been considered in the section that follows, which presents better features and performance.

IV. PREDICTION WITH RELATIVE ORBIT MOTION
During the last years, different researchers have focused on conceiving flight formation and rendevouz missions. For these missions, the relative position of a deputy satellite with respect to a chief satellite has been studied generating different relative orbit models. Those models are based on the Hill coordinate frame O [28], which is a cartesian coordinate system centered in the chief position. The orientation of this system is given by the unit vector {ô r ,ô θ ,ô h }, whereô r is in the orbit radius direction, whereasô θ is parallel to the orbit momentum vector, andô h completes the right-handed coordinate system. The relative position vector ρ of a deputy satellite with respect to the chief one is expressed on three components.
This vector represents thus the position of the deputy satellite from the point of view of the chief, determining directly the distance between them.
The definition of these three components becomes key to determine the distance. Different researchers have used nonlinear models to characterize these components. However, in [23] Schaub presents linear expressions using the orbit element differences δe [29].
This representation simplifies the definition because only the δM component is time variant. Applying (17), the Clohessy-Wiltshire (CW) relative equations [30], which determine the relative motion with respect to a circular chief satellite orbit, have the convenient analytic solution.
x(t) = A 0 cos(nt + α) + x off (18a) where A 0 , B 0 , α, β, x off , and y off are the integration constants determined by the initial conditions, n the mean angular motion of the chief satellite, and t represents the time. These expressions present the relative motion of a deputy satellite with respect to a chief satellite, while this chief follows a circular orbit. EO and broadband communications missions are normally composed of satellites that follow circular orbit, and same altitude (i.e. e ≈ 0, δe ≈ 0, and δa ≈ 0). Applying these conditions into the definition of the integration constants presented in [23], they become x off = 0, (19e) where M 0 corresponds to the initial mean anomaly of the chief. As noted previously, this model has largely been used to estimate the distance between satellites. To the best of our knowledge, it has never been used before to estimate the time slot of a satellite contact. For this reason, the following work describes the formulation to estimate the contact time boundaries from the relative orbit motion equations. Considering an omnidirectional antenna pattern, an ISL is feasible if (and only if) the distance d is less (or equal) to the communication range d th .
d ≤ d th (20) Applying the values of (19) in (18), the ISL existence condition, presented in (20), becomes At the end, a more compact formulation of the ISL existence condition is presented in (22).
where σ = nt + β is a compressed expression of the time variant component, and γ determines ISL characteristics. Depending on γ , this condition is always or never satisfied. |d th | ≥ |y off | (23) Proof: From (22), γ can be complex C or real R. An ISL exists if γ belongs to the real set γ ∈ R. To achieve this situation, d 2 th − y 2 off ≥ 0, which concludes to (23). (24) is satisfied.

Theorem 2 (Sporadic Inter-Satellite Link): Let d th be the maximum range of an omnidirectional ISL transceiver, y off the integration constant in the in-track component of the deputy relative position, and a δ ω the integration constant in the cross-track component of the deputy relative position. An ISL with these features is time-invariant (i.e. infinite duration), if the condition
Otherwise, it becomes time-variant with a bounded lifetime. Proof: Equation (22) presents the condition in which the distance is less than a distance threshold, having a feasible ISL. The condition is determined by a cosinusoidal behavior, defined in the range [−1, 1]. For this reason, if γ ≥ 1 the condition is always satisfied, and the ISL is feasible. Taking the definition of γ , and this condition, the theorem is proven. Theorems 1 and 2 present two conditions that filters those contacts that are feasible. When the corresponding conditions are satisfied, a temporal contact between chief and deputy satellites exists. The time boundaries of this contact are given by the solution of (22). Because γ is defined by a square root,  two values exist that can satisfy (22). Let's denote γ + as the positive solution of the square root (i.e. γ + = |γ | ≥ 0) and γ − as the negative one (i.e. γ − = −γ + ). Following these definitions, (22) is presented with different formulations. Fig. 7 represents the values of a cosine in the axes and the σ space in which both conditions are satisfied. In particular, σ has to be bounded in the following space to satisfy (25a): arccos(γ + ) ≤ σ ≤ arccos(γ − ) = arccos(−γ + ) (26) Simultaneously, σ must be defined inside the following boundaries to satisfy (25b): Taking in consideration σ and the γ − definitions, the contact time is bounded by known limits, shown in (28).
where k represents the number of turns that a satellite has performed following the orbit. These boundaries determine the beginning and the end of a contact. Equations (28a) and (28b) determine the model of the predictor based on relative satellite motion, which is used to construct the global contact plan. Fig. 8 presents the execution of this predictor showing different estimated contacts between two satellites with an inclination difference of 80 • . The distance (blue line) evolves over time generating different regions in which it is lower than the distance threshold 2500 km (orange line). In this scenario, contacts are correctly estimated and scheduled, with an accuracy of 92.6%.
Extending the results, this predictor has been evaluated in scenarios with different orbits (as previously done in Table 1). Table 2 presents its performance in terms of punctuality (P), duration matching (L), and accuracy (Q). The results indicate that the accuracy of the predictor is high (i.e. larger than 90%) when the inclination is the main difference. Furthermore, the duration matching metric is the one which determines the accuracy in almost all the cases. This indicates that the predictions are scheduled correctly inside the real contact time slot, however the duration of the estimated contact does not match with the reality. This impacts the communications duration, but it ensures that in all the contacts the ISL is a fact. Additionally, the predictor is able to detect those cases in which a contact never appears, such as cases 5, 10, 12, 15, and 18. However, when the difference is based on the right ascending node, the argument of perigee, and the initial mean anomaly the prediction is less accurate. This behavior is induced by the linearization of the relative motion (i.e. the CW equations), which is valid when the orbit element difference vector is small. Fig. 9 presents an example of the accuracy improvement while using the non-linear definition  of the integration constant δ ω , presented in (29). cos(δ ω ) = cos(i) cos(i+δi)+sin(i) sin(i+δi) cos(δ ) (29) Despite this limitation, the predictor is still able to estimate contacts between real satellites. Furthermore, the predictor is able to identify those cases in which the satellites are always in line-of-sight, having a constant ISL. For instance, the accuracy of the predictor is 99.9% of cases 16, and 17. This performance is due to the inclusion of the orbital elements in the definition of the predictor, unlike the predictor based on the PCA.
In Table 2, case 19 presents the accuracy of the predictor while detecting contacts between Calipso and Terra NASA satellites. As both follow close orbits in terms of orbital elements, the predictor is able to estimate the contacts with 94.3% of accuracy. This accuracy is degraded when the prediction is computed between Terra and Aura satellites (case 20), because the orbital element difference increases. Finally, the predictor can estimate contacts with an accuracy of 22.2% between Calipso and Aura satellites (case 21). Although this reduction of the accuracy, all the estimations have 100% of punctuality (P), indicating that they are inside the real time slot (i.e. maximum punctuality).
Comparing the performance of the predictors based on the PCA and the relative satellite motion, the last one is able to estimate and schedule satellite contacts with higher accuracy.
However, the linear representation of the predictor decreases the duration matching of the estimations when the orbit element difference is considerable. Although this limitation, all the estimated contacts are punctual, scheduling them when an ISL is feasible, and thus ensuring the communications. For these features, this predictor has been selected to generate the contacts that are used to construct the local contact plan.

V. GENERATING THE CONTACT PLAN
Previous section has presented how the contact prediction is computed and the performance of the predictor. Implementing this predictor in each satellite provides the local contact plan, but not the global one. This section describes in details the complementary part of the algorithm to construct a global contact plan in each satellite.
The algorithm is based on exchanging different information between satellites. The former is the orbital elements which are required by the predictor to estimate the contacts, and to define the local contact plan. In this scenario, each satellite periodically transmits its Two Line Elements (TLE) which are received (in due course) by its neighbors. This transmission is performed using the User Datagram Protocol (UDP). After receiving a TLE, the satellite retrieves the orbital elements of the neighbor, and defines the contact sequence of this deputy satellite using the predictor. Thus, the predictor iteratively defines future contacts after finishing the last one. The number of contacts predicted per iteration can be configured depending on each scenario. Note that a tradeoff regarding the forward knowledge with respect to the memory usage needs to be conducted in order to define this parameter. This tradeoff has not been considered in the present work because it depends on the spacecraft characteristics, related to each mission.
With the accumulation of the different TLEs, the satellite defines the local contact plan. Then, the satellite schedules a request of the local contact plan to each neighbor. To do that, it establishes a Transport Control Protocol (TCP) connection with the neighbor satellite, and it iteratively updates the internal global contact plan with the received local contact plan. Note that the local contact plan of the neighbor is at the same time incomplete, because it only contains information that it learned. However, the global contact plan is being completed after a lapse of time.
This de-centralized solution becomes suitable for heterogeneous satellite systems, like the ones proposed in the IoSat paradigm. In this scenario, the satellites are operated by different operation entities, and follow different orbit trajectories. Subsection V-A evaluates the performance of the predictive algorithm with a hybrid and heterogeneous satellite system. Additionally, the de-centralized solution provides a scalable strategy that is suitable with large satellite constellations, such as the Mega-constellations. For that reason, the performance in terms of construction time has been evaluated in Subsection V-B. Both scenarios have been executed using the event-based simuation engine conveived for DSS architectures [31]. Based on Network Simulator -vesion 3 (NS-3), this simulator provides a satellite environment (e.g. satellite motion) to evaluate the predictive algorithm.

A. SCENARIO WITH A HYBRID SATELLITE SYSTEM
The de-centralized nature of the proposed predictive algorithm make it suitable for an heterogeneous and hybrid satellite system, which could be composed of different satellite systems. As the ONION project highlights in [1], this kind of system could provide different features that benefit current EO missions. Therefore, to evaluate the capacity to construct the global contact plan in this context, a case study with an hybrid satellite system is presented in this subsection. Table 3 shows the orbit characteristics of the hybrid satellite system, composed of a tropical satellite constellation and a polar one. Five satellites belong to the polar satellite constellation, which follow an obrit trajectory with 90 • of inclination. Each satellite plane is equidistant with δ of 50 • , ensuring a total coverage of the equatorial area. Additionally, five more satellites compose the tropical constellation with an orbit of 10 • inclination. As in the previous constellation, satellite planes are equidistant placed with δ = 50 • . This satellite system is a representation of a hybrid system composed of the Cyclone Global Navigation Satellite System (CYGNSS) [32], and Planet Labs satellites [33]. Fig. 10 presents the global contact plan generated from the physical trajectory of the satellites (without using any predictive algorithm) after 14000 seconds of simulation. The local contact plan is presented for each satellite, that contains the different contacts over time with the neighbors (represented as colored boxes). Note that each satellite of the tropical constellation has a unique neighbor, which belongs to the polar satellite constellation. In particular, this neighbor is the one that follows an orbit with the same argument of the ascending node . In addition to these contacts, the satellites of the polar constellation have contacts between them in the Earth poles. This kind of constellation is characterized by having high satellite density in these areas, which is interesting to exchange information. Although the satellites of the tropical constellation cannot establish ISLs among them, they can use the contacts with the satellites of the polar constellation to recover their contact plan. This is achieved with the proposed predictive algorithm. Fig. 11 presents three moments (indicated with the red line) of the global contact plan construction from the satellite sat-0 point of view. At the very beginning, sat-0 only has information about its neighbor satellite sat-5. For that reason, only the contacts with this neighbor are predicted in Fig. 11a. When this contact is feasible, it exchanges the global contact plan with sat-5, and vice versa. Due to the polar constellation nature, sat-5 has had the contacts with the other satellites of the constellation before 2000 seconds of simulation. This enabled it to learn about its new neighbors (predicting the corresponding contacts), and update its global contact plan with their information. For that reason, when sat-5 provides its global contact plan to sat-0, it can correctly schedule the different contacts of the other satellites, as shown in Fig. 11b.
This behavior is continuously repeated completing the global view of sat-0, presented in Fig. 11c.
Note that the number of predicted contacts by the predictor performed at each iteration limits the forward knowledge of the global contact plan. Depending on this knowledge, the performance of the routing protocol (e.g. the CGR) could be impacted. However, increasing this parameter would require more memory storage in each satellite. This study needs to be conducted for each satellite mission, which is out of scope of the presented work.
In conclusion, the presented results demonstrate that using the predictive algorithm each satellite can dynamically construct the global contact plan by its own. This de-centralized algorithm is suitable for hybrid and heterogeneous satellite systems, because it does not consider any predefined satellite architecture. The following subsection evaluates the performance of the predictive algorithm in a scenario with a large number of satellites.

B. SCENARIO WITH A MASSIVE SATELLITE CONSTELLATION
Due to the satellite motion and the algorithm itself, each satellite requires a lapse of time to complete the global contact plan of the network. In the previous analysis, the complete global contact plan has been achieved in a small heterogeneous satellite system. An additional analysis with hundreds of satellites composing a massive satellite constellation has been performed. This analysis highlights the impact of network disruption and network size on the performance of the proposed solution. Table 4 presents the characteristics of the Telesat constellation, detailed in [34]. It is composed of two groups of satellites: Type 1 and Type 2. The former is composed of 6 polar planes in which 12 satellites per plane follow a circular   orbit at 1000 km of altitude. The other group follows the same architecture with planes at 1200 km of altitude and 37.4 • of inclination. All the planes are equally spaced in terms of argument of the ascending node , and the satellites are also equally located in the orbit with an incremental argument of periapsis ω. Note that this architecture is similar to the hybrid constellation presented in Table 3, with a larger number of spacecraft. In total, 122 satellites are deployed to achieve global Earth coverage.
As explained before, using the predictive algorithm a satellite learns the contacts of other satellites, meanwhile it shares its own knowledge. Thus, the satellite discovers new remote spacecraft with which interact using a scheduled sequence of contacts. Fig. 12 presents the averaged discovered remote satellites over time, considering all the satellites that compose the constellation. The number of discovered satellites depends on the connectivity of the constellation, which is directly related to the ISL communications range. For that reason, Fig. 12 shows the results with different ranges. Note that the satellites start discovering only their neighbors, presenting a slow slope in all the cases. After 1000 seconds of simulation, the contact plan of each satellite increases, sharing the predicted contacts with their neighbors. Therefore, each satellite discovers remote satellites that cannot be directly linked (i.e. point-to-point). This behavior remains until reaching an inflection point (at 4500 seconds) in which the number of discovered satellites increases substantially. This event appears in all the scenarios because each satellite collects enough information to entail a considerable benefit to the others in their contact plan construction. The importance of this event is linked to the communications range, having a maximum at 2750 km. After that, each satellite has wide knowledge of the network that the exchange of contact plans does not signify an important improvement.
Some differences are detected depending on the connectivity of the network. In particular, Fig. 12b shows an important correlation: the number of discovered satellites at the end of the simulation varies depending on the communications range. The disruption of the network is the cause of this situation, limiting the possibility to collect information from other satellites. The satellites can discover the 95% of the entire network using ranges larger than 2750 km. However, in those scenarios that the ranges are smaller, the vision of each satellite is limited. This low range provokes disruption of the network that generates sub-networks, so-called communities, in which all the nodes are connected (at least once). Fig. 13 presents this existence of communities by plotting the integration of the different topologies over time. The integrated topology with 2000 km of communications range is displayed in Fig. 13a. Note that a large group of satellites are isolated becoming small white dots in the representation. However, 10 communities (represented by colors) are generated in this scenario. These communities are generated because satellites in adjacent planes with polar orbits have close approaches in polar areas. Additionally, those satellites have a proximity with the satellites with tropical orbits that have the same argument of periapsis ω. When the range increases, a community merges with other ones resulting to a unified and larger community. This is presented in Fig. 13b, that a large blue community has been created with 2250 km of range. The predictive algorithm is able to learn about these communities, and thus generate the contact plane of them. Therefore, it can adapt itself and learn the status of the network that is physically useful to communicate for each satellite, learning what can be used.

VI. CONCLUSIONS
This work has presented a de-centralized predictive algorithm that can be used in satellites to estimate contacts with other satellites, and construct the contact plan. This contact plan is a schedule of the contacts between all the satellites that compose the system. This capability would improve current and future satellite missions by providing satellite networks. Due to satellite motion, this kind of network is impacted by disruptions or fragmentations of the topology, limiting the communications between satellites. Constructing a contact plan enables to predict link status, and thus to overcome this network phenomena. Alternatively to centralized solutions, the proposed algorithm is flexible, modular, and scalable which make it suitable for heterogeneous/hybrid satellites constellations, and massive satellite constellations (e.g. Mega-constellations).
The proposed algorithm is based on a predictor that iteratively estimates satellite contacts, which are then shared with the other satellites. Two predictors, and a mechanism to evalaute their prediction accuracy have been defined in this work. The former estimates contacts using the Probability of Close Approach (PCA) definition, largely used for space object collisions. As highlighted previously, this predictor can correctly estimate the duration of the contacts, but the schedule is not properly done (reaching accuracies of 10% in some cases). For this reason, another predictor using relative orbit motion models has been defined. This new predictor correctly estimates contacts between two satellites following similar orbits. Although it schedules the estimated contact at the correct moment (punctuality larger than 90%), the duration of the estimated contacts between satellites following orbits with important differences is not perfectly done. This situation is due to the linearity of the model, which makes it valid for certain scenarios. Therefore, future efforts may be oriented to non-linear relative motion models.
In addition, this work has presented the communication protocol necessary to exchange the contact plan between the satellites, to construct a global one. With a combination of the publication of the Two Line Elements (TLE) and the connection establishment, the system is able to iteratively construct the contact plan in each satellite. A representative hybrid constellation combining CYGNSS and Planet Labs satellites is used to evaluate the performance of the proposed solution. These results demonstrate that each satellite is able to construct the global contact plan, knowing at the end feasible routes to communicate with remote satellites. However, it is highlighted that this final knowledge is achieved after a lapse of time in which the algorithm is discovering.
To evaluate the nature of this feature, it has been evaluated the algorithm in the Telesat mega-constellation. The results demonstrate that this time is directly related to the ISL communications range, which determines the disruption of the network impact. However, an inflection point in which all the satellites have the needed knowledge appears in all the cases. Additionally, the results also highlight that the proposed algorithm is able to discover those remote nodes with which phisically an interaction can be established. Therefore, depending on the disruption of the network, the satellite does not necessarily learn about the entire network, just about its community. This efficient behavior makes it suitable for this kind of massive satellite constellations.
For all the presented work, the benefits and the feasibility of applying the proposed algorithm has been demonstrated for future satellite constellations. As the presented work has been focused on the contact plan construction, a study has to be performed regarding how this iterative process can connect with the Contact Graph Routing protocol. Additionally, due to the little level of accuracy uncertainity, this solution could be evaluated as an integration of routing protocols that take in consideration this situation, such as in [35]- [37].
ANNA CALVERAS AUGÉ was born in Barcelona, Spain, in 1969. She received the Ph.D. degree in telecommunications engineering from the Universitat Politècnica de Catalunya, Spain, in 2000. She is currently an Associate Professor with the Wireless Networks Group (WNG), Computer Networks Department. Her research interests and expertise areas comprise the design, evaluation and optimization of communications protocols and architectures for cellular, wireless multihop networks, ad hoc networks, wireless sensor networks, the Internet of Things, and application domains, such as smart cities, building automation, and emergency environments. She has been involved in several National and International research or technology transfer projects. She has published in International and National conferences and journals.
ADRIANO CAMPS (Fellow, IEEE) was born in Barcelona, Spain, in 1969. He received the degree in telecommunications engineering and the Ph.D. degree in telecommunications engineering from the Universitat Politècnica de Catalunya (UPC), Barcelona, in 1992 and 1996, respectively. In 1991 to 1992, he was with the ENS des Télécommunications de Bretagne, France, with an Erasmus Fellowship. Since 1993, he has been with the Electromagnetics and Photonics Engineering Group, Department of Signal Theory and Communications, UPC, where he was first Assistant Professor, an Associate Professor, in 1997, and a Full Professor, since 2007. His research interests include microwave remote sensing, with special emphasis in microwave interferometric radiometry by aperture synthesis (ESA's SMOS mission), passive microwave remote sensing using signals of opportunity (GNSS-reflectometry), and more recently on the use of nano-satellites as cost-effective platforms to test innovative concepts for Earth Observation, such as the 3Cat-2, a 6U-class CubeSat mission launched, in August 2016, carrying onboard an innovative GNSS-R payload. He has published more than 180 journal articles in peer-reviewed journals, and more than 350 international conference presentations, and holds ten patents.