P Ilabaca
email: pablo.ilabaca@mail.concordia.ca.

S Montejo-Sánchez

J A Fraire

R D Souza

S Céspedes

Network Size Estimation for Direct-to-Satellite IoT Pablo Ilabaca, Samuel Montejo-Sánchez, Senior Member, IEEE, Juan A. Fraire, Senior Member, IEEE, Richard Demo Souza, Senior Member, IEEE, and Sandra Céspedes, Senior Member, IEEE Abstract-The worldwide adoption of the Internet of things (IoT) depends on the massive deployment of sensor nodes and timely data collection. However, installing the required ground infrastructure in remote or inaccessible areas can be economically unattractive or unfeasible. Cost-effective nanosatellites deployed in low Earth orbits (LEO) are emerging as an alternative solution: on-board IoT gateways provide access to remote IoT devices, according to direct-to-satellite IoT (DtS-IoT) architectures. One of the main challenges of DtS-IoT is to devise communication protocols that scale to thousands of highly constrained devices served by likewise constrained orbiting gateways. In this paper, we tackle this issue by first estimating the (varying) size of the device set underneath the (mobile) nanosatellite footprint. Then, we demonstrate applicability of the estimation when used to intelligently throttle DtS-IoT access protocols. Since recent works have shown that MAC protocols improve the throughput and energy efficiency of a DtS-IoT network when a network size estimation is available, we present here a novel and computationally-efficient network size estimator in DtS-IoT: our optimistic collision information (OCI) based estimator. We evaluate OCI's effectiveness with extensive simulations of DtS-IoT scenarios. Results show that when using network size estimations, the scalability of a frame slotted Aloha-based DtS-IoT network is boosted 8-fold, serving up to 4 × 10 3 devices, without energy efficiency penalties. We also show the effectiveness of the OCI mechanism given realistic detection ratios and demonstrate its low computational cost implementation, making it a strong candidate for network estimation in DtS-IoT.

Index Terms-Direct-to-Satellite, Internet of Things, Medium Access Control, Network Size Estimation, Satellite Communications.

I. INTRODUCTION

The Internet of things (IoT) market is thriving. By 2025, the annual revenue of IoT will reach 11 trillion dollars obtained from more than 30 billion connected devices [START_REF] Greig | IoT industry projected to have economic impact of $11 trillion by 2025 -TechRepublic[END_REF], [START_REF] Vailshery | Global IoT and non-IoT connections 2010-2025 Statista[END_REF]. Applications enabled by IoT networks include monitoring the environment, meteorological or environmental phenomena, animals or plants as well as tracking and positioning services. Applications for logistics and location are also available among many others. However, many applications require the IoT devices to be deployed in remote areas with minimal or no human presence. In cases where connectivity via terrestrial techniques might be unappealing or not possible, satellite technology emerges as the appropriate approach. Indeed, satellite IoT has recently been recognized as instrumental to ensure data collection over remote regions of the planet, enabling a truly global IoT solution [START_REF] Centenaro | A Survey on Technologies, Standards and Open Challenges in Satellite IoT[END_REF].

The ad-hoc nature of space missions driving the deployment of satellite infrastructure typically demands a full reengineering of the complete space and ground segments-a factor that has kept most countries out of satellite technology development. However, the CubeSat standard provided a common mechanical and electrical framework for the development of nanosatellites [START_REF] Akyildiz | The internet of space things/cubesats[END_REF]. Since its introduction, the cost of these spacecraft has dropped dramatically, allowing new space actors to access a more "democratic" space [START_REF] Diaz | New opportunities offered by cubesats for space research in latin america: The suchai project case[END_REF]. A 1U (1 unit) CubeSat is characterized by a weight of 1 kilogram and a cubic shape of 10×10×10 centimeters. CubeSat platforms of 3U, 6U, and 12U are also popular in the industry. These homogeneous characteristics reduce manufacturing, procurement, and launching costs but, at the same time, impose limiting conditions such as lower processing, energy, and storage capabilities. Thus, CubeSat-based missions tend to be heavily constrained in terms of resources. The research presented in the remainder of this work is inspired by CubeSat platforms devoted to providing direct-to-satellite IoT (DtS-IoT) services.

In the DtS-IoT architecture, actuators and sensors on the ground communicate directly with the passing-by low Earth orbit (LEO) satellite leveraging technologies such as LoRa, long range-frequency hopping spread spectrum (LR-FHSS), narrowband-IoT (NB-IoT) and the fifth generation (5G) mobile communication [START_REF] Boquet | LR-FHSS: Overview and Performance Analysis[END_REF]- [START_REF] Leyva-Mayorga | LEO Small-Satellite Constellations for 5G and Beyond-5G Communications[END_REF]. Since a ground-based IoT gateway is not required, this architecture is particularly appealing for efficiently servicing remote areas such as oceans, mountains, or inaccessible regions, as well as disaster zones [START_REF] Fraire | Direct-to-satellite iot -a survey of the state of the art and future research perspectives[END_REF].

The main challenge of DtS-IoT is to ensure good performance in the uplink shared channel as the number of nodes in coverage may be significantly higher than in traditional terrestrial networks. Indeed, link saturation is already achieved with a few thousand nodes even when state-of-the-art wide-area IoT technologies are implemented [START_REF] Boquet | LR-FHSS: Overview and Performance Analysis[END_REF]. Beyond this scalability point, packet collisions and energy inefficiencies undermine any viable DtS-IoT application. In this context, media access control (MAC) protocols for DtS-IoT are of paramount importance in massive deployments [START_REF] Ferrer | Review and evaluation of mac protocols for satellite IOT systems using nanosatellites[END_REF]. Furthermore, the network size under the footprint of the nanosatellite is a priori unknown, but an accurate network size estimation can be exploited by MAC protocols to optimize the shared uplink channel resources [START_REF] Buratti | An overview on wireless sensor networks technology and evolution[END_REF], [START_REF] Vogelgesang | Uplink transmission probability functions for LoRa-based direct-to-satellite IoT: A case study[END_REF].

Indeed, estimating a node population's size is a major factor for achieving scalability and energy efficiency on certain access control protocols [START_REF] Vogelgesang | Uplink transmission probability functions for LoRa-based direct-to-satellite IoT: A case study[END_REF], [START_REF] Zhao | Slotted aloha game for medium access control in satellite networks[END_REF]. Additionally, knowing the number of devices to be served by a CubeSat guarantees the quality of service for their applications. Depending on how overwhelmed a CubeSat will be during operation, it is possible to justify whether to implement more satellites into the constellation or lower service capacity. Nonetheless, to the best of our knowledge, efficient network size estimation mechanisms operating in constrained IoT gateways deployed in continuously-moving LEO CubeSats have not yet been addressed in the literature. The problem of an unknown and massive population of nodes trying to communicate with nanosatellites in applications that satellite technology supports is recent. Consequently, the importance in the democratization of low-cost global connectivity and the lack of studies in the literature motivate the proposal of network size estimation techniques for this particular scenario.

In the past, several network size estimation mechanisms have been proposed for terrestrial wireless sensor networks (WSN) [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF]. These mechanisms assume static base stations and not necessarily energy-constrained devices, making their direct application difficult in the DtS-IoT context. The work in [START_REF] Ali | Two algorithms for network size estimation for master/slave ad hoc networks[END_REF] is aimed at dealing with time-evolving topologies, but the method is based on persistent node connectivity, a feature not present in DtS-IoT. Conversely, in this work we explore the design of network size estimation mechanisms specifically suited for LEO CubeSats providing DtS-IoT services to thousands of IoT nodes. Our preliminary results of the applicability of the estimator proposed in [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF] to DtS-IoT scenarios are presented and analyzed in [START_REF] Ilabaca | Network size estimation in direct-to-satellite iot[END_REF].

The contribution of this work is twofold. First, we propose the optimistic collision information (OCI)based estimator: a network size estimation mechanism for DtS-IoT networks. OCI is a novel scheme with reduced computational complexity, a high adaptability to different DtS-IoT scenarios in terms of detection ratios, and an increased numerical stability compared to existing estimators. We demonstrate that OCI achieves high scalability and a low estimation error. Second, we assess the practical value of leveraging the estimations from OCI. To this end, we present extensive simulation results of a DtS-IoT network governed by a frame slotted Aloha (FSA) protocol, where device transmissions are throttled by OCI estimations. We show that when OCI is employed, the network achieves a throughput close or equal to the theoretical maximum in DtS-IoT regions with more than 4 × 10 3 devices; furthermore, we provide theoretical and empirical complexity analyses of the network size estimation mechanisms, using a typical CubeSat Average of m estimates of the contender nodes on a frame k architecture. The proposed estimation method runs in only 24.3% of the time required by its most relevant competitor.

