QoS-Driven Satellite Constellation Design for LEO Satellite Internet of Things

TL;DR

This paper introduces a cost-effective LEO satellite IoT constellation design algorithm that ensures QoS requirements, featuring a novel fitness function and efficient operators for faster convergence and lower costs.

Contribution

It presents a new algorithm with a unique fitness function for designing LEO satellite IoT constellations, optimizing cost and QoS, with proven theoretical convergence.

Findings

The proposed algorithm converges faster than baseline methods.

It achieves lower constellation construction costs under the same QoS.

Simulation results validate the algorithm's effectiveness in practical scenarios.

Abstract

Low Earth orbit (LEO) satellite Internet of Things (IoT) has been identified as one of the important components of the sixth-generation (6G) non-terrestrial networks (NTN) to provide ubiquitous connectivity. Due to the low orbit altitude and high mobility, a massive number of satellites are required to form a global continuous coverage constellation, leading to a high construction cost. To this end, this paper proposes a LEO satellite IoT constellation design algorithm with the goal of minimizing the total cost while satisfying quality of service (QoS) requirements in terms of coverage ratio and communication quality. Specifically, with a novel fitness function and efficient algorithm's operators, the proposed algorithm converges more quickly and achieves lower constellation construction cost compared to baseline algorithms under the same QoS requirements. Theoretical analysis proves the global and fast convergence of the proposed algorithm due to a novel fitness function. Finally, extensive simulation results confirm the effectiveness of the proposed algorithm in LEO satellite IoT constellation design.

Tables

Table 1. TABLE I: Key Limitations of Existing Techniques and Proposed Solutions

References	Optimization Algorithms	Devices’ Locations	Design Metrics
			Coverage	Communication quality	Cost
[20]	Non-dominated sorting genetic algorithm II	Uniform distribution	✓	✓
[21]	Genetic algorithm	Not given		✓
[22]	Progressive construction scheme	Uniform distribution		✓
[23]	Multilayer tabu search algorithm	Not given	✓	✓
[24]	Tabu search Genetic algorithm with opposite-based learning strategy	Uniform distribution	✓		✓
[25]	Non-negative least squares	Not given	✓
[26]	Two-step optimization method	Not given	✓	✓
[27]	Multi-objective optimization framework	Not given	✓	✓
[28]	Two-stage algorithm based on reverse design	Not given	✓	✓
This work	Improved genetic algorithm	Poisson point process	✓	✓	✓

Table 2. TABLE II: Computational Complexity Comparison

Proposed algorithm	$N_{\mathrm{it}}{N}_{\mathrm{totp}}$
Tabu Search [23]	$N_{\mathrm{it}}\left(N_{\mathrm{totp}}^2+N_{\mathrm{totp}}+l_{ts}\right)$
Sine Cosine algorithm [48]	$N_{\mathrm{it}}{N}_{\mathrm{totp}}$
Particle Swarm optimization [50]	$N_{\mathrm{it}}{N}_{\mathrm{totp}}$
Grey Wolf Optimization [49]	$N_{\mathrm{it}}{N}_{\mathrm{totp}}$
Classical GA [47]	$N_{\mathrm{it}}{N}_{\mathrm{totp}}$

Table 3. TABLE III: Main System Parameters For LEO Satellite IoT Constellation

Parameter	Value
Satellite orbit	LEO
Frequency band	L/S/C
Earth’s radius	6378.14 km
Target area	[${60}^{\circ }$S, ${60}^{\circ }$N]
Carrier frequency $f_c$	5 GHz
Carrier bandwidth $B_w$	250 MHz
Satellite antenna gain $G_{\mathrm{sat}}$	17 dBi
IoT device antenna gain $G_{\mathrm{dev}}$	3 dBi
Average rain fading $r_0$	-2.6 dB
Rician factor ${\lambda }_{m,k}$	10 dB
Noise variance ${\sigma }_0^2$	$-106$ dBm
Number of satellite antennas $M_a$	$16$
Activity probability $\epsilon$	$0.005$
Density of IoT devices $\lambda$	$8\times {10}^{-5}$
Average transmit power $\xi$	$3$ dBw
Sequence length $L$	$100$
Satellite coverage angle $\beta$	${45}^{\circ }$
Coverage ratio requirement ${\eta }_{th}$	90%
Backhaul capacity requirement $C_{th}$	$80$ Mbps
Satellite weight $W_{\mathrm{sat}}$	$227$ kg
Phase Factor $F$	1
Start time	1 Jan 2025 00:00:00.0
End time	1 Jan 2025 23:59:00.0
Time step	60 s

Table 4. TABLE IV: Optimized Parameters Range of LEO Satellite IoT Constellation

Parameter	Value
Minimum altitude $h_{\min }$	500 km
Maximum altitude $h_{\max }$	1800 km
Minimum number of planes $P_{\min }$	4
Maximum number of planes $P_{\max }$	20
Minimum number of satellites per plane $N_{\min }$	4
Maximum number of satellites per plane $N_{\max }$	20
Minimum inclination angle $i_{\min }^{\mathrm{sat}}$	${20}^{\circ }$
Maximum inclination angle $i_{\max }^{\mathrm{sat}}$	${60}^{\circ }$

Table 5. TABLE V: Simulation Parameters for the Proposed Algorithm with Classical GA

Parameters	Proposed Algorithm	Classical GA
Population size $N_{\mathrm{totp}}$	30	30
Number of iterations $N_{\mathrm{it}}$	50	50
Mutation operator threshold	0.3	0.3
Penalty factor ${\rho }_1$ and ${\rho }_2$	$\backslash$	1000, 1000
Hyperparameters ${\alpha }_1$ and ${\alpha }_2$	2, 1	$\backslash$

Equations

$$\left(\begin{array}{rl}{x}_{\mathrm{sat}}& \\ y_{\mathrm{sat}}& \\ z_{\mathrm{sat}}& \end{array}\right)=(h+R_e)\left(\begin{array}{rl}\cos (\omega +\upsilon )\cos \Omega & -\sin (\omega +\upsilon )\sin \Omega \cos i\\ \cos (\omega +\upsilon )\sin \Omega & +\sin (\omega +\upsilon )\cos \Omega \cos i\\ \sin (& \omega +\upsilon )\sin i\end{array}\right)$$

$$\nabla u=\frac{360}{(N\ast P)}\ast F.$$

$$\phi_{\max}=\arccos\big{(}\frac{R_{e}}{R_{e}+h}\cos\theta_{\min}\big{)}-\theta_{\min}.$$

$$\beta=2\cdot\arcsin\big{(}\frac{R_{e}}{R_{e}+h}\cos\theta_{\min}\big{)}.$$

$$S_{\mathrm{sat}}=2\pi R_e^2(1-\cos {\varphi }_{\max }).$$

$$o_{\hat{g}}(t)=\{\begin{array}{rl}& 1\mathrm{within}\mathrm{coverage}\\ & 0\mathrm{beyond}\mathrm{coverage}\end{array}.$$

$${\eta }_{\mathrm{sat}}(t)=\frac{1}{G_{\mathrm{map}}}\sum _{\hat{g}\in \mathcal{G}}{o}_{\hat{g}}(t).$$

$${\eta }_{\mathrm{sat}}=\min \{{\eta }_{\mathrm{sat}}(t)\},\forall t.$$

$${\mathbf{h}}_{m,k}$$

$$=g_{m,k}\cdot {\tilde{\mathbf{h}}}_{m,k},$$

$$g_{m,k}=\sqrt{(\frac{c}{4\pi f_c{d}_{m,k}}{)}^2\cdot G_{\mathrm{dev}}{G}_{\mathrm{sat}}\cdot \frac{1}{r_0}},$$

$$d_{m,k}=-R_e\sin {\theta }_{m,k}+\sqrt{(R_e\sin {\theta }_{m,k}{)}^2+h^2+2h{R}_e},$$

$${\mathbf{h}}_{m,k}^{\text{LoS}}=\sqrt{M_a}{\mathbf{a}}_{m,k},$$

$${\mathbf{a}}_{m,k}({\phi }_{m,k},{\varrho }_{m,k})=$$

$$\exp \{j{\phi }_{m,k}\cos ({\varrho }_{m,k}-{\eta }_1)\},\cdots ,$$

$$\displaystyle\exp\{j\varphi_{m,k}\mathrm{cos}(\varrho_{m,k}-\eta_{M_{a}-1})\}\big{]},$$