The remainder of this paper is organized as follows. Section II reviews the related work. Section III provides the description of the scenario of study and the design of the proposed network size estimator. Section IV analyzes the performance of the OCI mechanism and presents the results of exploiting size estimations in a DtS-IoT scenario. Section V delves into practical considerations while Section VII concludes the paper. Table I offers the definition of the different acronyms and Table II provides a list of variables used in this paper.

II. BACKGROUND

In this section, we first briefly discuss orbit characteristics for deploying space technology, and then provide a review of previous direct and indirect network size estimation mechanisms proposed in the literature.

A. Space Technology and Orbits Characteristics

Satellites can be deployed in different orbits with diverse characteristics that condition the exploitation of the hosted communication payload. The three main types of satellite orbits are the following:

• Geosynchronous Equatorial Orbit (GEO): This orbit has a rotation period equal to that of the Earth. A GEO satellite appears motionless to an observer on the ground. GEO satellites have an altitude of 35,786 km, which gives them a coverage of approximately 30% of the Earth's surface. This orbit offers services such as weather data and broadcast TV [START_REF] Ses | How orbital altitude impacts network performance in satellite data services[END_REF]. However, the long distance to the Earth surface results in high latency, on the order of 120 ms for the best case scenario [START_REF] Ferrer | Review and evaluation of mac protocols for satellite IOT systems using nanosatellites[END_REF]. Furthermore, high transmission power and large antenna gains are required [START_REF] Ses | How orbital altitude impacts network performance in satellite data services[END_REF]. Also, the process of deploying a GEO satellite is the most expensive among the three orbits. • Medium Earth Orbit (MEO): With an altitude ranging between 2,000 km and 3,586 km, MEO satellites experience a lower latency than GEO satellites and provide navigation services such as the Global Positioning System (GPS). The launching cost of MEO satellites is midway between the cost of the systems located in the other two orbits.

• Low Earth Orbit (LEO): Its altitude ranges between 300 km and 2,000 km, offering the lowest latency, on the order of tens of milliseconds [START_REF] Ferrer | Review and evaluation of mac protocols for satellite IOT systems using nanosatellites[END_REF]. Due to the smaller distance to Earth, LEO satellites require lower transmission power and launching costs. Although smaller and cheaper, LEO satellites move at ∼7 km/s, which may result in passes of <10 min over a given spot on the surface. Thus, LEO constellations are needed to provide persistent coverage. LEO satellites will play a key role in providing ubiquitous and low-cost coverage for IoT devices in the near future [START_REF] Leyva-Mayorga | LEO Small-Satellite Constellations for 5G and Beyond-5G Communications[END_REF]. One of the main issues in estimating the size of a network as seen from a satellite passing by is that the ground track (the projection of the satellite trajectory over the Earth surface) does not necessarily repeat. In other words, the satellite might observe a given set of nodes in one pass, and a subset of those (plus possible others) in subsequent passes. While such a condition challenges any network size estimation techniques, orbits that periodically cover the same ground track are often used in practical space applications (i.e., multiple observations of the same area). These orbits are known as Earth-repeat orbits [START_REF] Montenbruck | Satellite orbits: models, methods, and applications[END_REF], [START_REF] Lara | Searching for repeating ground track orbits: a systematic approach[END_REF] and are described in detail in Section V.

B. Direct Network Size Estimation Approaches

In [START_REF] Ali | Two algorithms for network size estimation for master/slave ad hoc networks[END_REF], the authors present two algorithms for estimating the size of controller/agent and ad hoc networks. These algorithms are adaptations of the existing random tour and gossip-based aggregation methods to fit controller/agent architectures and to be compatible with random changes in the topology. The random tour method is based on the transfer of a token, registering each new visited node to make the estimation. The gossip-based aggregation method is based on communication between adjacent nodes. The network is divided into clusters, with a controller in each cluster repeating a process in which it averages a value with the agent nodes, sends the updated value to all agents, and eventually obtains a converged value equal to 1 N , where N is the number of nodes. The two algorithms proposed in [START_REF] Ali | Two algorithms for network size estimation for master/slave ad hoc networks[END_REF] require communication among nodes, which is most likely not available in the DtS-IoT architecture.

In [START_REF] Bustos-Jimenez | Estimating The Size Of Peer-To-Peer Networks Using Lambert's W Function[END_REF], the authors developed a network size estimation mechanism focused on Peer-to-Peer (P2P) networks. This mechanism "marks" neighboring nodes that are discovered when notified by other nodes of the network in a Capture-Recapture (CR) logic. This approach is based on a Breadth-First Search (BFS) rooted in the node that wants to discover the size of the network. The sum of new neighbors discovered at each iteration of the BFS is fitted to a logarithmic function whose coefficients are empirically set. Then, this function is solved using Lambert's W function [START_REF] Corless | On the lambertw function[END_REF] to estimate the network size. An error of at most 10% is reached, allowing the BFS to iterate only up to the third level. Later, the algorithm of [START_REF] Bustos-Jimenez | Estimating The Size Of Peer-To-Peer Networks Using Lambert's W Function[END_REF] is modified to obtain a different logarithmic function, achieving a better performance in terms of accuracy versus the number of operations [START_REF] Bustos-Jiménez | Estimating the size of natural networks using local information[END_REF]. The improvement reaches an estimation error of at most 5% while visiting less than 10% of the network. Once again, network exploration is based on direct communication among peers, making the approaches non-feasible for DtS-IoT scenarios.

In [START_REF] Accettura | The Capture-Recapture approach for population estimation in computer networks[END_REF], a survey of the CR approach for population estimation in computer networks is offered. In general, CR solutions show low precision and the estimation process may be unsuitable for cases when the network size changes are fast, which is expected from the DtS-IoT scenario. A more reliable and faster estimate can be obtained by storing the identity of the elements (sensor nodes) caught in the different samples, but paying the price of higher storage and computational requirements to handle the sampling history, which could be difficult for a CubeSat in terms of memory. In [START_REF] Peng | Estimation of a Population Size in Large-Scale Wireless Sensor Networks[END_REF], the corner stone of the CR methods, the Lincoln-Petersen Index, is employed to estimate the size of a Wireless Sensor Network (WSN). In [START_REF] Achtzehn | IEEE 802.11 Wi-Fi access point density estimation with capture-recapture models[END_REF], this method is used to estimate the number of wireless access points in a closed study area, as well as the Jackknife method. However, both of these methods work well on a closed population, i.e. no element is entering or leaving the population during the sampling and estimate operations, which differs from the dynamic characteristic of DtS-IoT networks. In [START_REF] Accettura | The Capture-Recapture approach for population estimation in computer networks[END_REF], CR methods for open populations are also described, such as the Jolly-Seber model. This model measures several quantities such as the number of individuals caught in each sample, the number of individuals already marked or caught for the first time in each sample, the number of marked individuals released after a sample, among others. These quantities are then used to estimate the number of individuals in the population up to the last sample, and other parameters. However, this model would require a meticulous measurement and control of the number of nodes present in an area in the pass of a nanosatellite. Characteristics, such as the non-homogeneity of the conditions of each node, the need to store the node identities, and the fact that an orbit can take several cycles to repeat itself, make this method difficult (if at all possible) to implement in the scenario studied.

C. Indirect Network Size Estimation Approaches

Other mechanisms approach the estimation problem focusing on the observed channel load or the number of collisions (and other events) occurring during contention periods [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF], [START_REF] Wei | Modeling and analysis of random access channels with bursty arrivals in OFDMA wireless networks[END_REF]- [START_REF] Tello-Oquendo | Efficient random access channel evaluation and load estimation in LTE-A with massive MTC[END_REF]. To address the channel load estimation problem in a Long Term Evolution (LTE) network, the authors in [START_REF] Wei | Modeling and analysis of random access channels with bursty arrivals in OFDMA wireless networks[END_REF]- [START_REF] Tello-Oquendo | Efficient random access channel evaluation and load estimation in LTE-A with massive MTC[END_REF] provide a Probability Distribution Function (PDF) of the result of the access requests to ease the management at the evolved Node B (eNB) when the simultaneous activation of a large number of Machine-Type Communications (MTC)also called User Equipments (UEs) -occurs.

In [START_REF] Wei | Modeling and analysis of random access channels with bursty arrivals in OFDMA wireless networks[END_REF], the authors derive an explicit expression to compute the PDF of the number of successful transmissions. Later, the authors in [START_REF] Duan | D-ACB: Adaptive Congestion Control Algorithm for Bursty M2M Traffic in LTE Networks[END_REF] compute the PDF of the number of successful preamble transmissions using combinatorial analysis. Both contributions show numerical stability issues around 200 devices, the point at which they stop generating valid results. In [START_REF] Tello-Oquendo | Efficient random access channel evaluation and load estimation in LTE-A with massive MTC[END_REF], the authors calculate the joint PDF of the successful and collided access requests within a random access opportunity, increasing the effectiveness up to 300 devices. However, the method is applied to a scenario that differs from the DtS-IoT scenario since it requires terrestrial infrastructure (i.e., the eNB) capable of carrying the computations of PDF and the estimations.

In [START_REF] Sun | Detecting the Number of Active Users in IRSA Access Protocols[END_REF], the authors present an estimation technique to detect the number of active users in irregular repetition slotted ALOHA (IRSA) access protocols based on successive interference cancellation (SIC) at the receiver's side. Each device transmits multiple replicas of a packet in different slots, and the interference contribution of already identified devices is canceled, permitting the identification of new devices. But, this technique requires high processing and memory capacity at the receiver, which prevents its applicability in CubeSats [START_REF] Ferrer | Review and evaluation of mac protocols for satellite IOT systems using nanosatellites[END_REF].

The authors in [START_REF] Shahzad | Fast and Accurate Estimation of RFID Tags[END_REF] present the average run-based tag (ART) estimation scheme for FSA wireless networks in radio frequency identification (RFID). Collisions due to simultaneous tag responses may occur in the same way as in WSNs. ART is based on the average run length of 1 ′ s in the bit string received at the reader (i.e., receiver), considering 1 as the bit value of one or multiple transmissions received in a slot and 0 as the value of no transmission. These authors also provide a set of optimal parameters that save some calculations for wide ranges of node populations, but complex and costly numerical solutions are still needed at the end of each frame during the estimation stage. This makes ART inadequate for a CubeSat implementation.

Another scheme that uses the scalable minimum mean square error (sMMSE) is proposed as an estimation method for the RFID tag population and it actually outperforms the ART method in both the achieved estimation error and the estimation time [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF]. This approach bases its operation on the iterative adaptation of the FSA frame length to the node population to be estimated. Once the frame size is adapted, the estimation is calculated by numerically solving a minimum square error equation. Given the superior performance of sMMSE over ART and the possibility to adapt it to the DtS-IoT scenario studied, this scheme is considered one of the benchmarks of our proposal.

The other benchmark selected for its performance and simplicity is the collision set size estimator proposed by Zanella for the FSA protocol in RFID systems [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF]. A collision set is the group of nodes contending for a time slot in a given contention window of FSA. A collision occurs if two or more nodes in the set transmit concurrently. The receiver can perceive these collisions, keeping track of the number c of collisions that have occurred, as well as the number of successful transmissions s and idle slots i. These three numbers together compose a time slot observation. Then, Zanella applies a maximum likelihood estimator to these values, which returns the number n of transmitting nodes that maximizes the conditional probability of having that observation at the end of the frame. The number of transmissions in each time slot is considered to be an independent Poisson random variable.

Table III presents a summary with a qualitative assessment of the main estimation methods surveyed here and their applicability to DtS-IoT scenarios. In addition, Section III provides a detailed description of the selected benchmark schemes and the adaptations required for their application in the DtS-IoT scenario.

III. NETWORK SIZE ESTIMATION MECHANISM FOR DTS-IOT

In this section, we describe the proposed network estimation mechanism employed in a DtS-IoT network, considering the capabilities and limitations of an IoT gateway deployed on a CubeSat. We first detail the study scenario and assumptions and then describe OCI, a simple yet accurate network estimation mechanism.

A. System Model and Assumptions

The network scenario under study follows a DtS IoT architecture, with a CubeSat nanosatellite acting as the IoT gateway. The CubeSat, whose coverage footprint passes over a cluster of terrestrial IoT nodes, travels around the earth using a LEO orbit. The orbit in which the satellite is deployed has a height and speed of approximately 500 km and 7.5 km/s, respectively [START_REF] Cakaj | The coverage analysis for low earth orbiting satellites at low elevation[END_REF]. The satellite completes a turn around the world in 90 minutes.

For network size estimation purposes, the communication between the CubeSat and IoT nodes follows an FSA logic. The FSA communication windows are illustrated in Fig. 1; each window is divided for downlink and uplink transmissions. In the downlink, a beacon serves the satellite to announce the beginning of the frame, which corresponds to an uplink period divided into w time slots intended for contention. The beacon Random tour and gossip-based aggregation algorithms rely on the interaction between adjacent nodes of the network.

Unsuitable, requires interaction between nodes which is not assured in the DtS-IoT architecture . [START_REF] Bustos-Jimenez | Estimating The Size Of Peer-To-Peer Networks Using Lambert's W Function[END_REF] & [START_REF] Bustos-Jiménez | Estimating the size of natural networks using local information[END_REF] BFS approach based on CR logic for P2P networks. Nodes discover the network size by interacting with neighboring nodes Unsuitable, nodes require P2P communication with other network's nodes. [START_REF] Accettura | The Capture-Recapture approach for population estimation in computer networks[END_REF]- [START_REF] Achtzehn | IEEE 802.11 Wi-Fi access point density estimation with capture-recapture models[END_REF] Survey of CR approaches for population estimation in closed and open computer networks [START_REF] Accettura | The Capture-Recapture approach for population estimation in computer networks[END_REF], and CR approaches for closed wireless networks [START_REF] Peng | Estimation of a Population Size in Large-Scale Wireless Sensor Networks[END_REF], [START_REF] Achtzehn | IEEE 802.11 Wi-Fi access point density estimation with capture-recapture models[END_REF].

Unsuitable for closed and open populations due to the rapid network size changes, the node identity storage requirement, and meticulous control of all nodes present in a region. [START_REF] Wei | Modeling and analysis of random access channels with bursty arrivals in OFDMA wireless networks[END_REF]- [START_REF] Tello-Oquendo | Efficient random access channel evaluation and load estimation in LTE-A with massive MTC[END_REF] Estimation mechanism according to two PDFs of the number of successful and collided preambles during contention periods of LTE networks. Estimation effectiveness for up to 300 devices.

Unsuitable. The calculations to compute the PDF and the estimation are unfeasible considering CubeSat's limitation. Furthermore, the effectiveness is desired for bigger node population. [START_REF] Sun | Detecting the Number of Active Users in IRSA Access Protocols[END_REF] Scheme based on the irregular repeat slotted ALOHA (IRSA) protocol. User devices transmit multiple replicas of their packets in random slots. Their interference contribution is canceled by applying (SIC) techniques at the receiver, allowing the identification of new transmitter devices.

Unsuitable for CubeSats since SIC techniques require high processing and memory capacity at the receiver. Furthermore, the high expected channel congestion impairs the operation of this technique if a frame of adequate size is not used, which risks being excessively large. [START_REF] Shahzad | Fast and Accurate Estimation of RFID Tags[END_REF] Estimation of the tag population size of RFID networks based on the average run size of successful and collided slots in an FSA frame.

Unsuitable for a CubeSat serving multiple regions due to the high processing cost involved in performing complex numerical calculations at the beginning and end of each frame. [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF] Dynamically and iteratively adapting the FSA frame length to the RFID network size according to the channel load. Network size is estimated using a minimum mean square error function, whose result is obtained through numerical solutions.

Not directly applicable. This scheme requires modifications to bear the cost of adapting the frame for each coverage region.

Frame lengths risk being excessively large (e.g. 4096 slots for 1500 nodes with the proposed estimator). [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF] Collision set size estimator based on a Poisson simplification of the maximum likelihood estimator of contender nodes given the observation of the channel load in RFID systems.

Suitable. The channel observation is feasible in a DtS-IoT scenario, and the complexity of the calculations is viable on a CubeSat-type nanosatellite.