$${\mathbf{Y}}_m={\alpha }_k\underset{\text{Desired signal}}{\underbrace{\sqrt{{\xi }_k}{\mathbf{s}}_k{\mathbf{h}}_{m,k}^H}}+\underset{\text{Inter-device Interference}}{\underbrace{\sum _{\begin{array}{c}{k}^{\prime }\in {\mathcal{D}}_m\\ k^{\prime }\ne k\end{array}}{\alpha }_{k^{\prime }}\sqrt{{\xi }_k}{\mathbf{s}}_{k^{\prime }}{\mathbf{h}}_{m,k^{\prime }}^H}}+\underset{\text{Noise}}{\underbrace{{\mathbf{n}}_m}},$$

$$\{\begin{array}{l}\mathbb{P}\{{\alpha }_k=1\}={\epsilon }_k\\ \mathbb{P}\{{\alpha }_k=0\}=1-{\epsilon }_k.\end{array}$$

$${\Gamma }_{m,k}=\frac{||\sqrt{{\xi }_k}{\mathbf{s}}_k{\mathbf{h}}_{m,k}^H|{|}_F^2}{I_{m,k}+{\sigma }_0^2},$$

$$I_{m,k}=\sum _{\begin{array}{c}{k}^{\prime }\in {\mathcal{D}}_m\\ k^{\prime }\ne k\end{array}}{\alpha }_{k^{\prime }}||\sqrt{{\xi }_{k^{\prime }}}{\mathbf{s}}_{k^{\prime }}{\mathbf{h}}_{m,k^{\prime }}^H|{|}_F^2.$$

$$\mathbb{E}[I_{m,k}]=\xi \epsilon L{M}_a\lambda (\frac{c}{4\pi f_c}{)}^2\cdot G_{\text{sat}}{G}_{\text{dev}}\cdot \frac{1}{r_0}2\pi R_e^2\ln (2\frac{R_e}{h}+1).$$

$$R_{m,k}=B_w{\log }_2(1+{\Gamma }_{m,k}),$$

$$\begin{aligned} \Xi=\mathbb{E}[&\log_{2}(1+\Gamma_{m,k})]=\frac{\pi(R_{e}+h)}{S_{\mathcal{A}_{k}}\ln 2\cdot R_{e}}\Big{[}\ln(1+\Psi d_{m}^{-2})d_{m}^{2}\\ &+\Psi\ln(d_{m}^{2}+\Psi)-\ln(1+\Psi h^{-2})d_{m}^{2}-\Psi\ln(h^{2}+\Psi)\Big{]}\\ \end{aligned},$$

$$d_m=-R_e\sin {\theta }_m+\sqrt{(R_e\sin {\theta }_m{)}^2+h^2+2h{R}_e}.$$

$$\mathbb{E}[R_{m,k}]=B_w\cdot \mathbb{E}[{\log }_2(1+{\Gamma }_{m,k})]=B_w\Xi .$$

$$C_k=\sum _{m\in {\mathcal{A}}_k}{R}_{m,k}.$$

$$\mathbb{E}[C_k]=\sum _{m\in {\mathcal{A}}_k}\mathbb{E}[R_{m,k}]={\bar{N}}_k^{\mathrm{sat}}\cdot \mathbb{E}[R_{m,k}].$$

$${\varpi }_{\text{manu}}=0.00185\cdot W_{\mathrm{sat}},$$

$${\varpi }_{\text{launch}}=0.000166\cdot W_{\mathrm{sat}}\cdot (h/1.609{)}^{0.43}.$$

$${\varpi }_{\text{ins}}={\beta }_{\mathrm{ins}}\cdot ({\varpi }_{\text{manu}}+{\varpi }_{\text{launch}}),$$

$${\varpi }_{\text{space}}=N\cdot P({\varpi }_{\text{manu}}+{\varpi }_{\text{launch}}+{\varpi }_{\text{ins}}).$$

$$Q1:\underset{N,P,h,i^{\mathrm{sat}}}{\mathrm{minimize}}$$

Full text

I Introduction

Nowadays, a massive number of IoT devices have been deployed all over the world for various applications, e.g., intelligent transport, smart healthcare, and factory automation, which have notably enhanced people’s quality of life [1]. Hence, International Telecommunication Union (ITU) has identified ubiquitous connectivity as one of the usage cases of 6G wireless networks [2]-[4]. However, the remote and isolated locations of some IoT devices and unstable wireless communication environment limit the traditional terrestrial cellular networks. For instance, it is difficult to construct new base stations for IoT applications in mountains, oceans, deserts and other remote areas due to the great challenge for construction and the heavy cost for maintenance. In this context, the Direct-to-Satellite (DtS) communication scheme is utilized to establish massive access111According to [5], massive access refers to that a massive number of wireless devices access the system randomly. The tremendous growth in the number of connected devices requires a fundamental rethinking of the conventional multiple access technologies in favor of new schemes suited for massive random access. Specifically, massive access can potentially support the random access millions of devices in an area of one square kilometer. for IoT devices deployed in remote areas. Specifically, DtS does not need a terrestrial gateway or other relay stations to propagate the signals, but the devices connect directly to the LEO satellite [6]- [8]. It can lead to reduced latency, improved data transfer speeds, and enhanced constellation coverage, especially in remote or geographically challenging areas [9]. Additionally, DtS communication can offer greater flexibility and mobility, making it a valuable option for applications such as maritime and aerial communication, where continuous connectivity is essential [10]. However, one LEO satellite moves very fast and thus cannot provide persistent coverage to a specific area. Therefore, it is necessary to construct a LEO satellite network to provide seamless coverage, sufficient capacity and large-scale connectivity for IoT devices.

Generally, the construction of LEO satellite networks can be divided into several stages, in which one of the most important being satellite constellation design, as it directly affects the overall system performance. A LEO satellite constellation can be regarded as a number of satellites deployed into the pre-determined orbits to form a satellite network. To provide coverage for Earth, there have been several classical constellation architectures, such as polar orbit constellation [11] and Walker Delta constellation [12]. Both of these constellations are favored for their symmetry, which plays a crucial role in ensuring efficient coverage and system performance. Generally speaking, polar orbit constellation can provide excellent global coverage. As Earth rotates underneath the orbit, the satellites can cover the entire surface of Earth over time. The Walker Delta constellation is also a specific type of satellite constellation designed to provide predictable coverage within a latitude-band region due to its symmetry. For instance, the first stage of Starlink is constructed based on Walker Delta constellation [13]. It has 24 orbit planes with 66 LEO satellites on each plane at 550 km altitude. The orbit inclination is 53∘ and the phase factor is 11. Moreover, Flower constellation is a special type of Walker delta constellation with only one satellite in each orbital plane, aimed to achieve high global coverage with a small number of satellites [14], [15]. In general, quality-of-service (QoS) requirements of terrestrial IoT devices and deployment costs of a large number of satellites are two key factors of satellite constellation design. Particularly, the QoS requirements of the satellite constellation often include coverage ratio and communication quality, while the costs of the satellite constellation usually consist of three parts: manufacturing cost, launching cost and insurance cost [16], [17]. Therefore, the design of LEO satellite constellation should balance the QoS requirements of LEO satellite IoT and the costs of the LEO satellite deployment.

In fact, the aforementioned basic constellation models have already been applied to LEO satellite IoT [6]-[10]. With the developments of IoT, the optimization of satellite constellations for various applications has received significant attention in recent years. Researchers explored various approaches to enhance the performance, coverage, and reliability of satellite constellations [20]-[28]. For instance, the authors in [20] proposed a regional satellite constellation design scheme to maximize the target area, and non-dominated sorting genetic algorithm II (NSGA-II) is adopted to address the problem. Moreover, the authors in [21] proposed a scheme for satellite constellation design by minimizing the system response time, and an optimized algorithm based on genetic algorithm (GA) was proposed for the constellation optimization. Considering the communication delay, the authors in [22] proposed a progressive construction scheme for satellite constellation design to minimize the end-to-end delay. Considering the user demands, a quality-of-experience (QoE) aware scheme for satellite constellation optimization was proposed in [23], and a multilayer tabu search (MLTS) algorithm was designed for obtaining the constellation design scheme to achieve the best QoE. In [24], the authors proposed a tabu search genetic algorithm combined with an opposite-based learning strategy to maximize the network invulnerability benefit. Building upon the need for more unified coverage models, the authors in [25] transformed the satellite constellation design into a non-negative least squares problem and derived a unified formula for coverage evaluation. Additionally, literature [26] proposed a two-step optimization method to design an asymmetric satellite constellation with resolution requirements for hurricane monitoring in specific areas. Moreover, for emergency observation, the authors in [27] investigated the performance metrics of both responsiveness and coverage to approximate the actual emergency conditions, and then proposed a multi-objective optimization framework for satellite constellation design. Furthermore, literature [28] proposed a two-stage algorithm based on reverse design for a cross-domain fusion constellation integrating communication, navigation and remote sensing. In this work, we mainly consider the coverage ratio and backhaul capacity for the QoS requirements in LEO satellite constellation design. Specifically, the coverage ratio requirement is used to ensure the coverage of the target area and backhaul capacity requirement is used to guarantee the quality of communication.

However, there are still several important issues not being addressed in the previous researches. Firstly, the aforementioned works mainly focus on the aspects of the satellite deployment but ignore the terrestrial users and devices or simply assuming that the locations of terrestrial users and devices obey the uniform distribution, which is different from the real situation. Secondly, most existing studies mainly focused on one or two of the three key design metrics—cost, communication quality, and coverage. Thirdly, the optimization techniques adopted in the existing works have low convergence speed or cannot achieve good performance, which limit the overall performance of the constellation. Motivated by this, we intend to design a stable and effective satellite constellation optimization scheme with QoS provision. Specifically, the main contributions of this paper are listed as follows:

We provide an in-depth performance analysis of LEO satellite IoT constellation. Specifically, we derive out the coverage ratio and transmission rate as a function of constellation parameters and device distribution. 2. 2.

We propose a LEO satellite IoT constellation optimization algorithm by minimizing the total deployment cost while fulfilling QoS requirements in terms of coverage ratio and communication quality. 3. 3.

We analyze the convergence and computational complexity of our proposed algorithm. Moreover, we provide extensive simulation results to verify the effectiveness of the proposed algorithm.

For clarity, we summarize the recent works on constellation design, point out their key limitations, and present the solutions of this work in Table I.

The remainder of this paper is organized as follows. In Section II, we introduce the system model of the LEO satellite IoT constellation, including the orbit elements, constellation model, coverage model and channel model. In Section III, we derive out the backhaul capacity expressions and introduce the cost model. Based on these, a cost minimization problem of LEO satellite constellation is formulated and a satellite constellation optimization algorithm is proposed. Next, extensive simulation results are provided in Section IV to verify the effectiveness of the proposed algorithm. Finally, Section V concludes the paper.

Notations: We use bold lowercase letters and bold uppercase letters to denote the column vectors and matrices, respectively, $(\cdot)^{H}$ to denote conjugate transpose, $||\cdot||_{F}^{2}$ to denote Frobenius norm of the matrix, $\mathbb{P}\{\cdot\}$ to denote the event probability, $\mathbb{E}\{\cdot\}$ to denote the expectation. $\mathcal{CN}(\mu,\sigma^{2})$ denotes the complex Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$, $\mathbf{I}_{M}$ denotes the identity matrix of dimension $M$.

II System Model

We consider the scenario of a LEO satellite IoT constellation as illustrated in Fig. 1, where a large number of LEO satellites serve a massive number of IoT devices including terrestrial devices and low altitude devices distributed over the world. To satisfy QoS requirements in terms of wireless coverage and transmission rate with a low cost, the satellite constellation architecture is supposed to be designed according to the characteristics of IoT services. In general, the LEO satellites in the constellation are operated in $P$ planes of altitude $h$, denoted by $P=\{1,2,\cdots,P\}$, wherein the orbital plane $p$ consists of $N_{p}$ LEO satellites, denoted by $\mathcal{N}_{P}=\{1,2,\cdots,N_{P}\}$. In what follows, we first introduce the satellite orbit elements and the constellation model, then present the coverage model of one LEO satellite and the whole satellite constellation, respectively. Finally, we present the channel model in LEO satellite IoT constellation.

II-A Orbit Elements and Constellation Model

Firstly, we introduce the orbit elements in the LEO satellite IoT constellation. Following previous relevant works [30]-[32], we assume that Earth is an ideal sphere and the satellite orbit is circular222Accounting for Earth’s oblateness and the uneven mass distribution, an oblate spheroid model may be more accurate for Earth shape. For ease of analysis, this paper adopts the commonly used spherical Earth model. In the follow-up researches, a practically inclined elliptical orbit can be adopted for the orbit model.. As shown in Fig. 2, six orbital elements are defined in the Earth-Centered Inertial (ECI) reference frame to describe the size, shape, and orientation of a LEO satellite’s orbit as follows:

Semi-major $a_{s}$: Half of the length of the orbit’s major axis. It determines the size of the orbit and directly affects the satellite’s altitude. In circular orbits, it is equal to the radius of the orbit, i.e., $a_{s}=R_{e}+h$, where $R_{e}$ is the radius of Earth and $h$ denotes the satellite altitude. 2. 2.

Eccentricity $e$: The eccentricity of the orbital ellipse. It determines the shape of the orbit. Given the circular orbit, the eccentricity $e$ is 0. 3. 3.

Inclination angle $i$: Angle measured counterclockwise from the equatorial plane to the satellite orbital plane at the ascending node. It describes the tilt of the orbit relative to Earth’s equatorial plane. An inclination angle of more than $90^{\circ}$ indicates a retrograde orbit, where the satellite’s motion is opposite to that of Earth’s rotation. 4. 4.

Right ascension of ascending node $\Omega$: Angle measured counterclockwise from the vernal equinox node to the ascending node in the equatorial plane. It determines the orientation of the orbital plane in space. 5. 5.

Argument of the perigee $\omega$: Angle measured counterclockwise from the ascending node to the perigee, i.e., the point when the satellite operates closest to Earth, measured along the orbital plane. It specifies the location of perigee within the orbit. 6. 6.

True anomaly $\upsilon$: Angle measured counterclockwise from the perigee to the satellite position. It determines the satellite’s position along its orbit at a given time.

With these six orbit elements, the position of each satellite can be described accurately. Mathematically, the Cartesian coordinate of each satellite can be expressed as [34] 333This model is valid only under the assumption of a spherical Earth.

[TABLE]

For the constellation model, the specific topologies of classical constellations, e.g., the polar orbit constellation, the Walker-Delta constellation, and the Flower constellation, have been extensively studied for their capabilities to address particular coverage requirements while maintaining symmetries and making assumptions that aid analysis, deployment, and operation [35]. Moreover, they are used as parameter sets in both coverage analysis and network performance analysis [7], [37]. Therefore, we adopt Walker-Delta constellation as the constellation model for optimization in this work.

Specifically, for the Walker-Delta constellation, it can be described with three parameters $[N,P,F]$, where $N$, $P$ and $F$ denote the number of satellites per satellite plane, the number of satellite planes and the phase factor, respectively. For the phase factor $F$, it takes a value between 0 to $(P-1)$, which determines the phase difference $\nabla u$ between the corresponding satellites in two adjacent satellite planes, i.e.,

[TABLE]

II-B Coverage Model

Wireless coverage is a key factor in the design of LEO satellite IoT constellation. As illustrated in Fig. 3, the coverage model of a LEO satellite is determined by the altitude $h$ and the minimum elevation angle $\theta_{\min}$444 The coverage model used in this work assumes line-of-sight above a certain elevation threshold.. In this way, the angular radius $\phi_{\max}$, i.e., the angle from the point with the minimum elevation angle $\theta_{\min}$ to the projection point of the satellite, can be expressed as

[TABLE]

Similarly, the satellite coverage angle $\beta$ can be computed as

[TABLE]

For a given angular radius $\phi_{\max}$, the coverage area of a LEO satellite can be calculated as

[TABLE]

Therefore, the satellite at a higher altitude $h$ can provide wider coverage, while a larger elevation angle $\theta_{\min}$ implies the decrease of coverage.

In above, we obtained the coverage area of a LEO satellite. Actually, we are more concerned about the coverage ratio of the LEO satellite IoT constellation. To calculate the coverage ratio $\eta_{\mathrm{sat}}$ of the target area served by the LEO satellite IoT constellation, the grid method is adopted as illustrated in Fig. 4. Specifically, it grids the target area into observation points $\mathcal{G}=\{1,2,\cdots,G_{\mathrm{map}}\}$. Moreover, the total observation time $T_{\mathrm{total}}$ is divided into multiple time slots. Then, we use the binary variables $o_{\hat{g}}(t)$ to denote whether the observation grid $\hat{g}\in\mathcal{G}$ is covered or not by the satellite IoT constellation at the time slot $t\in\mathcal{T}=\{1,2,\cdots,T_{\mathrm{total}}\}$, i.e.,

[TABLE]

Based on the location of a satellite under Cartesian coordinates in (1) and the coverage area of a satellite in (5), we can indicate whether an observation grid is covered as conducted in (6). Then, the coverage ratio of a LEO satellite IoT constellation in time slot $t$ can be calculated as

[TABLE]

Without loss of generality, we use $\eta_{\mathrm{sat}}$ to denote the minimum coverage ratio of the LEO satellite IoT constellation in all time slots, i.e.,

[TABLE]

Based on the above analysis, the calculation of coverage ratio is related to the grid division precision and time division. Specifically, the coverage ratio is more accurate with higher grid precision and shorter time steps. However, increasing grid precision and shortening time steps inevitably increase computation time [36]. Thus, it makes sense to select the proper grid precision and time steps. Fortunately, Satellite Tool Kit (STK) is a helpful tool for us to get the precise numerical results including the coverage ratio for satellite constellations.

II-C Channel Model

According to the signal propagation characteristics of LEO satellite communications [38], the LEO satellite channel usually includes two components, namely line of sight (LoS) and non-line of sight (NLoS), as shown in Fig. 5. To be more specific, the channel from the device $k$ to the satellite $m$ can be expressed as555In this paper, the Doppler frequency shift caused by the motion of LEO satellites is pre-compensated in advance. As the velocity of the satellite can be known from the Ephemeris information, the Doppler frequency shift caused by satellites can be calculated and compensated according to [39].

[TABLE]

where $g_{m,k}$ denotes the channel large-scale fading factor, which can be calculated as

[TABLE]

where $(\frac{c}{4\pi f_{c}d_{m,k}})^{2}$ is the free path loss with $c$ being the light speed, $f_{c}$ being the carrier frequency, and $d_{m,k}$ being the distance between the satellite $m$ and the device $k$ given by

[TABLE]

where $\theta_{m,k}$ denotes the elevation angle of the device $k$ with respect to the satellite $m$. Moreover, $G_{\mathrm{dev}}$ denotes the transmit antenna gain of IoT devices, $G_{\mathrm{sat}}$ denotes the receive antenna gain of satellites, $r_{0}$ denotes the average rain attenuation coefficient. For the small-scale channel gain $\tilde{\mathbf{h}}_{m,k}$, $\lambda_{m,k}$ denotes the Rician factor of the channel between the satellite $m$ and the device $k$. The LOS component $\mathbf{h}^{\text{LoS}}_{m,k}$ can be modeled as

[TABLE]

where $M_{a}$ denotes the number of satellite antennas, and $\mathbf{a}_{m,k}$ is the array responding vector of antennas. Without loss of generality, each LEO satellite is equipped with a uniform circular array which is widely adopted in current LEO satellite constellations. Thus, the array responding vector $\mathbf{a}_{m,k}$ can be expressed as

[TABLE]

where $\eta_{j}=2\pi j/M_{a},j=0,1,\cdots,M_{a}$, $\varphi_{m,k}=\frac{\pi d_{s}f_{c}}{c}\sin\vartheta_{m,k}$ with $d_{s}$ being the diameter of circular antenna array equipped on the satellite and $\vartheta_{m,k}$ being the off-axis of the satellite $m$ boresight to the device $k$, and $\varrho_{m,k}$ is the azimuth angle of the device $k$ to the array center of satellite $m$. In addition, $\mathbf{h}^{\text{NLoS}}_{m,k}$ represents the NLoS component of the LEO satellite channel which follows the complex Gaussian distribution, i.e., $\mathbf{h}^{\text{NLoS}}_{m,k}\sim\mathcal{CN}(\mathbf{0},v_{m,k}^{\mathrm{NLoS}}\mathbf{I}_{M_{a}})$ with $v_{m,k}^{\mathrm{NLoS}}$ being the variance.

III LEO Satellite Constellation Design

In this section, we first introduce the communication model and then develop a cost model for the LEO satellite IoT constellation. Based on these, we formulate a constellation optimization problem by minimizing the cost while providing QoS guarantee in terms of communication quality and coverage ratio, and thereby derive an efficient constellation design algorithm.

III-A Communication Model for LEO Satellite IoT

In LEO satellite IoT, due to the bursty characteristics of IoT services, only a small portion of IoT devices are active. In a certain time slot, active IoT devices simultaneously transmit their data sequences of length $L$ over the same time-frequency resource block [40], [41]. Thus, the received signal $\mathbf{Y}_{m}\in\mathbb{C}^{L\times M_{a}}$ at the satellite $m$ can be expressed as

[TABLE]

where $\mathcal{D}_{m}$ denotes the interference area of the satellite $m$, $\alpha_{k}$ is the activation indicator of the $k$-th device with $\alpha_{k}=1$ if the device $k$ is active, otherwise, $\alpha_{k}=0$. Let $\varepsilon_{k}$ represent the active probability of the device $k$, then we have

[TABLE]

Moreover, $\xi_{k}$ denotes the transmit power of the device $k$, $\mathbf{s}_{k}$ is the transmitted data sequence, which is Gaussian distributed with unit norm, and $\mathbf{n}_{m}$ denotes the additive white Gaussian noise (AWGN) with variance $\sigma_{0}^{2}$. Hence, the signal to interference plus noise ratio (SINR) for the link between satellite $m$ and device $k$ can be computed as

[TABLE]

where $I_{m,k}$ is the interference for the link, which can be calculated as

[TABLE]

In the satellite constellation design, we are more concerned with the overall performance of the system rather than that of a single device or single satellite. Hence, we analyze the overall performance according to the distribution of IoT devices. Without loss of generality, the IoT devices are assumed to be distributed as a homogeneous Poisson point process (HPPP) $\Phi$ with device density $\lambda$ devices/km2 [42]. Based on this, we can obtain the following theorem:

*Theorem 1: * The average interference of the link between the device $k$ and the satellite $m$ can be calculated as

[TABLE]

Proof:

See Appendix A. ∎

Further, with the SINR in (16), the transmission rate for the link from the device $k$ to the satellite $m$ is given by

[TABLE]

where $B_{w}$ is the channel bandwidth. To calculate the average data rate $\mathbb{E}[R_{m,k}]$, we provide the following theorem:

*Theorem 2: * The term $\mathbb{E}[\log_{2}(1+\Gamma_{m,k})]$ can be computed as

[TABLE]

where $d_{m}$ denotes the maximum communication distance of satellite $m$. According to (11), $d_{m}$ can be given by

[TABLE]

Proof:

See Appendix B. ∎

Thereby, the average transmission rate for the link between device $k$ and satellite $m$ can be given by

[TABLE]

In the LEO satellite IoT constellation, one IoT device may be served by multiple satellites due to multiple-satellite coverage. To calculate the total backhaul capacity of device $k$, we define the set of serving satellites for the device $k$ as $\mathcal{A}_{k}=\{m|\theta_{m,k}\geq\theta_{\min}\}$, where $\theta_{m,k}$ is the elevation angle from the device $k$ to the satellite $m$. Therefore, the total backhaul capacity of the device $k$ can be calculated as

[TABLE]

It is assumed that the average number of satellites in $\mathcal{A}_{k}$ is $\bar{N}^{\mathrm{sat}}_{k}$, then the average total backhaul capacity $\mathbb{E}[C_{k}]$ can be calculated as

[TABLE]

III-B Cost Model

The construction cost of LEO satellite constellations is crucial for enabling efficient planning, budgeting, and the promotion of sustainable operations. In general, the costs associated with the space segment significantly contribute to the overall expenses of LEO satellite constellations. Therefore, we first analyze the composition of the costs of the space segment [16]-[19].666We adopt the cost model proposed in [16] and [17] and update the coefficients of the model according to [18] and [19]. Notice that the proposed cost model is an approximation rather than a generalized rule, the coefficients in the model can be replaced according to the costs of different satellite companies. Specifically, the costs of the space segment can be divided into manufacturing cost, launch cost, and insurance cost, respectively. For the manufacturing cost of a LEO satellite, it can be calculated as

[TABLE]

where $W_{\mathrm{sat}}$ denotes the weight of a LEO satellite. For the launch cost, it can be modeled as

[TABLE]

Finally, for the insurance cost, it can be estimated based on the manufacturing cost and launch cost as

[TABLE]

where $\beta_{\mathrm{ins}}$ denotes the insurance ratio factor. Therefore, the total costs of space segment for a LEO satellite constellation are the summation of these three parts as follows

[TABLE]

It is clear that the costs of the satellite constellation are influenced by the number of LEO satellites and the altitude of the deployed constellation. In this context, we aim to design a LEO satellite IoT constellation by minimizing the total cost under QoS constraints.

III-C Problem Formulation

To balance the construction cost and the practical utilization of the desired LEO satellite IoT constellation, we formulate a cost minimization problem while satisfying the QoS requirements under the Walk-Delta architecture. Particularly, the design of LEO satellite IoT constellation can be formulated as the following optimization problem

[TABLE]

where $C_{k}^{t}=\mathbb{E}[C_{k}](t)$ denotes the average backhaul capacity $\mathbb{E}[C_{k}]$ in time slot $t$, $C_{th}$ denotes the backhaul capacity threshold, $\eta_{\mathrm{sat}}$ denotes the coverage ratio, and $\eta_{th}$ denotes the coverage ratio requirement. Specifically, the constraint (29g) is the coverage ratio constraint, the constraint (29g) is the communication quality constraint, and the constraints (29g)-(29g) are optimization ranges for the constellation parameters of constellation altitude $h$, number of orbit planes $P$, number of satellites per plane $N$ and inclination angle $i^{\mathrm{sat}}$, respectively. Moreover, as phase factor $F$ is mainly used to avoid satellite collision, the phase factor is not optimized in the problem [43].

For such an optimization problem, it is difficult to handle the constraints (29g) and (29g) directly due to the coverage probability and coverage repetitions. Moreover, these terms are intractable to get a closed-form solution. To this end, according to (24) and constraint (29g), the minimum serving satellites $\bar{N}^{\mathrm{sat}}_{th}$ for the device $k$ can be calculated as

[TABLE]

Similarly, the minimum serving satellites $\bar{N}^{\mathrm{sat}}_{\min}$ can be obtained by measuring the coverage ratio via STK [20]. Then the constraint (29g) is transformed to satisfy $\bar{N}^{\mathrm{sat}}_{\min}\geq\bar{N}^{\mathrm{sat}}_{th}$. Since the formulated problem is a non-linear, non-convex, and discrete complex combined optimization problem, classical optimization techniques based on continuous convex function cannot be employed directly. Fortunately, intelligent optimization techniques can be applied to handle these problems for its low computational complexity and fast convergence speed. Therefore, we develop an intelligent optimization algorithm based on a GA framework to solve the above optimization.

III-D Algorithm Design

In this subsection, we design a constellation optimization algorithm based on the GA framework to handle the problem. In general, the GA framework includes the following steps: 1) Chromosomes coding, which constructs a data structure for the optimization parameters. 2) Population initialization, which initializes the chromosomes for the initial population, a good method of population initialization can accelerate the convergence speed. 3) Fitness function design, which determines the direction of optimization. 4) Operator design, which generates possible solutions for the best fitness value. 5) Next generation selection, which chooses the best chromosomes to form the next population for the algorithm. In the following, we will introduce these steps in detail.