is an invite to the IoT nodes against which the only possible response is a transmission attempt, i.e., nodes do not remain idle. Consequently, during the estimation process, every node competes in each available frame for which it receives a beacon. We assume a cluster is stable for the duration of a frame, as depicted in Fig. 1. Time synchronization between nodes and the satellite is assumed to properly operate the FSA logic. Although the synchronization aspect is out of the scope of this work, recent studies show it is feasible to maintain synchronization in DtS-IoT under certain constellation designs [START_REF] Fraire | Sparse Satellite Constellation Design for LoRa-based Direct-to-Satellite Internet of Things[END_REF]. During the estimation process, the information sent by the IoT nodes may be successfully decoded by the CubeSat or not. The failures in the decoding process can be due to a) collisions, i.e. simultaneous reception of signals from different nodes, or to b) insufficient received signal power due to path loss and other communication channel conditions. When the power of the received signal is less than the sensitivity of the receiver, successful decoding is not possible; but note that very weak signals will be undetected by the receiver on the satellite. In such cases, the CubeSat will erroneously assume that time slot to be idle. Moreover, in this work, we call effective detection the valid detection of actual transmissions from the IoT nodes, while erasure happens when the CubeSat erroneously detects an idle time slot. The latter resembles an on-off channel model [START_REF] Tondo | Direct-to-Satellite IoT Slotted Aloha Systems with Multiple Satellites and Unequal Erasure Probabilities[END_REF], [START_REF] Munari | Multiple-relay slotted aloha: Performance analysis and bounds[END_REF].

Although all nodes under the footprint choose a slot upon beacon reception, there can be slots that remain idle if no node chooses that slot for a transmission attempt. Consequently, the satellite perceives three outcomes for a given slot: successful, collided, or idle, and can track the number of failures and successes within a frame.

We consider the activity of the CubeSat divided into two phases: a) an estimation phase and b) an operational phase. During the estimation phase, the CubeSat travels its orbit, announcing the beginning of new frames and decoding the received node transmissions at each time slot. This work aims to estimate the number of contender nodes n k , using the information acquired during the contention of each frame k. The estimated number of contender nodes is noted y k . We consider each frame k to be related to a specific region of the globe, and thus n k corresponds to the number of nodes present in that region. For this, we assume the application of Earthrepeat orbits, further discussed in Section V. Multiple passes over the same region may be needed to estimate n k . In that case, we consider each pass an iteration of the process, and the estimate of frame/region k at iteration m is denoted y (k,m) . M denotes the total number of iterations for the estimation phase. This means that after M passes, the CubeSat advances to the operational phase, using the estimates to throttle the communication with the ground nodes.

The estimation phase may need to be repeated periodically to serve two purposes: 1) to update the network's estimated size at each zone so as to account for dynamics in the network population and 2) to adjust communication parameters according to the changing requirements of the DtS-IoT network. However, determining an optimal periodicity for the updates is out of the scope of this work.

B. Network Size Estimator Design

Following the mathematical development presented by Zanella in [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF], let v = ⟨c, s, u⟩ be one observation of the set V = ⟨C, S, U ⟩ which stores the number of collided, successful, and idle slots in the frame. The number of observations of the set V is numerically equal to the number of passes of the nanosatellite over the same zone considered in the estimation. Then, the conditional probability of observing v given that n nodes transmit can be expressed as

P n (v) = P [V = v|n] (1) = n! c j=0 c j (-1) c-j c-j l=0 c-j l j n-s-l (n-s-l)! w n w s w-s c -1
, where w indicates the length of the FSA frame in terms of slots. The development of this probability mass distribution is further explained in [START_REF] Riordan | An Introduction to Combinatorial Analysis[END_REF].

The classical maximum likelihood estimator is the one that returns the value of n that maximizes the probability shown in [START_REF] Greig | IoT industry projected to have economic impact of $11 trillion by 2025 -TechRepublic[END_REF]. However, such a calculation has several problems, such as the numerical instability of the computation of binomial terms, along with a high computational cost, which makes the maximum likelihood estimator inappropriate for implementation in constrained devices and much less in a nanosatellite. A simplification, suggested in [START_REF] Kodialam | Fast and reliable estimation schemes in RFID systems[END_REF], is to consider the number of transmissions in each slot as an independent Poisson random variable of mean µ = n w . The conditional probability, P [V = v|µ], of observing v given µ can then be expressed as

P µ (v) = µ s e -µw (e µ -1 -µ) c . (2)
Now, to find the value that maximizes that expression, the derivative of (2) in µ is set to zero, which admits only one non-negative solution

µw -s c = µ(e µ -1) e µ -1 -µ , (3)
which can be determined by the bisection search method [START_REF] Eiselt | Linear programming and its applications[END_REF].

The simplification made to get the expression of equation ( 2) allows the computational and energy cost to be adequate for a low-end device. However, the approximation restrains the mechanism to Poisson-type traffic. Consequently, our proposed estimation mechanism shares with [START_REF] Kodialam | Fast and reliable estimation schemes in RFID systems[END_REF] the design objective of using computationally simple calculations but overcoming the Poisson-type traffic restriction over the network traffic. We introduce our proposed estimator in the following section.

C. Optimistic Collision Information-Based Estimator (OCI)

Considering that n k nodes are within the CubeSat footprint at region k, if every node transmits in the frame choosing a slot with probability 1 w , we argue that a simple yet naive way of estimating n k is by the consideration ϕ = s + 2c, where ϕ is the estimation of n k , assuming no more than two nodes are involved in every collision, and that the only way a transmission fails is due to collisions. This assumption is sufficiently accurate as long as n k ≤ w; otherwise, the probability of triple, quadruple, and higher order collisions becomes larger, and thus ϕ becomes an underestimation.

To assess the effectiveness of the naive estimator, we simulated the communication between the CubeSat and nodes using MATLAB. In this simulation, nodes are randomly distributed within the nanosatellite footprint, and a total of K regions are simulated with varying node population sizes. A 100% effective detection is considered; that is, the outcome of a transmission attempt can only be a success or a collision. We employ a step-size of 10 to account for network sizes n k from 10 nodes to 2000 nodes , reaching a total of K = 200 different regions with different values of node population. For these simulations, k = 1, 2, ..., K is a variable designating fictitious regions of increasing size.. The frame length is set to w = 128 as in [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF]. At each step, the contention of n k nodes over the w slots of an FSA frame is simulated, and the results in terms of c, s, and u are stored. The naive estimates ϕ = s + 2c are calculated at each step and stored in the array Φ of length K. The orange curve shown in Fig. 2 reflects the elements of Φ when plotted as a function of the number of nodes contending for a slot. The curve converges to a value equal to twice the number of slots in a frame due to high channel congestion for greater network sizes. However, setting a longer frame length, with w = 512, shows a better result, as depicted by the blue curve in Fig. 2. As long as triple and higher order collisions are not considered, this naive estimation always converges to 2w, limiting its estimation capabilities for large network sizes.

A way to solve this issue, without knowing the order of the collisions, is to create a function that relates each point of the Naive estimation w Fig. 2: Network size estimation with a naive approach (orange and blue curves) and with OCI estimator (green curve) for which the correction coefficients are [7.024e-09,-1.056e-05,0.006,-0.036,41.705] naive estimation curve to the corresponding real value. For this to work correctly, bijectivity must be ensured between the naive estimation and the real value of nodes. Using a polynomial fit solves this matter as long as the frame length w is long enough in terms of slots so that the naive estimation has no repeated values for the simulated domain, thus ensuring injectivity. In a practical scenario, this function would be provided to the CubeSat, in such a way that it only uses the naive estimate method and then calculates the image using the given function, keeping the estimation process as simple as possible. Accordingly, we propose to use a first polynomial fit of order seven to smooth the naive estimation curve, ensuring injectivity. Then, we calculate the coefficients required to relate each point of the smooth curve to the real number of nodes by using a second polynomial fit; the latter maps the elements of the naive estimates ϕ to the value of the real number of nodes. The order of the second polynomial fit has been chosen to be high enough with the purpose that the correction oscillates less around the true values (i.e., the red curve in Fig. 2). In the simulation, the order of the second polynomial fit is set to 4 by inspection; however, the impact of this parameter is further discussed in Section V.

The coefficients of a polynomial of degree q that best fit the estimation y(x), y(x) = p 0 x q + p 1 x q-1 + • • • + p q x 0 , in a least-squares sense, can be calculated from

X • P = Ñ (4)
where X represents the K ×q+1 Vandermonde matrix formed from the elements of the naive estimations array Φ, on a geometric progression. P is a row vector of length q + 1 containing the polynomial coefficients in descending powers to be calculated, and N is a K length vector containing the real value of nodes n k . Thus, the detailed composition of each of these elements is

     ϕ q 1 ϕ q-1 1 • • • 1 ϕ q 2 ϕ q-1 2 • • • 1 . . . . . . . . . . . . ϕ q K ϕ q-1 K • • • 1           p 0 p 1 . . . p q      =      n 1 n 2 . . . n K      (5)
In order to calculate the coefficients in P, the matrix X is decomposed on a orthogonal-triangular manner, also known as QR decomposition, such that X = Q • R, where the R factor is an K-by-(q + 1) upper triangular matrix and Q is an K-by-K orthogonal unitary matrix, such that

Q ⊤ • Q = I.
Then, P is calculated as

P = R \ (Q ⊤ • N). ( 6
)
This process is valid if the polynomial degree q is less than the length of X. Once the coefficients are determined, they can be stored in CubeSat's ROM memory. Due to the processing and energy costs of these calculations, a ground-based station is responsible for the polynomial regression. Architectures with down-link and up-link channels have been developed for nanosatellites such as CubeSats [START_REF] Rahim | Estimation of communication link on Ground Sensor Terminal (GST) system for nanosatellite (UiTMSAT-1) store-and-forward mission[END_REF]- [START_REF] Verspieren | Innovative Approach to Data Gathering in Remote Areas Using Constellations of Store & Forward Communication Cubesats[END_REF], allowing these lowcost devices to download and store data from an earth station. A practical example of the latter is Swarm Technologies, which uses nanosatellite constellations for IoT communication using a store and forward design [START_REF]Swarm -Low cost, global satellite connectivity for IoT[END_REF].

These ground-based stations are also responsible for studying the conditions of the channel and readjusting the coefficients in case of changes. The recalculation of the coefficients is done by integrating the new channel characteristics to the simulation, followed by the polynomial regression. This work represents these channel characteristics using the effective detection percentage.

The CubeSat's detailed process at the end of each frame is depicted in Algorithm 1. Step 1 is the calculation of the naive estimation ϕ using input parameters s and c. After declaring y in step 2, steps 3 to 5 calculate y as the evaluation of ϕ by the polynomial whose coefficients are present in vector P, in decreasing degree order. The calculation corresponds to Horner's method of evaluating a polynomial, which has proven to be optimal in terms of arithmetic operations [START_REF] Pan | METHODS OF COMPUTING VALUES OF POLYNOMI-ALS[END_REF] as it only requires some additions and multiplications equal to the order of the polynomial. Step 6 returns the estimate y (k,m) , which is the equally-weighted average of the estimate of m passes if M > 1, or y if M = 1. The design variable M (i.e., the number of satellite passes devoted to the estimation stage) that governs the iterations m, is finite and should be vastly smaller than the duration of the operation stage of the network. In step 7, the CubeSat overwrites y (k,m-1) with the new estimate y (k,m) .

The result of applying Algorithm 1 to the naive estimates associated with w = 512 and M = 1 are illustrated by the green dots in Fig. 2.

We call our proposed estimation mechanism the Optimistic Collision Information (OCI) based estimator. The initial OCI's estimation curve in Fig. 2 has a noisy behavior, which can be refined by carrying out the estimation in multiple passes (i.e m > 1 in algorithm 1) and progressively averaging the estimations of each pass of the CubeSat over the same zone.

Algorithm 1: OCI algorithm

Input: s, c: number of successful and colliding slots within FSA frame k on iteration m Stored Data: p[q + 1]: array of size q + 1 of the correction coefficients in decreasing degree order M : the number of satellite passes for the estimation y (k,m-1) : estimated number of contenders for frame k at the previous iteration (only if M ≥ 2) Output: y (k,m) : estimated number of contenders for frame k at iteration m 1 ϕ = s + 2c; 2 y = 0; 3 for a = 1 to (q+1) do

D. Additional Benchmark Scheme

In [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF], sMMSE is presented as a two-phase algorithm: a) an adaptation phase where the FSA frame size is dynamically and iteratively increased to be adapted to the number of devices and b) the estimation phase with a minimum square error function. The adaptation phase consists of increasing the frame length if the percentage of slots where transmission was received is larger than a predefined threshold. The process is repeated until the response rate RR is lower than the threshold, meaning that the frame length is adapted to the node population size to be estimated. An optimal set of parameters is provided in [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF] so that the frame is doubled when 40% of its slots are occupied. The starting frame size of the algorithm is 128 slots in [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF], resulting, as an example, in a length of 4096 slots for a 1500 node population, representing a considerable resource efficiency disadvantage in comparison to other schemes with a fixed frame size. Another downside is the number of iterations needed to adapt the frame, which implies an additional adaptation cost. To adapt sMMSE to the requirements of our scenario, the frame is started at an initial length of 512 slots, which reduces the number of iterations. The explanation for this choice is found in Section III.D.

The estimation phase starts right at the end of the adaptation phase. This phase consists of calculating the estimated number of nodes n that minimize the distance between the observed tag responses v f and the expected value V f , according to

H sM M SE (v f ) = arg min n∈N P n (v f ), (7)
which, using a binomial distribution of the probability of node occupancy and Chebyshev's inequality [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF], is expressed as

P n (v f ) = ||v f -V f || 2 (8) = (RR -1 + ψ n ) 2 +(1 -RR -ψ n ) 2 ,
where RR represents the response ratio, ψ = 1 -1 kL , k and L are parameters such that kL represents the FSA window size, and n is the number of devices to be estimated.

IV. PERFORMANCE EVALUATION

In this section, we provide an analysis of the performance of OCI. We selected the Zanella (henceforth Zan) collision set size [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF] and the sMMSE estimators [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF], discussed in Sections II and III, for comparison purposes. Furthermore, we evaluate the impact of using network size estimators to provide feedback to the MAC protocol in the scenario of study.

A. Network Size estimator Performance

To measure the performance of the estimators, we calculate the Root-Mean-Square Error (RMSE), commonly used to rate how a regression model fits a dataset. For this, we compare the elements of the array of the real number of nodes N present in the satellite footprint to the ones of the array of the estimated size with the equation

ε RM SE = K h=1 (y (h,M ) -n h ) 2 ) K (9)
RMSE is heavily penalized when an estimate deviates far from the actual observed value. For this reason, it is used in cases where significant errors are undesirable. Zan and sMMSE estimations are calculated for benchmarking purposes to compare to the performance of OCI. In previous work, we have adapted and evaluated the Zan estimator to DtS-IoT scenarios [START_REF] Ilabaca | Network size estimation in direct-to-satellite iot[END_REF]. All simulations are developed using MATLAB R2020a, with a fresh implementation of OCI and sMMSE, and using Zan's code as provided by its author [START_REF]dgt Dei Gruppo Telecomunicazioni[END_REF]. To ensure fair comparisons, sMMSE's frame size is set at 512 as its starting point.

To evaluate the performance, the first set of simulations considers n k nodes distributed randomly within the nanosatellite footprint, ranging from 10 to 2000 with steps of 10, for a total of K = 200 different node populations. At each size-step, n k nodes contend for a random slot in the FSA frame, and the estimations of Zan and OCI are calculated accordingly, using the record of successes and collisions, s and c. Effective detection of 100% is considered. Furthermore, for each k = 1, 2, ..., K, we consider multiple satellite passes M , ranging from 1 to 200. For each pair (k, M ), the estimation results of OCI and Zan are denoted as y (k,M ) and z (k,M ) , respectively. Note that the estimate y (k,M ) is equal to the average of all estimates y (k,m) for m = 1, 2, ..., M . The matrices Ȳ and Z, expressed in [START_REF] Ferrer | Review and evaluation of mac protocols for satellite IOT systems using nanosatellites[END_REF], are constructed using the aforementioned OCI and Zan estimates, respectively.

y (k,m) = m h=1 y (k,h) m z (k,m) = m h=1 z (k,h) m Ȳ =         y (1,1) • • • y (k,1) • • • y (K,1)
. . . . . . . . . . . .

y (1,m) • • • y (k,m) • • • y (K,m) . . . . . . . . . . . . y (1,M ) • • • y (k,M ) • • • y (K,M )         Z =         z (1,1) • • • z (k,1) • • • z (K,1) . . . . . . . . . . . . z (1,m) • • • z (k,m) • • • z (K,m) . . . . . . . . . . . . z (1,M ) • • • z (k,M ) • • • z (K,M )         ( 10
)
Each row of Ȳ and Z compose a different estimation phase of m passes. We calculate the RMSE between rows and the real values of nodes that were set up in the simulations, resulting in the 200 points plotted for each mechanism in Fig. 3. The adaptation phase of sMMSE is not considered for the count of satellite passes. This means that the first satellite pass is considered to be carried out when the adaptation has already been made in all regions. Because of this, the different estimates of each pass result from the numerical solution of ( 7) with a fixed frame size but a different response rate (RR) at the end of the communication window.The simulation is then repeated, considering three different effective detection percentages: 100%, 95%, and 75%.