1) Chromosomes Coding: Upon the GA framework, the first step is to code the chromosomes. According to the characteristics of the formulated optimization problem, the coding vector for a chromosome is set as

[TABLE]

In addition, we adopt the real coding for the chromosome coding.

2) Population Initialization: For the initial population generation, we use the random generation method to generate the initial population. Specifically, the chromosomes for the initial population are randomly generated between the low bound $\mathbf{l}=[h_{\min},P_{\min},N_{\min},i^{\mathrm{sat}}_{\min}]$ and the upper bound $\mathbf{u}=[h_{\max},P_{\max},N_{\max},i^{\mathrm{sat}}_{\max}]$.

3) Fitness Function Design: To evaluate the quality of potential solutions, it is necessary to design a fitness function based on both objective and constraints. For addressing constraint, the penalty function method is widely applied, with a prevalent model being the amalgamation of the original objective function and the penalty function, which can be expressed as

[TABLE]

where $f(x)$ is the objective function, $p(x)$ is the constraint function, and $r$ is the penalty factor. It can be easily found in (32) that choosing a proper penalty factor for the problem is important. Specifically, if the penalty factor is too large, the problem will fail into a local optimal solution. Else if the penalty factor is too small, the constraints may not be satisfied. To overcome this challenge, we propose a novel penalty function design method. Firstly, we introduce the concepts of objective satisfaction function and constraint satisfaction function as below:

Definition 1: The objective satisfaction function $f_{s}(x)$ is defined as

[TABLE]

where $f_{\mathrm{max}}$ and $f_{\mathrm{min}}$ are the maximum value and minimum value of $f(x)$, respectively.

Definition 2: The constraint satisfaction function $p_{s}(x)$ is defined as

[TABLE]

where $p_{\mathrm{max}}$ is the maximum value of $p(x)$.

With these two satisfaction functions, we can evaluate the quality of the solutions. In particular, we define a new penalty function based on the two kinds of satisfaction functions as

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are hyperparameters to adjust the level of punishments, $m_{\mathrm{fe}}$ and $N_{\mathrm{totp}}$ denote the number of infeasible solutions and total solutions of genetic algorithm, respectively, and $n_{\mathrm{it}}$ is the current number of evolution. Comparing the proposed model (35) with the classical model (32), we adjust the penalty factor by utilizing the information of the current population rather than a statistic constant. Specifically, we introduce the ratio of the number of the infeasible solutions to the total solutions $\frac{m_{\mathrm{fe}}}{N_{\mathrm{totp}}}$ in the penalty factor. When there are more infeasible solutions, the model will increase the penalties. If there are more feasible solutions, the model will reduce the penalties and correspondingly increase the value of objective function. Hence, the penalty factor is readily manageable. More importantly, the designed penalty factor promotes the search towards a better solution, leading to a faster convergence speed. However, if we adopt the classical model in (32) for penalty, it is difficult to distinguish the feasible solutions and the infeasible ones, but the proposed penalty model can address this challenge by utilizing the product form of two satisfaction functions.

To apply the proposed penalty function in problem $Q1$, the objective function $f_{1}(\mathbf{x})$ and constraint functions $p_{1}(\mathbf{x})$ and $p_{2}(\mathbf{x})$ should be written as

[TABLE]

[TABLE]

[TABLE]

where $\Phi(x)=\max\{x,0\}$ is used to measure the degree of violation of the constraints. To calculate the corresponding satisfaction functions, we introduce $f_{1}^{\min}$, $f_{1}^{\max}$, $p_{1}^{\max}$, and $p_{2}^{\max}$ to denote the minimum and maximum values of objective function and constraint functions (36) - (38) in the population, respectively. As a result, the satisfaction functions of (36)-(38) can be constructed as

[TABLE]

[TABLE]

[TABLE]

In this way, the proposed penalty function for problem $Q1$ can be expressed as

[TABLE]

which serves as the fitness function to evaluate the quality of potential solution.

4) Operation Design: The fourth step of the GA framework is to design the GA operators, including crossover operator and mutation operator. To improve the ability of the algorithm to find optimal solutions, we add the information of the optimal individual in the crossover operator. Specifically, we use $\mathbf{x}_{\mathrm{best}}$ to denote the best individual in the $n_{\mathrm{it}}$-th population $\mathcal{P}(n_{\mathrm{it}})$. Then, we choose $N_{1}$ individuals according to roulette wheel selection with their fitness to form a parent set $\mathcal{P}_{pa}$. For the parents $\mathbf{x}_{p,1}$ and $\mathbf{x}_{p,2}$, two new individuals are breed by the following crossover operator