Results in Fig. 3 show that OCI outperforms Zan and sMMSE in scenarios with more realistic effective detection rates, lower than 100%. Zan's RMSE increases when the effective detection percentage decreases, since the Poisson approximation of the traffic is no longer valid. The Zan mechanism underestimates the number of nodes, given that a lower effective detection percentage means that a portion of nodes does not manifest their presence in the frame. Something similar happens with sMMSE as the response rate (RR) is decreasing due to reduced detections, and thus sMMSE underestimates the number of nodes. On the contrary, OCI's correction coefficients implicitly have this information, which explains why the resulting RMSE remains at a low value with lower effective detection rates. Additionally, OCI's RMSE decreases as the effective detection ratio is reduced. When the collision probability decreases, the probability of multiple collisions (triple and further) also decreases, rendering OCI's approximation more accurate. On average, Zan's RMSE is 4 times higher than OCI's considering a 95% effective detection scenario, and 38.7 times higher on a 75% effective detection scenario. Zan's RMSE is close to sMMSE's when imperfect detection ratios are considered.

When a higher number of satellite passes is considered for the estimation phase, the RMSE decreases progressively, but at a decreasing rate for the three estimators. Nevertheless, none of the RMSEs reached a point of convergence, as shown in -+-OCI, 75% EDR 10 °-'-r----------�-�----------�-------< 10 0 10' 10 2 2xl0 2 Satellite passes Fig. 3: RMSE (log scale) obtained from the simulations of the three estimators under different effective detection rates (EDR) and up to 200 satellite passes in the estimation phase. Initial frame length w = 512 slots for Zan, sMMSE (assuming the adaptation phase is finished), and OCI. Fig. 3 where the RMSE lowers asymptotically. Overall, the results on RMSE are good for all the estimators since, when confronted with a large number of nodes, the error between the estimate and the real number is slight.

These results show the adaptability of the OCI mechanism over different scenarios, maintaining an RMSE lower than Zan and sMMSE when there is a nonzero percentage of signal loss. Space-ground channels are more prone to packet errors and low signal-to-interference-plus-noise ratio (SINR) values than ground IoT links. With this in mind, OCI's effectiveness under imperfect detection ratios makes it a better candidate for DtS-IoT.

An important aspect to compare is the efficiency of channel resources. In an environment such as the DtS-IoT, it is desirable to minimize the number of wasted channel resources, as they may negatively impact the network's performance. Fig. 4 shows a comparison between the frame length of each scheme. Note that the low RMSE in sMMSE's estimation goes along with low channel resource efficiency, as it has to increase the frame length to keep an occupancy rate lower than 40%. Furthermore, note that we have omitted the adaptation phase of sMSSE in Fig. 3; however, when the two phases are considered, the FSA frame adaptation phase of sMMSE requires four satellite passes to reach a 512-frame length as the starting point for the estimation. In such a case, the green curve of sMMSE in the figure would shift to the right, rendering a performance similar to Zan's. In addition, considering the two phases, sMMSE could only provide an estimate at the 5 th pass, whereas OCI and Zan can start the estimation process from the first pass.

B. FSA with Network Size Feedback

In [START_REF] Zhao | Slotted aloha game for medium access control in satellite networks[END_REF], the Slotted Aloha Game is proposed to control access to the medium in a satellite network. The proposal focuses on determining the transmission probabilities of a group of Fig. 4: Frame length comparison between sMMSE, Zan, and OCI versus the number of estimated nodes. independent terminals (nodes) that share the communication channel and transmit via the FSA protocol. The probability of transmission of each node depends on the number of active nodes, in such a way that if this number is less than or equal to the number of slots in the frame, they are allowed to transmit with a probability equal to 1. However, if the number of active nodes exceeds the number of slots, the probability of transmitting decreases in such a way that, on average, only a number of nodes equal to the size of the frame transmit. The evaluation of the Slotted Aloha Game shows to allow the communication of a large number of nodes, maintaining an approximately constant throughput, with a low access delay and a low power consumption [START_REF] Zhao | Slotted aloha game for medium access control in satellite networks[END_REF]. Similarly to [START_REF] Vogelgesang | Uplink transmission probability functions for LoRa-based direct-to-satellite IoT: A case study[END_REF], the authors assume that the number of active nodes is perfectly known.

In this section, we use the network size estimates to set the threshold employed by an FSA protocol following the aforementioned Slotted Aloha Game. Two metrics are calculated to evaluate the performance: throughput and normalized energy efficiency. Both pure FSA and FSA using the Slotted Aloha Game approach fed with network size feedback are implemented in the simulation for comparison purposes. First, a scenario with 100% effective detection is simulated, where Fig. 5a shows the throughput calculated as the number of successful transmissions in a frame divided by the time length of the frame (i.e. w = 512), considering one pass of the satellite for the estimation (M = 1). Then, Fig. 5b shows the normalized energy efficiency, which is calculated considering that when a sensor node transmits, it consumes one arbitrary unit of energy as in [START_REF] Zhao | Slotted aloha game for medium access control in satellite networks[END_REF]; then, the energy consumed by successful transmissions is divided by the total energy consumed during the entire process for transmissions. The results represent the averaged results of 30 repetitions for each estimator. The results show the performance for twice the number of nodes that OCI was designed for. This is to appraise the scalability of the designed estimation mechanism, and the viability of letting the node population grow without re-adjusting the polynomial coefficients and frame length. and5b, it can be seen that Pure FSA's throughput decreases once the number of nodes contending exceeds the number of slots in the frame and its energy efficiency decreases with the number of nodes. When either of the three estimators' feedback is applied, these two metrics have the same behavior until the network's size exceeds the frame length. From that point, throughput and energy efficiency oscillate close to 0.37. However, with a number of nodes higher than 2000, OCI starts to show slightly worse performance than sMMSE and Zan, while the latter shows discontinuities when the network's size is close to 4000 nodes. These results allow us to affirm that the feedback from the three estimators benefit the MAC protocol performance even if the number of nodes for which its parameters were adjusted is exceeded. Zan and sMMSE show the best result for this scenario.

On the one hand, the discontinuities that Zan presents (seen by the red dots at abscissa 0), which translate into 0 throughput, and the slightly worse performance of OCI are explained by the number of nodes for which the estimators' parameters were chosen. The algorithm designed by Zanella in [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF] returns an estimation equal to infinity once the number of nodes increases past the point where all the slots in a frame collide. On the other hand, OCI's polynomial coefficients are adjusted for a setting with a maximum of 2000 nodes, for which a frame length of w = 512 slots is suitable. Once the number of nodes increases considerably beyond the original maximum, the estimator progressively losses estimation capability until it converges to an estimate of 2415 nodes for a population size of 4000 nodes. Nevertheless, the throughput is maintained close to the theoretic bound by allowing 21.2% of the present nodes to transmit [START_REF] Zhao | Slotted aloha game for medium access control in satellite networks[END_REF], which is the percentage that w = 512 represents to a 4000 node population. To improve the performance for larger network sizes, the frame length should be increased, and the coefficients for OCI's correction recalculated for the new network size. However, depending on the characteristics of the application, it might be unsuitable to adopt a number of slots past a certain length. sMMSE's estimation needs first an adaptation period which explains why its behavior is maintained close to the theoretic bound even for large node population sizes.

We also carried out a simulation with up to 200 satellite passes; the results showed that a single satellite pass devoted to the estimation process is enough to improve the performance, and thus a plot of these results is not shown in this manuscript. Making use of more passes narrows the resulting points in the graph but does not bring any substantial improvements.

To account for more realistic scenarios, we simulated with a 75% effective detection at the satellite. Fig. 6a shows a slightly worse performance, in terms of throughput, of FSA with Zan's and sMMSE's feedback compared with using OCI's feedback. Due to the 25% of ineffective detection, the frame is less congested, which translates to higher probabilities of successful messages. This latter dampens the effects of the underestimation of Zan and sMMSE and explains the difference in performance. When comparing energy efficiency, Fig. 6b shows a bigger difference in favor of OCI's feedback since the better estimate allows better control on node transmissions, mitigating energy waste. After the node count exceeds 2000, the number for which OCI was designed, the performance decreases towards matching Zan's and sMMSE's feedback performance at 4000 nodes. Overall, the results of Figs. 5a and 6a show that the peak throughput is achieved at a network size equal to the frame length. Nevertheless, when using OCI, the throughput is maintained at a high level for large network sizes of up to 4000 nodes. Hence, the supported network size is eight times the number of nodes where the throughput started to decrease for the Pure FSA scheme.

V. PRACTICAL CONSIDERATIONS In this section, we discuss relevant practical aspects involved in the estimators' implementation, such as (i) the applicability

A. Arrival Process

One advantage of the OCI estimator is its independence from the nature of the data arrival process at the satellite. For the Zan estimator to be applicable in restricted IoT gateways, it is necessary to use the simplification in which the number of transmissions in each slot of the frame is considered to be an independent Poisson random variable. Such a restriction is not present in OCI. However, the OCI estimator is limited to be predefined for a precise frame length used during the estimation process. To mitigate this limitation, when the number of slots in the contention window is expected to change, a lookup table can be stored in CubeSat's ROM, containing the coefficients needed for the algorithm adjusted to different window sizes.

B. Computational Cost

Computational cost is an area where differences in estimators can be observed. In the case of the Zan estimator [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF], it increases with the frame length w, according to log 2 w. Instead, OCI's computational cost is independent of the frame length but grows with the order of the polynomial function used for the correction process depending on the polynomial evaluation method. For instance, Horner's method [START_REF] Pan | METHODS OF COMPUTING VALUES OF POLYNOMI-ALS[END_REF] has a complexity O(q), where q represents the order. sMMSE's complexity depends on the numerical method used to solve its equation. To this end, and to validate the complexity analysis, OCI and Zan estimators were empirically tested on a Raspberry Pi Zero W board, which is currently used as the on-board computer (OBC) of CubeSat implementations such as in [START_REF]OBC: On Board Computer -ECE3SAT[END_REF]. The hardware characteristics of the emulated OBC are specified in Table V. The test is based on calculating the average time needed to compute 10 5 repetitions of estimations using both methods and considering a contention window w = 512 slots. The experiment concluded with both estimators being suitable for a CubeSat in terms of computational cost: OCI shows a time requirement of 8.13 • 10 -5 s while Zan shows a requirement of 3.35 • 10 -4 s. The reported times demonstrate the advantage of OCI in terms of computations required over Zan, as OCI needed only 24.3% of the time demanded by Zan's method while achieving very similar or better performance. However, the sMMSE simulations could not be run on the same hardware platform, so to fairly compare the three estimators in terms of computation time, the same experiments were run on a computer. Now, the average time used to compute 10 5 repetitions of estimations was calculated using Python 3.10.7. The results, summarized in Table IV, show 1.22 • 10 -3 s on average for sMMSE, 5.48 • 10 -6 s for Zan, and 7.59 • 10 -7 s for OCI. For this test, only the estimation phase was considered for sMMSE, using k = 128 as done in [START_REF] Arjona | Scalable RFID Tag Estimator With Enhanced Accuracy and Low Estimation Time[END_REF] and different values for RR and L calculated by independent simulations of the first phase. These results confirm OCI as the fastest scheme since OCI needed only 13.9% of the time required by Zan's method and 0.06% of the time required by the sMMSE method. To further compare Zan and OCI in terms of time requirements, we carried out a simulation varying the frame length from 128 to 2048 slots with steps of 128. The results are shown in Fig. 7 where time is normalized to the maximum value reached by the Zan estimator when a frame of 2048 slots is used. As can be seen, Zan's computation cost increases along with the frame length, whereas OCI's time requirement remains at the same values except for small random variations. Furthermore, the Zan estimator is proven to be more time efficient than the best-performing estimators for RFID tags up to the year of its publication [START_REF] Zanella | Estimating collision set size in Framed Slotted Aloha wireless networks and RFID systems[END_REF], which means that OCI also outperforms methods like [START_REF] Vogt | Efficient Object Identification with Passive RFID Tags[END_REF] in terms of computational cost. In [START_REF] Bui | A Comprehensive Distributed Queue-Based Random Access Framework for mMTC in LTE/LTE-A Networks with Mixed-Type Traffic[END_REF], a look-up table is proposed to store all the possible estimations for different combinations of FSA frame results (only two among s, c and u due to the constraint of s + c + u = w) of Vogt's method [START_REF] Vogt | Efficient Object Identification with Passive RFID Tags[END_REF] to improve time delays. For OCI, a look-up table could also be implemented. This would imply negligible delays but there is a need for storing the polynomial evaluations of all possible frames, which results with a fixed w.

We also compared the memory needed for both estimators. Since the satellite travels at approximately 7.5 km/s in its orbit, it goes around the world in 90 minutes, i.e., 5400 seconds. It is common in the literature that the length of a slot in a communication protocol such as FSA lasts an amount of time equal to the time of packet transmission. In the existing literature, applications with packet lengths of around 1000 bits can be found [START_REF] Su | Randomly Pre-coded Packets based Random Access Scheme for IoT-Oriented Satellite Networks[END_REF]; the packet transmission times vary depending on the IoT application. If the transmission of a packet takes for example 5 ms [START_REF] Li | Age-Optimal HARQ Design for Freshness-Critical Satellite-IoT Systems[END_REF], then a frame takes 512 • 5 = 2560 ms, that is 2.56 seconds. Then, with that frame length, there could be 2109 contention windows, i.e. positions, in a turn around the world. In the first satellite pass, OCI needs to store two integers of 3 digits maximum (the number of successes and collisions among the frame slots) and the For OCI to work properly, the set of coefficients required for the correction is to be stored in memory. These coefficients are numbers with up to 15 digits that may require up to 64 bits each, which then implies that 64 • C bits are required, with C being the polynomial order. Fig. 8 shows the impact that the order of the polynomial correction has on the RMSE of OCI under the same experiment carried out in Section III. In Fig. 8 each label represents the error percentage of each point with respect to the previous order. It can be seen that a coefficient of order 4 is sufficient to have a low RMSE compared to lower orders. If a 4 order is chosen, then 256 bits would be needed.

C. Repeating Orbits

As mentioned in Section II, one important issue in estimating the size of the network beneath a passing by satellite is that the ground track does not necessarily repeat. In the literature, orbits that repeat the exact same ground track are known as Earth-repeat orbits [START_REF] Montenbruck | Satellite orbits: models, methods, and applications[END_REF], [START_REF] Lara | Searching for repeating ground track orbits: a systematic approach[END_REF]. Such Earth-repeat orbits exploit the nodal precession effect (a phenomenon provoked by the non-spherical nature of the Earth) to shift the orbital plane such that the ground track matches a previous one [START_REF] Vtipil | Determining an earth observation repeat ground track orbit for an optimization methodology[END_REF]. In other words, the nodal precession is used to balance out the offset of the revolution of Earth. Thus, the same ground track can be covered after a controlled number of orbital revolutions in a given time interval [START_REF] Montenbruck | Satellite orbits: models, methods, and applications[END_REF]. Although the efficient determination of Earth-repeat orbits demands specific computation and orbit maintenance techniques [START_REF] He | High-precision repeatgroundtrack orbit design and maintenance for earth observation missions[END_REF], it unlocks the possibility of learning from reiterative passes over the same region of the network. The longitudinal rotation after a certain period of time of a planet is

∆L 1 = -2π T T E , (11)
while the quantization of the effect of the nodal precession can be expressed as

∆L 2 = - 3πJ 2 R 2 e cos(i) a 2 (1 -e 2 ) 2 , ( 12
)
where T E a sidereal day, J 2 is the Earth's second dynamic form factor, R e is the Earth's radius, i is the orbital inclination, a is the orbit's semi-major axis, e is the orbital eccentricity, and T represents the elapsed time, which can be expressed as a function of the standard gravitational parameter for the Earth µ, like T = 2π a 3 µ . Therefore, a repeat orbit sees these two effects canceling each other after a set of orbital revolutions j and sidereal days k, such that j|∆L 1 + ∆L 2 | = k2π. Then, the required condition for the orbit to repeat can be written as j -2π