[TABLE]

[TABLE]

where $r_{ij}(i=1,2;j=1,2,3)$ are random values in the range 0 to 1, $\mathbf{x}_{c,1}$ and $\mathbf{x}_{c,2}\in\mathcal{P}_{1}(n_{\mathrm{it}})$ are the child solutions and $\mathcal{P}_{1}(n_{\mathrm{it}})$ is the child generation set.

For the mutation operator, the algorithm is expected to have a relatively large search range at the beginning of the evolution to preserve the diversity of the population, while at the later stages of the evolution it should be guaranteed that there are enough feasible solutions for the algorithm to converge. Motivated by this, we design the mutation operator with a mutation indicator $\theta_{\mathrm{mut}}\in[0,1]$ as follows: if $\theta_{\mathrm{mut}}\leq(\frac{1}{n_{\mathrm{it}}+1}+1)\varsigma$, the chromosome after mutation $\mathbf{x}_{mu}$ is determined by

[TABLE]

where $\varsigma$ is the threshold of mutation operator, and $\Delta\mathbf{m}=(\Delta m_{1},\Delta m_{2},\cdots,\Delta m_{n})^{T}$ with $\Delta m_{i}$ following the Gaussian distribution with zero mean and $\sigma_{i}^{2}$, i.e., $\Delta m_{i}\sim\mathcal{N}(0,\sigma_{i}^{2})$. Otherwise, if $\theta_{\mathrm{mut}}>(\frac{1}{n_{\mathrm{it}}+1}+1)\varsigma$, the chromosome after mutation $\mathbf{x}_{mu}$ is set as

[TABLE]

where $r_{3,i}$ and $r_{4,i}$ are random numbers range in $[0,1]$, $(x_{mu})_{i}$ and $x_{i}$ respectively denotes the $i$-th element of $\mathbf{x}_{mu}$ and $\mathbf{x}$, $u_{i}$ denotes the $i$-th element of upper bound of search space $\mathbf{u}=[h_{\max},P_{\max},N_{\max},i^{\mathrm{sat}}_{\max}]$, and $l_{i}$ denotes the $i$-th element of lower bound of search space $\mathbf{l}=[h_{\min},P_{\min},N_{\min},i^{\mathrm{sat}}_{\min}]$. Consequently, at the beginning of the evolution, the population will be more likely to mutate according to (45). While at the end of the evolution, the population will be more likely to mutate according to (46) as to perform local search around the mutated chromosome. For convenience, we use $\mathcal{P}_{2}(n_{\mathrm{it}})$ to denote the population after mutation.

5) Next Generation Selection: Finally, we adopt the elite preservation method to form the population of next generation. In other words, we choose the best $N_{\mathrm{totp}}$ individuals in the set $\mathcal{P}(n_{\mathrm{it}})\cup\mathcal{P}_{1}(n_{\mathrm{it}})\cup\mathcal{P}_{2}(n_{\mathrm{it}})$ to form the $\mathcal{P}(n_{\mathrm{it}}+1)$. The whole algorithm is summarized in Algorithm 1. For clarity, we describe the algorithm flow in Fig. 6.

 

Algorithm 1 : Design of LEO Satellite IoT Constellation

 

1: Input: $f_{c}$, $G_{\mathrm{dev}}$, $G_{\mathrm{sat}}$, $B_{w}$, $T$, $r_{0}$, $\beta$, $M$, $\varepsilon_{k}$, $L$, $\lambda$, $W_{\mathrm{sat}}$, $\eta_{th}$, $C_{th}$, $\mathbf{u}$, $\mathbf{l}$, $N_{\mathrm{totp}}$, $n_{\mathrm{it}}$, $\xi_{k}$

2: Initialization: generate $\mathcal{P}(1)=[\mathbf{x}_{1},~\mathbf{x}_{2},~\cdots,~\mathbf{x}_{N_{\mathrm{totp}}}]$ randomly with upper bound $\mathbf{u}$ and lower bound $\mathbf{l}$, iteration index $n_{\mathrm{it}}=1$;

3:  While $n_{\mathrm{it}}\leq N_{\mathrm{it}}$

4: Calculate the fitness value of individuals in $\mathcal{P}(n_{\mathrm{it}})$ according to (36) - (42);

5: Parents Selection: Choose the best $N_{\mathrm{totp}}$ individuals in $\mathcal{P}(n_{\mathrm{it}})$ by roulette wheel selection with the fitness function calculated by step 4 to form the parents’ set $\mathcal{P}_{pa}(n_{\mathrm{it}})$.

6: Crossover: Randomly choose two individuals in $\mathcal{P}_{pa}(n_{\mathrm{it}})$ and the best individual to perform the crossover operators according to (43) and (44) until getting $N_{\mathrm{totp}}$ children to form the population $\mathcal{P}_{1}(n_{\mathrm{it}})$;  

7: Mutation: For each individual individual in $\mathcal{P}_{1}$, perform the mutation operator as: For $s=1:N_{\mathrm{totp}}$:      Randomly choose $\theta_{\mathrm{mut}}\in[0,~1]$;      If $\theta_{\mathrm{mut}}\leq(\frac{1}{n_{\mathrm{it}}+1}+1)\varsigma$         perform the mutation operator for $\mathbf{x}(s)$ according to (45);      Else if $\theta_{\mathrm{mut}}>(\frac{1}{n_{\mathrm{it}}+1}+1)\varsigma$         perform the mutation operator for $\mathbf{x}(s)$ according to (46); The individuals after mutation form the population $P_{2}(n_{\mathrm{it}})$;

8: Next Generation Selection: Calculate the fitness function in the population $\mathcal{P}(n_{\mathrm{it}})\cup\mathcal{P}_{1}(n_{\mathrm{it}})\cup\mathcal{P}_{2}(n_{\mathrm{it}})$, and select the best $N_{1}$ individuals to form the population $\mathcal{P}(n_{\mathrm{it}}+1)$;

9: $n_{\mathrm{it}}=n_{\mathrm{it}}+1$;

10: End while;

11: Output: the minimum space costs $\varpi_{\mathrm{space}}^{\mathrm{opt}}$, the optimized number of satellites per plane $N_{\mathrm{opt}}$, the optimized number of satellite planes $P_{\mathrm{opt}}$, the optimized altitude $h_{\mathrm{opt}}$, and the optimized inclination angle $i^{\mathrm{sat}}_{\mathrm{opt}}$.

III-E Algorithm Analysis

To gain further insights of the proposed algorithm, we analyze its convergence properties and computational complexity.

Convergence Analysis: The proposed algorithm has fast convergence speed due to the special design even with very large constellation parameters. The proof of convergence please refers to Appendix C.

Complexity Analysis: The proposed algorithm is a kind of heuristic algorithm, which trades off the computational precision for reduced complexity. Moreover, for the candidate solutions generated by the algorithm, parallel computing techniques are adopted to reduce the computational complexity. For each iteration of our proposed algorithm, the main computational complexity is caused by the calculation of coverage ratio by STK. Therefore, the complexity of our proposed algorithm is $O(N_{\mathrm{it}}N_{\mathrm{totp}})$. To verify the effectiveness of the proposed algorithm, we compare it with five typical algorithms which are widely adopted in satellite constellation optimization in Table II. The comparison results show that our proposed algorithm does not increase the computational complexity.

IV Simulation Results

In this section, we provide extensive simulation results to testify the performance of the proposed algorithm for LEO satellite IoT constellation design. According to TR 38.821 [46], the main simulation parameters of the satellite constellation are summarized in Table III. Moreover, the range of optimization parameters are set in Table IV.

To visualize the coverage of constellations, Fig. 7 has provided the two-dimensional (2D) simulation result of the proposed LEO satellite IoT constellation design algorithm under the default simulation parameters of Table III and IV. Specifically, the optimized number of planes is $P_{\mathrm{opt}}=6$, the optimized number of satellites per plane is $N_{\mathrm{opt}}=8$, the optimized altitude is $h_{\mathrm{opt}}=1589$ km, and the optimized inclination angle $i^{\mathrm{sat}}_{\mathrm{opt}}=41^{\circ}$. It can be seen from Fig. 7 that the proposed LEO satellite IoT constellation design algorithm has provided sufficient coverage for the target area. Moreover, it is presented in Fig. 8 that the three-dimensional (3D) model of the wireless coverage of the designed LEO satellite IoT constellation, which also confirms the strong coverage performance of the proposed algorithm.

Furthermore, to verify the convergence and effectiveness of the proposed algorithm, we compare it with five baseline intelligent algorithms, which are widely adopted in satellite constellation optimization, including classical GA [47], Sine Cosine algorithm (SCA) [48], Grey Wolf Optimization (GWO) [49], Particle Swarm optimization (PSO) [50], and Tabu Search (TS) algorithm [23]. For the classical GA in [47], the crossover operation adopted in classical GA does not take into account information about the current optimal solution, and the mutation operator remains unchanged during the iterations. The simulation parameters of the proposed algorithm and classical GA are listed in Table V. Moreover, SCA and GWO are two novel swarm intelligent optimization algorithms. For PSO, it is a widely adopted optimization algorithm in which each particle retains memory of its previous best position, making it applicable for many optimization problems. Finally, TS is a metaheuristic local search algorithm by checking its immediate neighbors of the current solution in search of finding an improved solution. To avoid contingency, statistical analysis is performed for each algorithm. Specifically, the average total cost of each algorithm is evaluated for 200 trials. As the penalty function in (42) contains several hyperparameters of GA algorithm, it is not suitable for the baseline algorithms to utilize (42) as the fitness function. Therefore, a classical penalty function is adopted for the baseline algorithms as the fitness function, which is expressed as

[TABLE]

where $\rho_{1}$ and $\rho_{2}$ are two fixed penal factors. In Fig. 9, we show the convergence performance of the proposed algorithm alongside five baseline intelligent algorithms. Evidently, the proposed algorithm demonstrates the fastest convergence among all the algorithms. Noticeably, due to the special design of the fitness function and algorithm operators, the proposed algorithm has the lowest cost. Therefore, the proposed algorithm is effective in LEO satellite IoT constellation optimization.

Moreover, it is essential to verify whether the proposed algorithm satisfies the both QoS requirements. Hence, we investigate the coverage ratio performance of the proposed algorithm under different backhaul capacity requirements. As illustrated in Fig. 10, the minimum coverage ratio of the proposed algorithm exceeds the desired coverage ratio for a given backhaul capacity. Besides, the average coverage ratio over all time slots grows as the backhaul capacity requirements $C_{th}$ increases. This is because a higher backhaul capacity requirement for IoT devices always needs more LEO satellites to provide coverage for a specific area. Particularly, the simulation results at 30 Mbps and 40 Mbps demonstrate identical performance, indicating that the proposed constellation is the most cost-effective while meeting the coverage requirements. Additionally, it is illustrated in Fig. 11 that the proposed algorithm can satisfy the backhaul capacity requirements under different coverage ratio requirements. Therefore, the proposed algorithm can provide high QoS for the target area.

Subsequently, we explore the relationship between the total cost and the QoS requirements. As depicted in Fig. 12, a higher backhaul capacity requirement or coverage ratio requirement will increase the total cost of the LEO satellite IoT constellation. Furthermore, it is also found that the total cost keeps nearly unchanged when $C_{th}$ increases from $50$ Mbps to $60$ Mbps at the coverage ratio requirement of $100\%$, since the constellation with seamless coverage has achieved the backhaul capacity more than 60 Mbps. In this way, we should balance the total cost and the QoS requirements when designing the LEO satellite IoT constellation.

LEO satellite IoT constellation needs to provide service for diverse IoT devices in different scenarios, e.g., ocean, desert, forests and mountains. Generally, these IoT devices may have quite different device densities and activity probabilities. For this purpose, we investigate the relation between the total cost of the LEO satellite constellation and the IoT device densities and activity probabilities. It can be seen in Fig. 13 that for a given activity probability $\varepsilon$, as the device density increases, the total cost increases. Moreover, the total cost is added as the activity probability increases. This is because more activated IoT devices need to be supported by LEO satellites in a certain area, resulting in the increment of the number of LEO satellites in a satellite constellation. Additionally, when the product of device density and activity factor is equal, the total cost of the designed constellation has the same value, since the value of the activity factor times the device density is directly related to the average interference and backhaul capacity, as indicated in (19) and (20).

Generally, the coverage performance is essential for the design of LEO satellite constellation. Therefore, the relationship between the total cost of the satellite constellation and the coverage ratio requirements has been shown in Fig. 14 under different coverage areas. It can be found that $100\%$ coverage for the target area will result in a sharp increment in costs. Moreover, it is worth noticing that the proposed algorithm is applicable in the scenarios with different coverage areas. Specifically, for a given required coverage ratio, the cost is added as the coverage area is broadened.

For transmit power of IoT devices, it will also affect the performance of the LEO satellite IoT constellation. Fig. 15 shows that increasing the transmit power of IoT devices can effectively reduce the total cost of the constellation. However, the total cost will come to a convergence as the transmit power of devices $\xi$ increases. Fortunately, as illustrated in Fig. 15, a feasible way to reduce the total cost of LEO satellite IoT constellation is to equip the satellite with sensor antennas which has large coverage angles. In this way, designing sensor antennas with large coverage angles for LEO satellite is an important task in the future of satellite constellation.

Finally, we compare the proposed LEO satellite constellation with three classical LEO satellite constellations under identical communication backhaul capacity and coverage ratio requirements [51]-[53]. It can be seen from Fig. 16 that the proposed LEO satellite constellation is cost-competitive compared to other constellations, meeting our desired objectives. Furthermore, it is worth pointing out that the results in Fig. 16 do not mean that other constellations are useless. Actually, all the constellations are designed for different tasks. Specifically, Global Star is to provide seamless coverage for users except North and South Poles. Iridium aims to provide global communication service. Telesat wants to provide cost-effective, fiber-like, global network connectivity for business, government and individual users. Fortunately, the proposed algorithm can minimize the deployment cost while guaranteeing the coverage ratio and communication quality.