2π a 3 µ T E - 3πJ 2 R 2 e cos(i) a 2 (1 -e 2 ) 2 = k2π. ( 13
)
Replacing the variables with the parameters of the Earth and using similar parameters to those used in a known Chilean CubeSat named Suchai [START_REF] N2yo | SUCHAI Satellite details 2017-036Z NORAD 42788[END_REF], [START_REF] Chris Peat | SUCHAI -Orbit[END_REF], it is possible to identify different repeat orbits, as shown in Fig. 9. In this figure, dots of different colors represent a set of orbital inclination and orbital height parameters that holds the equation [START_REF] Zhao | Slotted aloha game for medium access control in satellite networks[END_REF]. Darker colors represent shorter repetition periods, and bigger markers indicate higher accuracy in the resolution of the equation [START_REF] Zhao | Slotted aloha game for medium access control in satellite networks[END_REF]. It should be noted that the parameters of the Suchai nanosatellite are misaligned since its orbit is not repetitive. Up to this point, the estimation process has been considered to be held by a single CubeSat. However, several studies on store and forward technology and its application to CubeSats have been done over the years [START_REF] Rahim | Estimation of communication link on Ground Sensor Terminal (GST) system for nanosatellite (UiTMSAT-1) store-and-forward mission[END_REF], [START_REF] Abbasi | Store & forward communication payload design for leo satellite systems[END_REF], [START_REF] Addaim | Design of store and forward data collection low-cost nanosatellite[END_REF]. With store and forward, a CubeSat can store its estimate on a ground terminal which can later be sent to another CubeSat of the same constellation. In this way, the estimation process done over multiple passes of a single CubeSat can be carried out by a single pass of multiple CubeSats.

VI. CONCLUSIONS

The massive and worldwide deployments of IoT networks face connectivity challenges related to the scalability, coverage, resilience, and ubiquity of these implementations. The Direct-to-Satellite IoT paradigm allows cost-effective solutions such as CubeSat nanosatellites in Low-Earth Orbit that serve as orbiting gateways to remote IoT nodes. However, achieving efficient use of temporal and spectral resources to overcome the exponential connectivity demands remains one of the main challenges of DtS-IoT. This efficient allocation of transmission resources is strongly favored when MAC protocols are aware of the number of contending nodes, but to the best of our knowledge, network size estimation mechanisms operating in IoT gateways deployed in continuously-moving LEO CubeSats have not yet been addressed in the literature. To tackle this issue, in this work we proposed OCI, a novel and computationally efficient network size estimator for DtS-IoT.

The results of the experiments carried out showed that OCI has low estimation error and computational costs, and that the estimate can be used to assist a MAC protocol, achieving a stable throughput serving up to 4000 nodes without energy efficiency penalties. Channel impairments, such as erasures, can be integrated into OCI's correction coefficients to adjust the estimates to the specific scenario for which it will be implemented. This provides OCI with adaptability in low detection ratio scenarios absent in the state-of-the-art estimation mechanisms used for benchmarking, further positioning OCI as a better candidate for DtS-IoT.

Regarding future work, we intend to design new approaches that favor an equitable distribution of channel resources in terms of medium access. Given the nature of a DtS-IoT network with Frame Slotted Aloha, the work in progress aims to provide terrestrial nodes with the capability to decide in which frame to transmit, given the estimated potential competition in each zone obtained with OCI estimations.

Fig. 1 :

 1 Fig. 1: Characteristics of the DtS communication.

4 y

 4 = y • ϕ + p[a]; 5 end 6 Return y (k,m) = y (k,m-1) • m-1 m + y • 1 m ; 7 Replace in memory: y (k,m-1) by y (k,m)

  -sMMSE, 100% EDR •••••• Zan, 100% EDR --+-OCI, 100% EDR _...,_ sMMSE, 95% EDR •• •••• Zan, 95% EDR _...,_ OCI, 95% EDR -+-sMMSE, 75%EDR•••••.,,, ... ..... Zan, 75%EDR

Fig. 5 :

 5 Fig. 5: Comparison of FSA and FSA with size estimation feedback, considering the average of 10 passes and a 512 slots window size for Zan and OCI Analyzing the 100% effective detection scenario, in Figs. 5aand 5b, it can be seen that Pure FSA's throughput decreases once the number of nodes contending exceeds the number of slots in the frame and its energy efficiency decreases with the number of nodes. When either of the three estimators' feedback is applied, these two metrics have the same behavior until the network's size exceeds the frame length. From that point, throughput and energy efficiency oscillate close to 0.37. However, with a number of nodes higher than 2000, OCI starts to show slightly worse performance than sMMSE and Zan, while the latter shows discontinuities when the network's size is close to 4000 nodes. These results allow us to affirm that the feedback from the three estimators benefit the MAC protocol performance even if the number of nodes for which

Fig. 6 :

 6 Fig. 6: Comparison of FSA and FSA with size estimation feedback, considering the average of 10 passes and windows of 512 slots for Zan and OCI

Fig. 7 :

 7 Fig. 7: Normalized computation time required for both estimators varying the frame length.

Fig. 8 :

 8 Fig.8: RMSE of OCI as a function of the order of the polynomial fit used for the correction.

Fig. 9 :

 9 Fig.9: Parameters of orbits that hold equation (13) for circular orbits with varying height and inclination. The current parameters for the Chilean Suchai satellite are highlighted.

TABLE I :

 I List of Acronyms

	Acronym	Definition
	5G	Fifth Generation mobile communication
	BFS	Breadth-First Search
	CR	Capture-Recapture
	DtS-IoT	Direct-to-Satellite IoT
	eNB	evolved Node B
	FSA	Frame Slotted Aloha
	GEO	Geosynchronous Equatorial Orbit
	GPS	Global Positioning System
	IoT	Internet Of Things
	LEO	Low Earth Orbit
	LR-FHSS Long Range-Frequency Hopping Spread Spectrum
	LTE	Long Term Evolution
	MAC	Medium Acces control
	MEO	Medium Earth Orbit
	MTC	Machine-Type Communication
	NB-IoT	NarrowBand-IoT
	OCI	Optimistic Collision Information based estimator
	P2P	Peer-to-Peer
	PDF	Probability Distribution Function
	RFID	Radio Frecuency Identification
	RMSE	Root-Mean-Square Error
	ROM	Read Only Memory
	S&F	Store and Forward
	SINR	signal-to-interference-plus-noise ratio
	UE	User Equipment
	WSN	Wireless Sensor Network

TABLE II :

 II List of Variables

	q	Degree of the polynomial
	s	Number of successful slots in a FSA frame
	u	Number of idle slots in a FSA frame
	w	Frame length as number of slots
	y	Estimate correction of the contender nodes by polynomial
		evaluation
	y (k,m)

Variable Definition Φ Array of the naive estimates ϕ ϕ naive estimate ϕ = s + 2c K Maximum number of FSA frames considered M Number of iterations considered for the estimation phase N Number of nodes deployed in the cluster ⃗ N Array of the node populations n k , for k = 1, 2, ..., K ⃗ P Array of size q formed by the polynomial coefficients T Arbitrary time X Vandermonde matrix of dimensions j × q + 1 formed by the elements of Φ c Number of collided slots in a FSA frame k Frame related to a given CubeSat coverage zone 0 < k ≤ K m Iteration number, 0 < m ≤ M n k Number of nodes present at a given CubeSat coverage zone k, 0 < n ≤ N p Array of polynomial coefficients in decreasing degree order

TABLE III :

 III State of the Art

	Reference	Strategy	Suitability for Dts-IoT
	[15]

TABLE IV :

 IV Computational Costs

	Scheme	Avg. Compute Time	Ratio
	sMMSE	1.22 • 10 -3 s	100 %
	Zan	5.48 • 10 -6 s	0.45%
	OCI	7.59 • 10 -7 s	0.06%

TABLE V :

 V Main characteristic of the emulated CubeSat onboard computer

	Component Characteristic
	Processor	Broadcom BCM2835 ARM1176
	Architecture	ARMv6 (32bit)
	Speed	1.0 GHz
	RAM	512 MB
	OS	Raspbian 10 (buster) Linux 5.10.17

This work has been supported by the Project ANID Fondecyt Regular No. 1201893, the ANID Basal Project FB0008, the ANID FONDECYT Iniciación No. 11200659, FONDEQUIP-EQM180180, Project STARS STICAMSUD 21-STIC-12 Folio STIC2020003, and by CNPq Brazil (402378/2021-0, 305021/2021-4). This work started out while P. Ilabaca and S. Céspedes were with the Department of Electrical Engineering and with the NIC Chile Research Labs, Universidad de Chile, Santiago, Chile.