V Conclusion and Future Work

In this paper, we provided a QoS-driven satellite constellation design scheme for LEO satellite IoT. In particular, a LEO satellite IoT constellation optimization algorithm was proposed to minimize the deployment cost while guaranteeing the coverage ratio and communication quality of IoT devices. Theoretical analysis proved that the proposed algorithm can achieve the global convergence. Finally, extensive simulation results validated the effectiveness of the proposed algorithm. Specifically, the proposed algorithm can converge in $10$ iterations and obtain a lower satellite constellation cost than the baseline optimization algorithms.

As for future work, the paper suggests several extensions. These include integrating dynamic traffic demand models, exploring multi-objective optimization with real-time adaptive algorithms, and investigating the impact of inter-satellite links on constellation performance. These extensions would provide valuable insights for next-generation LEO satellite constellation design.

Appendix A The Proof of Theorem 1

Notice that the activation indicator $\alpha_{k}$, data sequence $\mathbf{s}_{k}$, large-scale fading factor $g_{m,k}$ and small-scale fading factor $\widetilde{\mathbf{h}}_{m,k}$ are mutually independent, such that the average interference $\mathbb{E}[I_{m,k}]$ can be decomposed into a product of the expectations of these variables as

[TABLE]

where $\mathbf{s}_{k}$ denotes the transmitted data sequence of device $k$, $L$ and $M_{a}$ represents the length of data sequence and the number of antennas, respectively. As the distribution of IoT devices follows HPPP with density $\lambda$, according to the Campbell’s theorem [54], the average interference $\mathbb{E}[I_{m,k}]$ can be derived as

[TABLE]

The area of $\mathcal{D}_{m,k}$ is a spherical crown with the angular radius $\varphi_{0}$, which motivates us to use the spherical segments as the partition of the integral region, then the average interference $\mathbb{E}[I_{m,k}]$ can be obtained as

[TABLE]

The proof is completed.

Appendix B The Proof of Theorem 2

For the expectation of term $\log_{2}(1+\Gamma_{m,k})$ in (19), it can be further expressed as

[TABLE]

where $\xi$ and $\varepsilon$ denote the average transmit power and average activity probability over all IoT devices, respectively, $S_{\mathcal{A}_{k}}$ denotes the area of visible satellite region of device $k$. To make the expression more clearly, we use $\Psi$ to denote the coefficient of the integral variable, i.e., $\Psi=\frac{\xi LM_{a}\cdot(\frac{c}{4\pi f_{c}})^{2}\cdot G_{\text{sat}}G_{\text{dev}}\cdot\frac{1}{r_{0}}}{\mathbb{E}[I_{m,k}]+\sigma_{0}^{2}}$. Then, the expectation can be rewritten as

[TABLE]

where in (a) the substitution $u=R_{e}^{2}+(R_{e}+h)^{2}-2R_{e}(R_{e}+h)\cos v$ is applied.

The proof is completed.

Appendix C The Proof of the global convergence of proposed algorithm

For the sake of convenience in the proof, we present the following important concepts.

Definition 3: Let $\{X_{k}:k\geq 1\}$ denote a sequence of random variables defined on a sample space $\Omega$. We say that $\{X_{k}:k\geq 1\}$ is convergent in probability to a random variable $X$ if and only if

[TABLE]

for any $\delta>0$.

Lemma 1 (Borel–Cantelli lemma, BCL): Suppose $\{A_{k}:k\geq 1\}$ denotes a sequence of random variables defined on a sample space $\Omega$. For simplicity, we let $p_{k}=\mathbb{P}(A_{k})$, if $\sum_{k=1}^{\infty}p_{k}<\infty$, then

[TABLE]

else if $\sum_{k=1}^{\infty}p_{k}=\infty$ and $\{A_{k}:k\geq 1\}$ are mutual independent, then

[TABLE]

The detailed proof of Borel–Cantelli lemma can be found in [45].

For $\forall\delta>0$, we define $Q_{1}$ which satisfies

[TABLE]

where $D$ denotes the feasible domain of the optimization problem, and $f_{\min}(x)=\mathrm{\min}\{f(x),x\in D\}$. Moreover, we define $B_{2}$ as the complementary set of $B_{1}$, i.e.,

[TABLE]

In this way, the population $\{P(t)\}$ processed by the proposed algorithm can be divided two states: $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$. If there exists one individual of $\{P(t)\}$ belonging to $B_{1}$, then it means that the population $P(t)$ is in state $\mathcal{S}_{1}$, i.e., $P(t)\in\mathcal{S}_{1}$. Else if all the individuals of $P(t)$ belong to $\mathcal{B}_{2}$, then it means that the population $\{P(t)\}$ is in state $\mathcal{S}_{2}$, i.e., $P(t)\in\mathcal{S}_{2}$. Based on it, we can derive the following lemma:

Lemma 2: Let $p_{ij}$ be the probability of $P(t)$ in $\mathcal{S}_{i}(i=1,2)$ while $P(t+1)$ in $\mathcal{S}_{j}(j=1,2)$. Then we have:

For $P(t)\in\mathcal{S}_{1}$, we have $p_{11}=1.$ 2. 2.

For $P(t)\in\mathcal{S}_{2}$, there exists $c\in(0,1)$ which satisfies $p_{22}\leq c$.

Proof:

In the proposed algorithm, we adopt elitist selection as the choice method for populations, if $P(k)\in\mathcal{S}_{1}$, then $P(k)\in\mathcal{S}_{1}$, i.e., $p_{11}=1$. Therefore, 1) is proved.

To prove 2), we assume that $x^{*}$ is the optimal solution. Then $\exists\gamma>0$, there exists

[TABLE]

for $\forall x\in D\bigcap N_{\gamma}(x^{*})$, where $N_{\gamma}(x^{*})=\{x\in S:|x-x^{*}|\leq\gamma\}$. Based on (58), we can obtain

[TABLE]

In the case of $P(t)\in\mathcal{S}_{2}$, we assume $\mathbf{x}_{mu}$ is the results of the proposed algorithm after step 4 with population $P(k)$. According to (45) and (46), the probability of $\mathbf{x}_{mu}=\mathbf{x}+\Delta\mathbf{m}$ can be expressed as

[TABLE]

Based on (60), the probability of $\mathbf{x}_{mu}\in N_{\gamma}(x^{*})\bigcap D$ can be calculated as

[TABLE]

where $x_{i}$ and $x_{i}^{*}$ denote the $i$-th element of $\mathbf{x}$ and $\mathbf{x}^{*}$, respectively. For clarity, we use $p^{*}(\mathbf{x})$ to denote $\mathbb{P}(\mathbf{x}+\Delta\mathbf{m}\in N_{\gamma}(x^{*})\cap D)$. As $N_{\gamma}(x^{*})\cap D$ is a bounded, non-void, closed set, and based on (63), we can obtain

[TABLE]

Moreover, as $D$ is a bounded closed set, we can obtain

[TABLE]

Based on (59) and (63), we can get

[TABLE]

Supposing $c=1-p^{*}(\hat{\mathbf{x}})$, we can obtain

[TABLE]

The Lemma 2.2) is proved. ∎

Based on Lemma 2, we denote $p_{t}=\mathbb{P}(|f(\mathbf{x}_{t}^{*})-f(\mathbf{x}^{*})|\geq\delta)$, where $\mathbf{x}_{t}^{*}=\min_{P(t)\cap D}f(\mathbf{x})$, then we can get

[TABLE]

According to Lemma 2, $\hat{p}_{t}$ can be estimated by $c$ as

[TABLE]

Therefore, we can obtain

[TABLE]

According to Lemma 1, we can obtain

[TABLE]

Finally, based on Definition 3, it is proved that the proposed algorithm achieves global convergence.

The proof is completed.