# A Hybrid Position/Force Controller for Joint Robots

Shengwen Xie1 and Juan Ren1 *This work has been published in the proceedings
of ICRA 2021 and we have added some extra explanations in this version.1
Department of Mechanical Engineering, Iowa State University, Ames, IA 50011,
USA swxie<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a hybrid position/force controller for operating
joint robots. The hybrid controller has two goals—motion tracking and force
regulating. As long as these two goals are not mutually exclusive, they can be
decoupled in some way. In this work, we make use of the smooth and invertible
mapping from the joint space to the task space to decouple the two control
goals and design controllers separately. The traditional motion controller in
the task space is used for motion control, while the force controller is
designed through manipulating the desired trajectory to regulate the force
indirectly. Two case studies—contour tracking/polishing surfaces and grabbing
boxes with two robotic arms—are presented to show the efficacy of the hybrid
controller, and simulations with physics engines are carried out to validate
the efficacy of the proposed method.

##  I Introduction

In many automation tasks in which the robot interacts with the environment,
both position and force of the end-effector need to be controlled to complete
the task, i.e., a hybrid motion/force controller is needed. For example, when
using a robotic arm to polish the car apron [1], both the position and contact
force should be regulated.

Controllers for joint robots with one control goal, e.g., force regulation or
position control, have been studied extensively. Position controller can be
designed in the joint space or task space [2, 3]. Designing force controller
is more challenging due to the interaction with the environment and the
interaction model is hard to determine. Impedance controllers and admittance
controllers can establish a mass-spring-damper model between the end-effector
and the environment, but force cannot be accurately controlled[4, 5]. By
incorporating force feedback in the position controller, accurate force
tracking can be achieved [6, 7]. When force sensor is not available,
sensorless approaches were developed [8]. Although various position or force
controllers are available, the choices for hybrid position/force controller
are still limited, and designing robust hybrid controller remains an active
research area.

One popular position/force controller is proposed in [9], in which the motion
and force controllers are applied in two complementary subspaces thus the
stability analysis of each controller can be inherited. However, the
controller is not robust since the constraint is hard to guarantee in the
entire process (e.g., contact loss can happen), and the orthogonality relies
on the choice of coordinate and metric [10, 2]. Marin et al. proposed to use
kinestatic filtering to overcome the limitation, but it still relies on the
constraints to hold, which makes the controller less robust for applications
without natural constraints[11]. Recently, the concept of energy tank based on
passivity has been used in robot controller design [12, 1, 13]. For example, a
Cartesian force/impedance control is developed in [1]. However, the controller
is not a position/force controller in the sense that the path to the setpoint
cannot be regulated. Iterative learning control (ILC) has been proposed to
increase the controller robustness when the nonlinear terms cannot be
accurately compensated, but ILC achieves the control results through
repetitions[14]. To deal with the lag of robot position control, Liu et al.
proposed to control the position and force with different control frequencies
[15]. Although the robustness has been improved, the bandwidth of the position
controller is limited. Therefore, robust hybrid position/force controllers
which can be used without natural constraints are still not available.

In this work, we present a hybrid position/force controller which does not
rely on the constraints and can regulate both the position and force applied
by the end-effector. It is assumed that the two control objectives do not
conflict, that is, achieving one objective does not prohibit achieving
another. Generally, there are 6 degree of freedoms (DOFs) in the task space.
If a smooth and invertible mapping exists such that the interaction force
depends on $n$ DOFs of the 6 DOFs, one can design force controller for the $n$
DOFs and motion controller for the rest $6-n$ DOFs. The motion controller
which utilizes a smooth and invertible mapping to transform the configuration
from the joint space to the task space is used [3]. Then, a force controller
is incorporated through manipulating the desired position corresponding to the
$n$ DOFs. The closed-loop stability of the hybrid controller is analyzed with
the linear interaction model between the end-effector and the environment.

To verify the effectiveness of the hybrid controller, two cases studies are
presented. The first case is inspired by the task in [1] which can be
generalized to a contour tracking problem with a desired contact force.
Compared to the approach in [1], the proposed controller is advantageous in
following the pre-designed trajectory while regulating the force in the
process. In the second case study, the direction of the force to be controlled
can be varying, but it always corresponds to 1 DOF among the 6 DOFs of the
end-effector in the task space. The key is to modify the smooth and invertible
mapping from the joint space to the task space to ensure that the force
controller is related with 1 DOF and thus the proposed hybrid controller
framework can be adopted. Although in all the case studies, the force
controller always corresponds to 1 DOF, it can be easily extended to more
DOFs. Moreover, the proposed controller can be used for mounting a screwdriver
onto a screw as reported in [19].

##  II Hybrid Controller Design

Let $R\in SO(3)$ denote the orientation of the end-effector and we can map $R$
to $\Theta=[\Theta_{x},\Theta_{y},\Theta_{z}]\in so(3)$[2] (p. 83). Thus, the
end-effector pose can be represented with
$[\Theta_{x},\Theta_{y},\Theta_{z},x,y,z]$ where $[x,y,z]$ is the position.
The hybrid controller design in this work is based on the assumption that the
interaction force in some direction can only be affected by the position in
that direction. For example, the torque around $x$ axis can only be changed
through varying $\Theta_{x}$ and the interaction force along $x$ axis can be
changed by manipulating position in $x$ direction. Note that the interaction
between the end-effector and the environment is very complex, the assumption
here is to simplify the theoretical analysis and the hybrid controller may
still work even when they are violated due to its robustness.

###  II-A Motion Control in the Task Space

We consider a robotic arm with 6 DOFs. The equation of motion is

$M(\bm{\theta})\ddot{\bm{\theta}}+C\dot{\bm{\theta}}+{\bm{g}}({\bm{\theta}})+J^{T}\pazocal{F}_{\rm
ext}=\bm{\tau}\text{,}$ (1)

where $\bm{\theta}\in\mathbb{R}^{6\times 1}$ is the joint angle vector, $J$ is
the arm’s Jacobian matrix, $M({\bm{\theta}})$,
$C({\bm{\theta}},\dot{\bm{\theta}})$, and $\bm{g}(\bm{\theta})$ are the mass
matrix, Coriolis and centrifugal vector, gravity vector, respectively, and
$\pazocal{F}_{\rm ext}=[\bm{m}^{T}_{\rm ext},\bm{f}^{T}_{\rm
ext}]^{T}=[m_{x},m_{y},m_{z},f_{x},f_{y},f_{z}]^{T}$ is the wrench applied by
the end-effector to the environment. The control input to the robot is
$\bm{\tau}$, the torque applied to the joint motors.

The motion controller in the task space from [3] (pp. 195-197) is used. Define
a smooth and invertible mapping from the joint space to task space
$f:\bm{\theta}\in\mathbb{R}^{6\times 1}\to\bm{x}\in\mathbb{R}^{6\times 1}$. In
this work, the joint space is mapped to the task space defined with the
position and orientation of the end-effector, i.e.,
$\bm{x}=[\alpha,\beta,\gamma,x,y,z]^{T}$. Note that this mapping is not
globally but locally invertible thus works in a region, so multiple mappings
may be used for the entire joint space. The Jacobian $J_{f}$ with respect to
$f$ is defined as

$\dot{\bm{x}}=J_{f}\dot{\bm{\theta}}\text{,}\quad\quad J_{f}=\frac{\partial
f}{\partial{\bm{\theta}}}\text{.}$ (2)

Thus,

$\ddot{\bm{x}}=\dot{J}_{f}\dot{\bm{\theta}}+J_{f}\ddot{\bm{\theta}}\quad\Rightarrow\quad\ddot{\bm{\theta}}=J_{f}^{-1}(\ddot{\bm{x}}-\dot{J}_{f}\dot{\bm{\theta}})\text{.}$
(3)

The motion controller can be expressed as

$\bm{\tau}=C\dot{\bm{\theta}}+{\bm{g}}({\bm{\theta}})+MJ_{f}^{-1}(\ddot{\bm{x}}_{d}+K_{v}\dot{\bm{x}}_{e}+K_{p}{\bm{x}}_{e}+K_{i}\int_{0}^{t}{\bm{x}}_{e}ds-\dot{J}_{f}\dot{\bm{\theta}})\textbf{,}$
(4)

where ${\bm{x}}_{e}={\bm{x}}_{d}-{\bm{x}}$; ${\bm{x}}_{d}$ and ${\bm{x}}$ are
the desired and actual configuration, respectively; $K_{v}$, $K_{p}$, and
$K_{i}$ are the controller parameters and should be chosen properly to ensure
stability [3] (pp. 198).

To introduce the hybrid controller, consider a task of cleaning a table with a
robotic arm. The table surface is a plane parallel to the $x\text{-}y$ plane
in world frame with the $z$ axis pointing upward. Therefore, for such a task,
a hybrid controller is needed to manipulate the end-effector to move from
point $A$ to point $B$ on the table while maintaining a downward force (i.e.,
in -$z$ direction).

Define the mapping
$f:{\bm{\theta}}\to{\bm{x}}=[\alpha(t),\beta(t),\gamma(t),x(t),y(t),z(t)]^{T}$,
where $x$, $y$, and $z$ are coordinates in world frame and the rest DOFs are
orientation angles. We can decouple the 6 DOFs in the task space, i.e., the
motion controller will be used to regulate
$\\{x(t),y(t),\alpha(t),\beta(t),\gamma(t)\\}$ to track the reference
trajectory and the force controller regulates the contact force in $z$
direction through manipulating $z(t)$. Let
${\bm{x}}_{d}=[\alpha_{d}(t),\beta_{d}(t),\gamma_{d}(t),x_{d}(t),y_{d}(t),z_{d0}]^{T}$
be the reference trajectory in the task space with $z_{r0}$ the initial $z$
position of the tip.

Plug Eqs.(4) and (3) into Eq. (1), we can obtain

$\ddot{\bm{x}}_{e}+K_{v}\dot{\bm{x}}_{e}+K_{p}{\bm{x}}_{e}+K_{i}\int_{0}^{t}{\bm{x}}_{e}ds-
J_{f}M^{-1}J^{T}\pazocal{F}_{\rm ext}=0\text{.}$ (5)

Let $[\cdot]_{i}$ denote the $i$th elements of a vector and
$u=J_{f}M^{-1}J^{T}\pazocal{F}_{\rm ext}$. Eq. (5) can be decomposed into two
parts as follows.

${[\ddot{\bm{x}}_{e}]}_{i}+K_{v}{[\dot{\bm{x}}_{e}]}_{i}+K_{p}{[{\bm{x}}_{e}]}_{i}+K_{i}\int_{0}^{t}{[{\bm{x}}_{e}]}_{i}ds-[u]_{i}=0\text{,}1\leq
i\leq 5$ (6)

$[\ddot{\bm{x}}_{e}]_{6}+K_{v}{[\dot{\bm{x}}_{e}]}_{6}+K_{p}{[{\bm{x}}_{e}]}_{6}+K_{i}\int_{0}^{t}{[{\bm{x}}_{e}]}_{6}ds-[u]_{6}=0$
(7)

Eq. (6) represents the motion error dynamics while Eq. (7) will be used for
the force controller design. Note that $\pazocal{F}_{\rm ext}$ is not canceled
in Eq. (5) since it is difficult to accurately measure the term
$\pazocal{F}_{\rm ext}$ and there will be a delay in canceling the term
resulting in some unexpected behaviors.

Remark 1. The force term $\pazocal{F}_{\rm ext}$ can be canceled (then $u=0$)
if accurate sensors are available and the conclusion will not be affected. We
are going to discuss the effect if this term is not canceled or just partially
canceled.

###  II-B Hybrid Control: Part 1-Motion Control

Eq. (6) is the motion error dynamics after applying the motion controller Eq.
(4). Here we explain how the term $\pazocal{F}_{\rm ext}$ affects the
controller performance. If $\pazocal{F}_{\rm ext}$ is not canceled then
$[u]_{i}=[J_{f}M^{-1}J^{T}\pazocal{F}_{\rm ext}]_{i}$; if $\pazocal{F}_{\rm
ext}$ is not fully canceled, $[u]_{i}$ represents the residual term after
canceling, then we can obtain the transfer function of the system with input
$u_{i}(t)$ and output $y_{i}(t)=[X_{e}]_{i}$ from Eq. (6) as

$G_{i}(s)=\frac{Y_{i}(s)}{U_{i}(s)}=\frac{s}{s^{3}+K_{v}s^{2}+K_{p}s+K_{i}}\text{.}$
(8)

If $K_{v}$, $K_{p}$, and $K_{i}$ are chosen such that Eq. (8) is stable, it
yields a low pass filter which has a very small gain in low frequency region.
Thus, $u_{i}$ does not affect the tracking error $y_{i}$ much. In particular,
when $s\to 0$, $G_{i}(s)\to 0$ which implies that constant $u_{i}$ cannot
affect the output at all.

###  II-C Hybrid Control: Part 2-Force Control

Let ${z}_{e}={z}_{d}-{z}_{t}={[{\bm{x}}_{e}]}_{6}$ where ${z}_{d}$ is the
desired $z$ position and $z_{t}$ is the actual $z$ position of the end-
effector. Let $f_{e}=f_{d}-f_{z}$, with $f_{d}$ the desired downward force and
$f_{z}$ the actual downward force. We integrate another controller Eq. (10)
with Eq. (7) as follows (Note that (9) is obtained by replacing
${[\bm{x}_{e}]}_{6}$ with ${z}_{e}$ in (7)).

$\displaystyle\ddot{z}_{e}+K_{v}\dot{z}_{e}+K_{p}{z}_{e}+K_{i}\int_{0}^{t}{z}_{e}ds-(\epsilon
f_{z}+\eta)=0$ (9) $\displaystyle
u_{c}=\dot{z}_{d}=K_{p1}f_{e}+K_{i1}\int_{0}^{t}f_{e}ds\text{~{},}$ (10)

where $(\epsilon f_{z}+\eta)=[J_{f}M^{-1}J^{T}\pazocal{F}_{\rm ext}]_{6}$,
$u_{c}$ is the input for the force controller. In Eq. (9), $f_{z}$ is
extracted and the rest terms in $[J_{f}M^{-1}J^{T}\pazocal{F}_{\rm ext}]_{6}$
are accounted for with $\epsilon$ and $\eta$. Again, if $\pazocal{F}_{\rm
ext}$ can be canceled, $\epsilon=0$, $\eta=0$. The reason why $f_{z}$ is
extracted here is that $f_{z}$ can be very large thus dominates over other
terms in $\pazocal{F}_{\rm ext}$ (i.e., $f_{x}$, $m_{x}$, etc.) when $f_{d}$
is large.

Controller Eq. (10) essentially indirectly regulates the contact force through
adjusting the desired $z$ position. The block diagram of the force controller
is shown in Fig. 1. We use a spring model to approximate the interaction
dynamics between the tip and the environment. As seen in Fig. 2, let $q_{z}$
be the $z$ position of the contact point before contacting (before contacting
there is no deformation at the contact point), the contact force can be
computed with $f_{z}=\delta(k)({q_{z}}-z_{t})$, where

$\begin{aligned} \delta(k)=\begin{cases}&k\text{~{}, if}~{}~{}q_{z}>z_{t}\\\
&0\text{~{}, if}~{}~{}q_{z}\leq z_{t}\end{cases}\end{aligned}\text{.}$ (11)

Sometimes the environment may be rigid and there is no deformation, the model
can still work due to the flexibility of the robotic arm. When implementing
the controller, $f_{z}$ can be read from the F/T sensor.

Figure 1: Control diagram for the force controller. Figure 2: Interaction
model between the end-effector and the environment.

To obtain the closed-loop force error dynamics, the system state
$W=[w_{1},\cdots,w_{6}]^{T}$ is defined as follows. For simplicity, it is
assumed that $\dot{f}_{d}=0$, but this assumption can be removed and similar
analysis can be carried out.

$\begin{aligned} \begin{cases}w_{1}&=\int_{0}^{t}z_{e}ds\\\ w_{2}&=z_{t}\\\
w_{3}&=\dot{z}_{e}\\\ w_{4}&=z_{d}\\\ w_{5}&=\int_{0}^{t}f_{e}ds\\\
w_{6}&=f_{e}\end{cases}\end{aligned}\Rightarrow\begin{aligned}
\begin{cases}\dot{w}_{1}&=z_{e}=z_{d}-z_{t}=w_{4}-w_{2}\\\
\dot{w}_{2}&=\dot{z}_{t}=\dot{z}_{d}-\dot{z}_{e}=u_{c}-w_{3}\\\
&=-w_{3}+K_{i1}w_{5}+K_{p1}w_{6}\\\
\dot{w}_{3}&=\ddot{z}_{e}=-K_{v}w_{3}-K_{p}(w_{4}-w_{2})\\\
&~{}~{}-K_{i}w_{1}+(\epsilon f_{z}+\eta)~{}(\text{by Eq.
(\ref{controller2_1})})\\\
\dot{w}_{4}&=\dot{z}_{d}=u_{c}=K_{i1}w_{5}+K_{p1}w_{6}\\\
\dot{w}_{5}&=f_{e}=w_{6}\\\
\dot{w}_{6}&=\dot{f}_{e}=-\dot{f}_{z}\end{cases}\end{aligned}$ (12)

The system state-space equation is

$\begin{aligned} \dot{W}&=A_{w}W+BU_{w}\\\ =&\begin{bmatrix}0&-1&0&1&0&0\\\
0&0&-1&0&K_{i1}&K_{p1}\\\ -K_{i}&\epsilon\delta(k)+K_{p}&-K_{d}&-K_{p}&0&0\\\
0&0&0&0&K_{i1}&K_{p1}\\\ 0&0&0&0&0&1\\\
0&0&-\delta(k)&0&\delta(k)K_{i1}&\delta(k)K_{p1}\end{bmatrix}W\\\
+&\begin{bmatrix}0&0\\\ 0&0\\\ 1&0\\\ 0&0\\\ 0&0\\\
0&1\end{bmatrix}\begin{bmatrix}\eta-\epsilon\delta(k)q_{z}\\\
-\delta(k)\dot{q}_{z}\end{bmatrix}\end{aligned}\text{.}$ (13)

Next, we analyze the case when $\epsilon(t)=0$ and $\eta(t)=0$, i.e.,
$\pazocal{F}_{\rm ext}$ can be canceled. If $\epsilon(t)\neq 0$, the system
Eq. (13) will be a time-varying model which is more challenging to analyze its
closed-loop stability.

Non-contact Mode: If $q_{z}\leq z_{t}$ for $t\geq t_{0}$, there is no contact,
thus $w_{6}=0$, $f_{e}=f_{d}$, the controller Eq. (10) becomes

$\displaystyle u_{c}(t)$
$\displaystyle=K_{p1}f_{d}+K_{i1}\int_{0}^{t_{0}}f_{e}ds+K_{i1}f_{d}(t-t_{0})$
(14) $\displaystyle=\dot{z}_{d}(t_{0})+K_{i1}f_{d}\cdot(t-t_{0})\text{,}$

where $t\geq t_{0}$. Eq. (13) degenerates to

$\dot{[W]}_{1:3}=\begin{bmatrix}0&-1&0\\\ 0&0&-1\\\
-K_{i}&+K_{p}&-K_{d}\end{bmatrix}[W]_{1:3}+\begin{bmatrix}z_{d}(t_{0})+\int_{t_{0}}^{t}u_{c}(t)ds\\\
u_{c}(t)\\\ 0\end{bmatrix}\text{.}$ (15)

As seen in Eq. (14), $u_{c}(t)$ will monotonously decrease or increase, which
may result in very large $u_{c}(t)$ before reaching the surface again if
$\|z_{d}-q_{z}\|$ is very large. To avoid this and make the controller more
robust, we can add a saturation function to the controller. Specifically,
assume the estimated contact point position $\hat{q}_{z}$ satisfies
$-\Delta<\hat{q}_{z}-q_{z}<0$ ($\Delta>0$ can be seen as the estimation error
bound), the controller Eq. (10) is modified to

$\displaystyle
u_{c}=\begin{cases}&S(K_{p1}f_{e}+K_{i1}\int_{0}^{t}f_{e}ds)\text{~{}~{}if}~{}~{}f_{z}=0\\\
&K_{p1}f_{e}+K_{i1}\int_{0}^{t}f_{e}ds\text{~{}~{}~{}~{}~{}~{}if}~{}~{}f_{z}\neq
0\end{cases}$ (16)

with

$\displaystyle S(z)=\begin{cases}&0\text{~{}~{}if}~{}~{}z_{d}<\hat{q}_{z}\\\
&z\text{~{}~{}otherwise}\end{cases}\text{.}$ (17)

Since $K_{p}$, $K_{v}$ and $K_{i}$ can be chosen to make Eq. (15) stable, the
end-effector will eventually reach the surface again, then the controller will
switch to contact mode as explained next.

Contact Mode: If $q_{z}>z_{t}$, choose $K_{i1}$ and $K_{p1}$ such that all the
non-zero eigenvalues of $A_{w}$ have negative real parts. However, since
$\text{rank}(A_{w})=5$, there will be a zero eigenvalue for $A_{w}$. We have
the following proposition.

Proposition 1. In contact mode, suppose $\pazocal{F}_{\rm ext}$ can be
canceled, i.e., $\epsilon=0$ and $\eta=0$, choose $K_{i1}$ and $K_{p1}$ such
that all the non-zero eigenvalues have negative real parts, then,

  1. 1.

if $q_{z}=0$ and $\dot{q}_{z}=0$, $f_{e}$ converges to 0, i.e., for the system
$\dot{W}=A_{w}W$,
$\lim\limits_{t\to\infty}f_{e}=\lim\limits_{t\to\infty}w_{6}=0$;

  2. 2.

if $q_{z}\neq 0$ and $\dot{q}_{z}\neq 0$,
$\lim\limits_{t\to\infty}\|f_{e}\|=\lim\limits_{t\to\infty}\|w_{6}\|\leq\Gamma\|U_{w}\|_{\infty}$,
where $\Gamma$ is a constant.

Proof. 1) Here we consider the case in which $A_{w}$ has distinct eigenvalues,
if $A_{w}$ has repeated eigenvalues, as long as they have negative real parts,
the similar analysis can be performed based on linear system theory without
changing the conclusion [16].

The solution of $\dot{W}=A_{w}W$ is
$W(t)=\sum\limits_{i}c_{i}e^{\lambda_{i}t}{\bm{v}}_{i}$ ($c_{i}$ is a constant
scalar depending on the initial state, ${\bm{v}}_{i}$ is the eigenvector
corresponding to the eigenvalue $\lambda_{i}$),
$\lim_{t\to\infty}c_{i}e^{\lambda_{i}t}{\bm{v}}_{i}=0$ if
$\text{real}(\lambda_{i})<0$. Since there is only one zero eigenvalue
$\lambda_{6}$, thus $\lim_{t\to\infty}W(t)=c_{6}{\bm{v}}_{6}$. The eigenvector
corresponding to the zero eigenvalue is ${\bm{v}}_{6}$, then
$A_{w}{\bm{v}}_{6}=0$, we have

$\begin{aligned}
\begin{cases}-[{\bm{v}}_{6}]_{3}+K_{i1}[{\bm{v}}_{6}]_{5}+K_{p1}[{\bm{v}}_{6}]_{6}=0\\\
K_{i1}[{\bm{v}}_{6}]_{5}+K_{p1}[{\bm{v}}_{6}]_{6}=0\\\
[{\bm{v}}_{6}]_{6}=0\end{cases}\end{aligned}\Rightarrow\begin{aligned}
\begin{cases}[{\bm{v}}_{6}]_{3}=0\\\ [{\bm{v}}_{6}]_{5}=0\\\
[{\bm{v}}_{6}]_{6}=0\end{cases}\end{aligned}\text{.}$ (18)

It follows that
$\lim\limits_{t\to\infty}f_{e}=\lim\limits_{t\to\infty}w_{6}=c_{6}[{\bm{v}}_{6}]_{6}=0$,
implying that the system is asymptotically stable in contact mode. Eq. (18)
also indicates that $\dot{z}_{e}\to 0$ and $\int_{0}^{t}f_{e}ds\to 0$ as
$t\to\infty$. $\blacksquare$

2) Let $y=f_{e}=CW=[0,0,0,0,0,1]W$ be the output of the system. Then, the
state-space model can be expressed as

$\begin{aligned} \dot{W}&=A_{w}W+BU_{w}\\\ y(t)&=CW\end{aligned}\textbf{.}$
(19)

The solution of Eq. (19) is known as [16]

$y(t)=Ce^{A_{w}t}W(0)+C\int_{0}^{t}e^{A_{w}(t-\tau)}BU_{w}(\tau)d\tau\text{,}$
(20)

where $W(0)$ denotes the state at $t=0$. As $A_{w}$ has distinct eigenvalues,
$A_{w}$ can be decomposed as $A_{w}=VDV^{-1}$, where
$D=\text{diag}\\{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5},0\\}$
and real$(\lambda_{i})<0$, the $i$th column of $V$ is an eigenvector of
$A_{w}$ corresponding to eigenvalue $\lambda_{i}$. Therefore,
$e^{A_{w}t}=Ve^{Dt}V^{-1}$,
$e^{Dt}=\text{diag}\\{e^{\lambda_{1}t},\cdots,e^{\lambda_{5}t},1\\}$. Let
${\bm{\gamma}}_{6}$ denote the 6th row of $V$, then

$\displaystyle\lim\limits_{t\to\infty}Ce^{A_{w}t}W(0)$
$\displaystyle=\lim\limits_{t\to\infty}CV\text{diag}\\{e^{\lambda_{1}t},\cdots,e^{\lambda_{5}t},1\\}V^{-1}W(0)$
$\displaystyle={\bm{\gamma}}_{6}\text{diag}\\{0,0,0,0,0,1\\}V^{-1}W(0)=0\text{,}$

where we have used the fact that $[{\bm{\gamma}}_{6}]_{6}=0$ which is proven
by Eq. (18). Let $V^{-1}=[{\bm{\alpha}}_{1},\cdots,{\bm{\alpha}}_{6}]$, where
${\bm{\alpha}}_{i}$ is a column vector, we have

$\begin{aligned} &\int_{0}^{t}Ce^{A_{w}(t-\tau)}BU_{w}(\tau)d\tau\\\
&=\int_{0}^{t}CV\text{diag}\\{e^{\lambda_{1}(t-\tau)},\cdots,e^{\lambda_{5}(t-\tau)},1\\}V^{-1}BU_{w}(\tau)d\tau\\\
&=\int_{0}^{t}{\bm{\gamma}}_{6}\text{diag}\\{e^{\lambda_{1}(t-\tau)},\cdots,e^{\lambda_{5}(t-\tau)},1\\}[{\bm{\alpha}}_{3},{\bm{\alpha}}_{6}]U_{w}(\tau)d\tau\\\
&=\int_{0}^{t}[[{\bm{\gamma}}_{6}]_{1}e^{\lambda_{1}(t-\tau)},\cdots,[\bm{\gamma}_{6}]_{5}e^{\lambda_{5}(t-\tau)},0][{\bm{\alpha}}_{3},{\bm{\alpha}}_{6}]U_{w}(\tau)d\tau\end{aligned}\text{.}$

Since
$\lim\limits_{t\to\infty}\int_{0}^{t}e^{\lambda_{i}(t-\tau)}d\tau<\infty$ for
$\lambda_{i}<0$, there exists a constant $\epsilon_{i}>0$ such as
$\lim\limits_{t\to\infty}\|\int_{0}^{t}e^{\lambda_{i}(t-\tau)}d\tau\|<\epsilon_{i}$.
Therefore,

$\begin{aligned}
&\lim\limits_{t\to\infty}\|f_{e}\|=\lim\limits_{t\to\infty}\|y(t)\|\\\
&\leq\left\lVert\int_{0}^{t}[[\bm{\gamma}_{6}]_{1}e^{\lambda_{1}(t-\tau)},\cdots,[\bm{\gamma}_{6}]_{5}e^{\lambda_{5}(t-\tau)},0][\bm{\alpha}_{3},\bm{\alpha}_{6}]U_{w}(\tau)d\tau\right\rVert\\\
&\leq[|[\bm{\gamma}_{6}]_{1}|\epsilon_{1},\cdots,|[\bm{\gamma}_{6}]_{5}|\epsilon_{5},0][|\bm{\alpha}_{3}|,|\bm{\alpha}_{6}|][1,1]^{T}\cdot\|U_{w}\|_{\infty}\\\
&=\Gamma\|U_{w}\|_{\infty}\end{aligned}\text{,}$

where $\gamma$ is a constant, and $|\cdot|$ on a vector is element-wise, i.e.,
taking absolute value for each entry of a vector.

If $A_{w}$ has repeated eigenvalues, $D$ will be a Jordan matrix and similar
arguments can be used to prove 2) with the linear system theory [16].
$\blacksquare$

Remark 2. 2) also implies that the system is bounded-input-bounded-output
(BIBO)[17].

Hybrid Mode: Usually, the system will start from non-contact mode and remain
in contact mode, but it is also possible that it switches between the two
modes when $q_{z}$ is changing. In such case, it is difficult to show the
asymptotically stability, however, we can show that the system output is
bounded with the proposed controller applied: if the system stays in contact
mode for infinite time, the system is asymptotically stable thus bounded based
on the above discussion; if it switches to non-contact mode at $t_{0}$, then
the output at $t_{0}$ must be a finite value thus bounded, once entering non-
contact mode the output is still bounded due to the saturation function in the
modified controller (Eq. (16)). Therefore, the output $z_{t}$ is always
bounded with the proposed force controller.

Remark 3. Eq. (10) is a PI controller, it can also be upgraded to a PID
controller as
$u_{c}=\ddot{z}_{d}=K_{v1}\dot{f}_{e}+K_{p1}f_{e}+K_{i1}\int_{0}^{t}f_{e}ds$.
The analysis can be carried out in the same way by adding an extra state
$w_{7}=\dot{z}_{d}$ to Eq. (12). Note that this analysis is quite
conservative. In simulation (e.g., in case study A, the topography $q_{z}$ is
changing), it turns out that $\lim_{t\to\infty}f_{e}=0$ when $\dot{q}_{z}$ is
not large.

Remark 4. In Eq. (12), $A_{w}$ is rank-deficient due to the linear interaction
model. The rank deficiency is not guaranteed if other interaction models were
used.

Remark 5. In non-contact mode, extra information is needed to make sure that
end-effector can reach the object surface, i.e., entering contact-mode again.
In this case, this extra information is that by moving in -$z$ direction, it
is guaranteed to reach the surface. Note that the signs of $K_{p1}$ and
$K_{i1}$ are determined with such extra information.

##  III Case Studies

###  III-A Contour Tracking/Polishing Parts

This case is similar to the task presented in [1]. The goal is to explore
everywhere on the object surface (position control) with constant downward
force (force control). We can directly apply the hybrid controller here with
force in $z$ direction maintained. The process is also equivalent to contour
tracking with constant force as shown in Fig. 3. In the simulation the object
is an ellipsoid as shown in Fig. 3.

Figure 3: (a) Setup for Case A. (b) Initial state. (c) Intermediate state
during contour tracking.

###  III-B Grabbing a Box with Two Robotic Arms

Using two arms to grab a box is a common task for human beings. To mimic this
with robots, consider the scenario in which two robotic arms are used to grab
a box as shown in Fig. 4 (a).

During the process, the end-effector of two robots have to track the pre-
designed trajectories, and maintain a force in the $z$ direction of the frame
with origin at $C_{a}$ as shown in Fig. 4 (b). Therefore a motion/force
controller is needed. The proposed hybrid controller can be applied in this
case through modifying the mapping $f$ in Eq. (3). Note that the hybrid
controller will only be applied to robot B as shown in Fig. 4 (a) in order to
establish an invertible and smooth mapping from the joint space to the task
space as explained next. Robot A is controlled with a motion controller.

As seen in Fig. 4 (a), the frame ($C_{a}$ is the frame origin) attached to the
end-effector of robot A can be determined with
$\phi(t)=\\{(x_{a}(t),y_{a}(t),z_{a}(t)),[V_{x}(t),V_{y}(t),V_{z}(t)]\\}$,
where $(x_{a}(t),y_{a}(t),z_{a}(t))$ is the coordinate of $C_{a}$ in world
frame and $[V_{x}(t),V_{y}(t),V_{z}(t)]^{T}\in\mathbb{R}^{3\times 1}$
represents the directions of the $x$, $y$, and $z$ axes of the attached frame
in world frame. Note that once the trajectory of end-effector of robot A is
given, $\phi(t)$ is fixed at time $t$ assuming the motion controller works
effectively. Recall that the aforementioned mapping $f:{\bm{\theta}}\in
Q\to\bm{x}\in\mathbb{R}^{6\times 1}$ with
$\bm{x}=[\alpha,\beta,\gamma,x,y,z]^{T}$. Now for the end-effector of robot B,
define a new mapping $h:\bm{x}\to\bar{\bm{x}}$ where
$\bar{\bm{x}}=[\alpha,\beta,\gamma,\bar{x},\bar{y},\bar{z}]^{T}$ and
$\bar{x},\bar{y},\bar{z}$ can be computed as follows.

$\begin{bmatrix}\bar{x}\\\ \bar{y}\\\
\bar{z}\end{bmatrix}=\begin{bmatrix}V_{x}^{T}\\\ V_{y}^{T}\\\
V_{z}^{T}\end{bmatrix}\Bigg{(}\begin{bmatrix}x\\\ y\\\
z\end{bmatrix}-\begin{bmatrix}x_{a}\\\ y_{a}\\\ z_{a}\end{bmatrix}\Bigg{)}$
(21)

It can be verified that $h$ is a smooth and invertible mapping. Therefore,
$h\circ f$, the composition of two smooth and invertible mappings, is still a
smooth and invertible mapping. Essentially, $h$ maps the coordinate of $C_{b}$
to the frame originated at $C_{a}$. Thus the force in the $z$ direction of
frame $C_{a}$ is determined by $\bar{z}$ which is the same with the hybrid
controller we discussed before. Therefore, the same hybrid controller can be
applied by modifying the mapping from $f$ to $h\circ f$.

Note that this does not mean that the box is guaranteed to be picked up and
follow the desired trajectory, the force-closure conditions are required to
hold [2]. Here we assume that the force-closure conditions hold by adjusting
properly the initial orientation and position of the two end-effectors as
shown in Fig. 4 (a).

Figure 4: Case B: grabbing with two robotic arms. Figure 5: (a) and (c)
Desired positions during control for constant and varying $q_{z}$,
respectively. (b) and (d) Force errors for constant and varying $q_{z}$,
respectively. Figure 6: The force profile in $z$ direction for case A. Figure
7: Position control errors for case B.

##  IV Simulation Results

###  IV-A Force Controller Simulation Results

To simulate the force controller Eq. (9)-(10), parameters of controller Eq.
(6) were set as $K_{v}=35$, $K_{p}=405$, and $K_{i}=1500$. The desired force
was $f_{z}=5$N. The PI force controller parameters are $K_{p1}=-0.05$ and
$K_{i1}=-0.01$. We simulated the contact model for three different spring
constants: $k=300$N/m, $k=1500$N/m, and $k=4500$N/m. At the beginning, the tip
position was $z_{t}=0.05$m. The initial desired position was $z_{d}=-0.001$m.

If $\epsilon=0$ and $\eta=0$, it can be checked that all the nonzero
eigenvalues of $A_{w}$ have negative real parts, thus the above analysis can
be applied. Fig. 5 (a) and (b) show the performance of the force controller
for different $k$s when $q_{z}=0$. It can be seen that the system starts from
non-contact mode and then switches to contact model and remains there for the
rest time. The force can be accurately regulated for different $k$s and
converges in about 0.1s.

To test the controller performance for time-varying $q_{z}$, we simulated the
controller with $q_{z}(t)=0.01\sin(2\pi 0.2t)+0.01$m. The results are shown in
Fig. 5(c) and (d). This can happen when the tip is sliding on a curved
surface. It can be seen that the converged force error is very small for
different $k$ even when the contact point changes from 0.02m to 0m in about
5s.

The robotic arm used is UR5 which has the size of a human arm and the
simulation platform is MuJoCo [18]. For both two cases, the controller
parameters are $K_{v}=35$, $K_{p}=405$, $K_{i}=1500$, $K_{p1}=0.008$ and
$K_{i1}=0.002$. The video is available in the supplementary material.

###  IV-B Case A: Contour Tracking/Polishing Parts

In this case, the topography $q_{z}$ is changing as shown in Fig. 3. The force
curve during the polishing process is shown in Fig. 6 (b) and the downward
force is maintained to $-5$N.

###  IV-C Case B: Grabbing with Two Robotic Arms

In this case, the goal is to lift the box with mass of 1kg and move along the
vector $[0.27,0.1,0.27]$m within 5s in the world frame. Note that the hybrid
controller is only applied to robotic arm B while the other one is controlled
by a motion controller alone. The force is regulated to 20N with the friction
coefficient set to be 0.9 and the force control error is shown in Fig. 8 (a).
The motion control errors for the two robots are shown in Fig. 8 (b) and (c),
respectively. Note that for robotic arm B the motion error are expressed in
term of $\bar{x}$ and $\bar{y}$ based on the mapping in Eq. (21). Again, the
hybrid controller works as expected with both the motion and force regulated.

Note that for both cases the force error is very noisy as can be seen in Figs.
6 and 8 (a). This is partly due to the fact that the environment in simulation
is very hard. The contact is shortly broken and then established again due to
the controller resulting in spikes in the plots. On the other hand, the
algorithm for calculating the interaction force in MuJoCo may also contribute
to this. Computing the interaction force is transformed to solving an
optimization problem in MuJoCo making it very challenging to set parameters to
change the elasticity of the environment.

Figure 8: (a) Force control errors for case B. Motion errors for (b) robotic
arm A and (c) robotic arm B.

##  V Conclusion and Future Work

In this paper, a hybrid position/force controller is proposed. The position
and force controllers correspond to different DOFs in the task space with the
smooth and invertible mapping. The simulation results validated the proposed
approach.

The interaction model used in this work is linear, which may not be the real
case, more complex interaction models will be considered in the future. More
examples in which the force controller corresponds to more DOFs will be
studied. In addition, experiments with real robots will be conducted to verify
the proposed method.

Although we have not observed Zeno behavior of the hybrid controller in
simulations, current analysis cannot rule out the possibility of this
phenomena which needs further investigations.

## ACKNOWLEDGMENT

Shengwen would like to thank Professor Yan-Bin Jia for the technical
discussions which improved the quality of this work.

## References

  * [1] C. Schindlbeck and S. Haddadin, “Unified passivity-based cartesian force/impedance control for rigid and flexible joint robots via task-energy tanks,” in 2015 IEEE International Conference on Robotics and Automation, pp. 440–447, IEEE, 2015.
  * [2] K. M. Lynch and F. C. Park, Modern Robotics.  Cambridge University Press, 2017.
  * [3] R. M. Murray, Z. Li, S. S. Sastry, and S. S. Sastry, A mathematical introduction to robotic manipulation.  CRC press, 1994.
  * [4] N. Hogan, “Impedance control: An approach to manipulation: Part i—theory,” 1985\.
  * [5] C. Ott, R. Mukherjee, and Y. Nakamura, “Unified impedance and admittance control,” in 2010 IEEE International Conference on Robotics and Automation, pp. 554–561, IEEE, 2010.
  * [6] F. Lange, W. Bertleff, and M. Suppa, “Force and trajectory control of industrial robots in stiff contact,” in 2013 IEEE International Conference on Robotics and Automation, pp. 2927–2934, IEEE, 2013.
  * [7] T. Kroger, B. Finkemeyer, M. Heuck, and F. M. Wahl, “Adaptive implicit hybrid force/pose control of industrial manipulators: Compliant motion experiments,” in 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems(IEEE Cat. No. 04CH37566), vol. 1, pp. 816–821, IEEE, 2004\.
  * [8] S. A. M. Dehghan, M. Danesh, and F. Sheikholeslam, “Adaptive hybrid force/position control of robot manipulators using an adaptive force estimator in the presence of parametric uncertainty,” Advanced Robotics, vol. 29, no. 4, pp. 209–223, 2015.
  * [9] M. H. Raibert and J. J. Craig, “Hybrid position/force control of manipulators,” 1981.
  * [10] J. Duffy, “The fallacy of modern hybrid control theory that is based on “orthogonal complements” of twist and wrench spaces,” Journal of Robotic systems, vol. 7, no. 2, pp. 139–144, 1990.
  * [11] A. G. Marin and R. Weitschat, “Unified impedance and hybrid force-position controller with kinestatic filtering,” in 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 3353–3359, IEEE, 2016.
  * [12] V. Duindam and S. Stramigioli, “Port-based asymptotic curve tracking for mechanical systems,” European Journal of Control, vol. 10, no. 5, pp. 411–420, 2004.
  * [13] C. T. Landi, F. Ferraguti, C. Fantuzzi, and C. Secchi, “A passivity-based strategy for coaching in human-robot interaction,” in 2018 IEEE International Conference on Robotics and Automation, pp. 1–6, IEEE, 2018.
  * [14] S. Chen, Z. Wang, A. Chakraborty, M. Klecka, G. Saunders, and J. Wen, “Robotic deep rolling with iterative learning motion and force control,” IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 5581–5588, 2020.
  * [15] G. Liu, L. Fang, B. Han, and H. Zhang, “Frequency-division based hybrid force/position control of robotic arms manipulating in uncertain environments,” Industrial Robot: the international journal of robotics research and application, 2020.
  * [16] C.-T. Chen, Introduction to linear system theory.  Holt, Rinehart and Winston, 1970.
  * [17] H. K. Khalil and J. W. Grizzle, Nonlinear systems, vol. 3.  Prentice hall Upper Saddle River, NJ, 2002.
  * [18] E. Todorov, “Convex and analytically-invertible dynamics with contacts and constraints: Theory and implementation in mujoco,” in 2014 IEEE International Conference on Robotics and Automation, pp. 6054–6061, IEEE, 2014\.
  * [19] S. Xie, “Mounting a screwdriver onto a screw using hybrid control,” Master’s thesis, Iowa State University, USA, 2020.


5G-NR

    5G New Radio
3GPP

    3rd Generation Partnership Project
ABS

    aerial base station
AC

    address coding
ACF

    autocorrelation function
ACR

    autocorrelation receiver
ADC

    analog-to-digital converter
AIC

    Analog-to-Information Converter
AIC

    Akaike information criterion
ARIC

    asymmetric restricted isometry constant
ARIP

    asymmetric restricted isometry property
ARQ

    Automatic Repeat Request
AUB

    asymptotic union bound
AWGN

    Additive White Gaussian Noise
AWGN

    additive white Gaussian noise
PSK

    asymmetric PSK
AWRICs

    asymmetric weak restricted isometry constants
AWRIP

    asymmetric weak restricted isometry property
BCH

    Bose, Chaudhuri, and Hocquenghem
BCHSC

    BCH based source coding
BEP

    bit error probability
BFC

    block fading channel
BG

    Bernoulli-Gaussian
BGG

    Bernoulli-Generalized Gaussian
BPAM

    binary pulse amplitude modulation
BPDN

    Basis Pursuit Denoising
BPPM

    binary pulse position modulation
BPSK

    Binary Phase Shift Keying
BPZF

    bandpass zonal filter
BSC

    binary symmetric channels
BU

    Bernoulli-uniform
BER

    bit error rate
BS

    base station
BW

    BandWidth
BLLL

     binary log-linear learning
CP

    Cyclic Prefix
CDF

    cumulative distribution function
CDF

    Cumulative Distribution Function
CDF

    cumulative distribution function
CCDF

    complementary cumulative distribution function
CCDF

    complementary CDF
CCDF

    complementary cumulative distribution function
CD

    cooperative diversity
CDMA

    Code Division Multiple Access
ch.f.

    characteristic function
CIR

    channel impulse response
CoSaMP

    compressive sampling matching pursuit
CR

    cognitive radio
CS

    compressed sensing
CS

    Compressed sensing
CS

    compressed sensing
CSI

    channel state information
CCSDS

    consultative committee for space data systems
CC

    convolutional coding
COVID-19

    Coronavirus disease
DAA

    detect and avoid
DAB

    digital audio broadcasting
DCT

    discrete cosine transform
DFT

    discrete Fourier transform
DR

    distortion-rate
DS

    direct sequence
DS-SS

    direct-sequence spread-spectrum
DTR

    differential transmitted-reference
DVB-H

    digital video broadcasting – handheld
DVB-T

    digital video broadcasting – terrestrial
DL

    DownLink
DSSS

    Direct Sequence Spread Spectrum
DFT-s-OFDM

    Discrete Fourier Transform-spread-Orthogonal Frequency Division Multiplexing
DAS

    Distributed Antenna System
DNA

    DeoxyriboNucleic Acid
EC

    European Commission
EED

    exact eigenvalues distribution
EIRP

    Equivalent Isotropically Radiated Power
ELP

    equivalent low-pass
eMBB

    Enhanced Mobile Broadband
EMF

    ElectroMagnetic Field
EU

    European union
EI

    Exposure Index
eICIC

    enhanced Inter-Cell Interference Coordination
FC

    fusion center
FCC

    Federal Communications Commission
FEC

    forward error correction
FFT

    fast Fourier transform
FH

    frequency-hopping
FH-SS

    frequency-hopping spread-spectrum
FS

    Frame synchronization
FS

    frame synchronization
FDMA

    Frequency Division Multiple Access
GA

    Gaussian approximation
GF

    Galois field
GG

    Generalized-Gaussian
GIC

    generalized information criterion
GLRT

    generalized likelihood ratio test
GPS

    Global Positioning System
GMSK

    Gaussian Minimum Shift Keying
GSMA

    Global System for Mobile communications Association
GS

    ground station
GMG

     Grid-connected MicroGeneration
HAP

    high altitude platform
HetNet

    Heterogeneous network
IDR

    information distortion-rate
IFFT

    inverse fast Fourier transform
IHT

    iterative hard thresholding
i.i.d.

    independent, identically distributed
IoT

    Internet of Things
IR

    impulse radio
LRIC

    lower restricted isometry constant
LRICt

    lower restricted isometry constant threshold
ISI

    intersymbol interference
ITU

    International Telecommunication Union
ICNIRP

    International Commission on Non-Ionizing Radiation Protection
IEEE

    Institute of Electrical and Electronics Engineers
ICES

    IEEE international committee on electromagnetic safety
IEC

    International Electrotechnical Commission
IARC

    International Agency on Research on Cancer
IS-95

    Interim Standard 95
KPI

    Key Performance Indicator
LEO

    low earth orbit
LF

    likelihood function
LLF

    log-likelihood function
LLR

    log-likelihood ratio
LLRT

    log-likelihood ratio test
LoS

    Line-of-Sight
LRT

    likelihood ratio test
LWRIC

    lower weak restricted isometry constant
LWRICt

    LWRIC threshold
LPWAN

    Low Power Wide Area Network
LoRaWAN

    Low power long Range Wide Area Network
NLoS

    Non-Line-of-Sight
Li-Fi

    light-fidelity
LED

    light emitting diode
LABS

    LoS transmission with each ABS
NLABS

    NLoS transmission with each ABS
MB

    multiband
MC

    macro cell
MDS

    mixed distributed source
MF

    matched filter
m.g.f.

    moment generating function
MI

    mutual information
MIMO

    Multiple-Input Multiple-Output
MISO

    multiple-input single-output
MJSO

    maximum joint support cardinality
ML

    maximum likelihood
MMSE

    minimum mean-square error
MMV

    multiple measurement vectors
MOS

    model order selection
${M}$-PSK

     $M$-ary phase shift keying
${M}$-PSK

     $M$-ary asymmetric PSK
MP

     multi-period
MINLP

    mixed integer non-linear programming
$M$-QAM

     $M$-ary quadrature amplitude modulation
MRC

    maximal ratio combiner
MSO

    maximum sparsity order
M2M

    Machine-to-Machine
MUI

    multi-user interference
mMTC

    massive Machine Type Communications
mm-Wave

    millimeter-wave
MP

    mobile phone
MPE

    maximum permissible exposure
MAC

    media access control
NB

    narrowband
NBI

    narrowband interference
NLA

    nonlinear sparse approximation
NLOS

    Non-Line of Sight
NTIA

    National Telecommunications and Information Administration
NTP

    National Toxicology Program
NHS

    National Health Service
LOS

    Line of Sight
OC

    optimum combining
OC

    optimum combining
ODE

    operational distortion-energy
ODR

    operational distortion-rate
OFDM

    Orthogonal Frequency-Division Multiplexing
OMP

    orthogonal matching pursuit
OSMP

    orthogonal subspace matching pursuit
OQAM

    offset quadrature amplitude modulation
OQPSK

    offset QPSK
OFDMA

    Orthogonal Frequency-division Multiple Access
OPEX

    Operating Expenditures
OQPSK/PM

    OQPSK with phase modulation
PAM

    pulse amplitude modulation
PAR

    peak-to-average ratio
PDF

    probability density function
PDF

    probability density function
PDF

    probability distribution function
PDP

    power dispersion profile
PMF

    probability mass function
PMF

    probability mass function
PN

    pseudo-noise
PPM

    pulse position modulation
PRake

    Partial Rake
PSD

    power spectral density
PSEP

    pairwise synchronization error probability
PSK

    phase shift keying
PD

    power density
$8$-PSK

     $8$-phase shift keying
PPP

    Poisson point process
PCP

    Poisson cluster process
FSK

    Frequency Shift Keying
QAM

    Quadrature Amplitude Modulation
QPSK

    Quadrature Phase Shift Keying
OQPSK/PM

    OQPSK with phase modulator
RD

    raw data
RDL

    ”random data limit”
RIC

    restricted isometry constant
RICt

    restricted isometry constant threshold
RIP

    restricted isometry property
ROC

    receiver operating characteristic
RQ

    Raleigh quotient
RS

    Reed-Solomon
RSSC

    RS based source coding
r.v.

    random variable
R.V.

    random vector
RMS

    root mean square
RFR

    radiofrequency radiation
RIS

    Reconfigurable Intelligent Surface
RNA

    RiboNucleic Acid
RRM

    Radio Resource Management
RUE

    reference user equipments
RAT

    radio access technology
RB

    resource block
SA-Music

    subspace-augmented MUSIC with OSMP
SC

    small cell
SCBSES

    Source Compression Based Syndrome Encoding Scheme
SCM

    sample covariance matrix
SEP

    symbol error probability
SG

    sparse-land Gaussian model
SIMO

    single-input multiple-output
SINR

    signal-to-interference plus noise ratio
SIR

    signal-to-interference ratio
SISO

    Single-Input Single-Output
SMV

    single measurement vector
SNR

    signal-to-noise ratio
SP

    subspace pursuit
SS

    spread spectrum
SW

    sync word
SAR

    specific absorption rate
SSB

    synchronization signal block
SR

    shrink and realign
tUAV

    tethered Unmanned Aerial Vehicle
TBS

    terrestrial base station
uUAV

    untethered Unmanned Aerial Vehicle
PDF

    probability density functions
PL

    path-loss
TH

    time-hopping
ToA

    time-of-arrival
TR

    transmitted-reference
TW

    Tracy-Widom
TWDT

    TW Distribution Tail
TCM

    trellis coded modulation
TDD

    Time-Division Duplexing
TDMA

    Time Division Multiple Access
Tx

    average transmit
UAV

    Unmanned Aerial Vehicle
URIC

    upper restricted isometry constant
URICt

    upper restricted isometry constant threshold
UWB

    ultrawide band
UWB

    Ultrawide band
URLLC

    Ultra Reliable Low Latency Communications
UWRIC

    upper weak restricted isometry constant
UWRICt

    UWRIC threshold
UE

    User Equipment
UL

    UpLink
WiM

    weigh-in-motion
WLAN

    wireless local area network
WM

    Wishart matrix
WMAN

    wireless metropolitan area network
WPAN

    wireless personal area network
WRIC

    weak restricted isometry constant
WRICt

    weak restricted isometry constant thresholds
WRIP

    weak restricted isometry property
WSN

    wireless sensor network
WSS

    Wide-Sense Stationary
WHO

    World Health Organization
Wi-Fi

    Wireless Fidelity
SpaSoSEnc

    sparse source syndrome encoding
VLC

    Visible Light Communication
VPN

    Virtual Private Network
RF

    Radio Frequency
FSO

    Free Space Optics
IoST

    Internet of Space Things
GSM

    Global System for Mobile Communications
2G

    Second-generation cellular network
3G

    Third-generation cellular network
4G

    Fourth-generation cellular network
5G

    Fifth-generation cellular network
gNB

    next-generation Node-B Base Station
NR

    New Radio
UMTS

    Universal Mobile Telecommunications Service
LTE

    Long Term Evolution
QoS

    Quality of Service

#  Ultra-Dense LEO Satellite-based Communication Systems: A Novel Modeling
Technique

Ruibo Wang, Mustafa A. Kishk, and Mohamed-Slim Alouini Ruibo Wang, Mustafa A.
Kishk and Mohamed-Slim Alouini are with King Abdullah University of Science
and Technology (KAUST), CEMSE division, Thuwal 23955-6900, Saudi Arabia

###### Abstract

Low earth orbit (LEO) satellite plays an indispensable role in the earth
network because of its low latency, large capacity, and seamless global
coverage. For such an unprecedented extensive irregular system, stochastic
geometry (SG) is a suitable research method. The SG model can not only cope
with the increasing network scale but also accurately analyze and estimate the
network’s performance. Several standard satellite distribution models and
satellite-ground channel models are investigated in this paper. System-level
metrics such as coverage probability and their intermediates are introduced in
the non-technical description. Then, the influence of some factors on latency
and coverage probability is studied. Finally, this paper presents the possible
challenges and corresponding solutions for the SG-based LEO satellite system
analysis.

##  I Introduction

In the last few years, we have been witnessing a boom in the industry of LEO
satellite networks. Constellations such as Starlink, Telesat, and OneWeb are
planning to or have already deployed thousands of satellites [1]. With the
relatively good quality in terms of capacity and latency, LEO satellites
provide an affordable solution for seamless network coverage. However, a
strong and tractable mathematical framework is still needed to study and
analyze the performance of LEO satellite-based communication systems. In
particular, with a system composed of large numbers of LEO satellites deployed
at different altitudes from the earth, analysis will enable us to understand
the required numbers and altitudes of satellites to maximize system
performance. In addition, it will enable us to to learn how much earth
stations are needed in order to enable seamless operation of the LEO satellite
communication system.

Based on the above description, the structure of this paper is as follows.
First, we describe the advantages and disadvantages of low-orbit satellites in
more detail. Then we give the motivation for applying SG to LEO satellite
networks. Based on SG, this paper goes over the analysis framework of some
recent literature and gives a non-technical description of models of satellite
locations distribution, channel model, typical parameters, and main
performance metrics. In addition, a simulation setup is proposed to highlight
the potential gains of using SG in modeling and analyzing LEO satellite
communication systems with useful system-level insights concluded. Finally, we
predict future scenarios for the challenges of applying SG in satellite
networks.

Figure 1: System-level metric comparison of three systems and three future
scenarios.

###  I-A Advantages and Challenges of LEO Satellite System

In this section, the advantages and challenges of the LEO satellite system are
analyzed by comparing with the geosynchronous equatorial orbit (GEO) satellite
system and the terrestrial system [2]. The radar diagram in the lower-left
corner of Fig. 1 visually shows the qualitative comparison points.

LEO satellite system can provide seamless coverage of the earth and has
incomparable advantages in providing reliable coverage. Despite the enormous
coverage of a single GEO satellite, the system has blind areas at the poles.
The coverage of terrestrial systems is limited by terrain. In remote areas
that lack infrastructure, existing ground networks are challenged to provide
coverage because the cost of extending fiber optics to these areas is
expensive.

LEO satellite system has relatively larger capacity density and lower latency.
In the same region, LEO satellite systems provide much more capacity than GEO
satellite systems, but less than terrestrial systems. The capacity is
sufficient to support high rate communications in rural and remote areas, but
not enough to meet the needs of densely populated cities. The latency of LEO
satellites ranges from several milliseconds to tens of milliseconds, which can
meet the needs of most application scenarios. Application scenarios requiring
high real-time performance, such as the Internet of vehicles, are still highly
dependent on BSs. The long communication distance of the GEO satellite system
leads to a much longer delay than the others.

However, LEO satellite systems face significant challenges in terms of cost
and frequency coordination. A large LEO satellite system requires high
construction and maintenance costs. Fortunately, aerospace manufacturers have
brought down launch costs, and streamlined production is bringing down the
cost of producing satellites. Due to the low altitude of the terrestrial
system, a large number of base stations need to be deployed to ensure
coverage. While construction costs are higher, maintenance costs of the ground
stations are lower than satellites because they are located on the ground. The
smaller GEO satellite system certainly has a cost advantage. In terms of
frequency licensing, the International Telecommunication Union (ITU) gives
frequency priority to GEO satellites. Because of the first-claim, first-
possession principle, once these huge constellations are built, it is not easy
to allocate dedicated spectrum to the new constellations. There is also some
channel overlap between LEO satellite system and 5G ground base station
system, which further complicates the challenges of spectrum resource
allocation problems.

###  I-B Motivations of Applying SG in LEO Satellite Network

As a powerful mathematical tool, SG has been widely used in cellular networks.
It also has advantages in LEO satellite system modeling and analysis.

The existing satellite analysis methods are mainly to design the specific
constellation. But for such a large system, fine-grained modeling of each LEO
satellite would take a lot of effort. It is far more tractable to do SG-based
system level analysis than to define and study the behavior of each satellite.
SG modeling does not need to rely on orbital metrics and the shape of a
specific constellation. A universal and flexible random network is what the
explosive growth of dynamic networks requires.

For satellite systems, more satellites and lower altitudes mean larger
density. In this case, interference will be an important factor affecting
system performance. However, apart from SG, there is still a lack of practical
analytical tools to simulate the interference caused by satellites.

In the existing non-SG models, satellites are often deployed in regular
circular or hexagonal cells with the same coverage [3]. However, the
distribution of satellites depends on latitude, and the coverage of satellites
varies significantly in most cases. Therefore, the results of traditional
analysis methods may greatly differ from the actual situation. SG is suitable
for modeling and analyzing irregular topological networks [4]. In recent
studies, various types of point processes used in LEO satellite network allow
multiple methods to simulate the distribution.

Applying SG requires the assumption that satellites are independent and
therefore can only be used to study the lower bounds of system performance.
However, it has been proved that the lower bound of coverage probability and
average achievable rate in the SG model is as tight as the upper bound of the
regular mesh model [5]. For the LEO satellite system, the performance of the
SG model closely matches with that of the deterministic constellation in
coverage probability and average achievable rate [6].

##  II SG-based Analytical Framework

In this section, we go over the analytical framework based on SG in recent
works. Table I at the end of this section provides a brief summary of these
frameworks. Before introducing the model, we first introduce some definitions
in astronomy.

  * •

Latitude: Latitude is the line plane angle between the normal line on the
earth and the equatorial plane. It ranges from 0 degrees at the equator to 90
degrees at the poles.

  * •

Zenith Angle: Zenith angle is the angle between the incident signal and the
direction of the local zenith (The geographical vertical).

  * •

Horizon: The horizon is a two-dimensional plane tangent to the earth through
the perception point, which can be the location of the gateway or the user.
From the gateway’s perspective, any satellite below the horizon is beyond line
of sight since the earth acts as a blockage.

###  II-A Models of Satellite Locations Distribution

Before studying satellite distribution, it is assumed that satellites, as well
as users, satellite gateways (GW), and BSs are independent and uniformly
distributed in the existing literature. The following sections will discuss
the three-point process distribution of satellite networks.

For Poisson Point Process (PPP), sample points obtained by random sampling
obey uniform distribution, and the number of points in any given area is
Poisson distributed while their locations are uniformly distributed within
this area. As the most frequently used point process, PPP can fit the ground
network well in performance analysis such as coverage probability and
achievable data rate. Therefore, the positions of the LEO satellite are
modeled as a PPP on a spherical surfaces of fixed height [7]. One of the most
significant advantages of PPP is that it can predict the probability of a
certain number of points in a spherical region with a specific area measure.

Although PPP has a strong practicality, it is not suitable for modeling a
finite area network with limited nodes. Since the positions of satellites are
on a spherical surfaces, and for deterministic constellations, the number of
satellites is fixed. Therefore, modeling the position of the satellite as a
Binomial Point Process (BPP) [6], [8] is an effective solution. Specifically,
it places a fixed number of points on the fixed height sphere. The zenith
(elevation) and azimuth angle are uniformly distributed.

Non-Homogeneous Poisson Point Process (NPPP) improves the PPP model from
another perspective. Since the distribution of deterministic satellite
constellations is not uniform at different latitudes, that is, the number of
satellites at the limit of inclination is actually greater than that in the
equatorial region, the number of effective satellites per latitude in NPPP is
used to compensate for the uneven density [9].

Corresponding to the labels in the upper left section of Fig. 1, there are
three typical scenarios: GW-relayed system, terrestrial station and LEO
satellite hybrid network and satellite to user direct communication scenario.
The GW can provide directional transmission to the satellite. By deploying a
GW as a relay, satellite coverage can be significantly improved [10].
Satellite networks can also assist ground networks. In addition to further
increasing system capacity, it can also provide communication during
emergencies. Satellites can also communicate directly with users, which of
course requires special user equipment to enable direct communication. Since
the coverage radius of GW or BS is negligible compared to the surface area of
the earth, GWs and BSs can be regarded as points located in a large two-
dimensional plane, and PPP is suitable to model their locations.

###  II-B Channel Model

The linear expression of the received power expressed as
$\rho\,\eta\,G\,r^{-\alpha}$, where $\rho$ is used to describe the transmitted
power and antenna gain, $\eta$ and $G$ are for Large scale fading and Small
scale fading, and $\alpha$ is called the path-loss exponent.

There are two main models of signal power transmission. Take the downlink
transmission as an example. The first is Equivalent Isotropically Radiated
Power (EIRP) criteria [7], that is, a satellite has the same gain in all
directions. The second is a satellite that has a larger gain in a particular
direction, in which the GW or user can receive main-lobe power, which can be
realized by beam-forming technology [11]. Due to directional transmission, the
received signal from the interfering source tends to come from the side-lobe,
which can improve the signal-to-interference ratio (SIR) of the system. In
addition, in the current literature, satellites at the same altitude are
typically assumed to have equal and constant transmission power.

The path-loss exponent is a number between $-4$ and $-2$ usually. Consider
that satellite-to-ground link propagates in free space in most cases,
$\alpha=-2$ is satisfied for open areas. However, when the ground clutter tier
reflects and absorbs the power, the multi-path effect cannot be ignored,
resulting in $\alpha=-4$. An intuitive result is that when $\alpha=-2$, the
interference becomes the main factor limiting system performance, and when
$\alpha=-4$, the influence of noise will be far more than that of
interference.

Large scale fading is usually modeled as a log-normal shadowing to describe
additional gain. It is inversely proportional to the square of carrier
frequency [6], [7], [8]. Furthermore, some literature take the attenuation of
air absorption caused by the resonance of water vapor into consideration [11],
[8], which is called rain attenuation. It changes with the geographical
location, such as desert and rain forest.

Small scale fading is a random variable which denotes the channel fading power
gains. In ideal free space propagation, Small scale fading does not occur.
Non-fading helps to investigate the upper bound on system performance. On the
contrary, Rayleigh fading helps investigate the lower bound, which happens
when serious multi-path distortion occurs. Rician fading is suitable for the
general situation and is widely applied to make accurate performance
estimates. Shadowed-Rician (SR) fading is a Rician fading channel with
fluctuating LoS components. It is the most accurate channel fading model,
especially for space-to-ground links [8].

Handling Non-Line-of-Sight (NLoS) interference is one of the keys to the
accuracy of a channel model. Satellites below the horizon may not provide
reliable service, but the diffraction signal may still interfere with the
user. At the same time, the satellite may not establish a line-of-sight link
due to obstacles from the ground. The former is often regarded as zero or some
small constant rather than a random variable because the signal is relatively
weak. The latter can be described by SR fading. It has been proved by
simulation that SR fading can accurately simulate the channel model under NLoS
interference [12]. In addition, multiplying by a factor representing the loss
of redundant paths in Large scale fading is also an effective solution. This
factor is modeled as a Log-Normal mixture random variable with two components,
and the mixing ratio is given by the probabilities of Line-of-Sight (LoS) and
NLoS [7]. Notice that the probability that an LoS link is established is only
determined by the zenith angle.

###  II-C Contact Distance and Satellite Availability

Contact distance is the distance between the user/GW and the associated
(service-providing) satellite. Based on the strongest average received power
association strategy, the contact distance is the distance from the user to
the nearest satellite at a given altitude. As a random variable, the
complementary cumulative distribution function (CCDF) of contact distance at
$d_{0}$ is equal to the probability that there are no satellites in the
spherical cap at the intersection of spherical surface over which satellites
are distributed and the cone with side length $d_{0}$ centered at the location
of the user/GW region closer to the user.

Contact angle is defined as the minimum zenith Angle between the user and the
satellite. As the distance between user and satellite increases, its zenith
Angle also increases, and there is a unique correspondence between the two.
Therefore, the contact Angle can be regarded as the representation of the
contact distance in spherical coordinates.

Figure 2: LEO satellite system structure diagram.

From spatial distribution perspectives, satellites, GWs, and users are all
located on spherical surfaces. We need to determine the range of achievable
contact distance. The upper left region of Fig. 2 shows the range of possible
locations of a serving LEO satellite for a given location of the user/GW in
the downlink. This range is mainly defined, as explained earlier, by the
horizon, below which no LoS satellites are available [11]. When Beamforming
(BF) is used, the main lobe beam-width is also one of the limitations [13].
The range of achievable contact distances is the intersection region above the
horizon and in the beam’s main-lobe width.

Nearest neighbor is used to describe the shortest distance from a selected
satellite to another. It is an important part of multi-hop network analysis.
The difference between nearest neighbor and contact distance is that, the
user/GW does not belong to the satellite point process, while the selected
satellite in the nearest neighbor does [10].

When the satellites are deployed at different altitudes with different values
of transmission power, the concept of association probability needs to be
introduced [14]. Especially in a multi-tier satellite network, it is
reasonable to adopt different transmitting power for different altitude
satellites. The satellite locations are divided into multiple point processes
according to their altitudes and transmission power. For each point process,
the distribution of the nearest satellite to the user needs to be obtained,
which is exactly equivalent to obtaining the contact distance. The association
probability with a given tier is the probability that the nearest satellite
from this tier provides higher average received power than the that of the
nearest satellites from all other tiers.

Satellite availability means that at least one satellite is within the visible
range of the user/GW. In small, low-altitude constellations, it is common to
have no usable satellites above the horizon. In the existing literature, there
are three methods to deal with the satellite availability problem: (i) reduce
the probability of having no usable satellites below an acceptable threshold
by increasing satellite density [10], (ii) enhance the probability of having
available satellites by increasing the height and transmitting power
simultaneously [11],[13], and (iii) calculate the probability of no available
LoS satellite and set SINR to zero if no LoS satellites are available [6]. The
first two proposals provide guidelines for constellation designers, and the
last proposal provides the most accurate estimation of system performance for
a designed constellation.

TABLE I: Models used in different literature.

References | Scenarios | Distribution of satellites | Small scale fading | System metrics |  System types
---|---|---|---|---|---
[6], [9], [13] | Direct communication |  Non-homogeneous PPP BPP |  Non-fading Rayleigh fading Rician fading |  Coverage probability Average achievable rate |  Noise-limited system Interference-limited system Generic system
[7], [11] | Hybrid network | PPP | Rayleigh fading | Coverage probability |  Ideal system Noise-limited system Generic system
[8], [10] |  Multi-tier system GW-relayed | BPP | Shadowed-Rician fading | Coverage probability | Noise-limited system

###  II-D Coverage Probability and Average Achievable Rate

In this section, we study the application of SG in system performance
analysis, with emphasis on coverage probability. The definitions of the three
metrics for measuring system performance are given below.

  * •

The  signal-to-interference plus noise ratio (SINR) is used to measure the
communication quality of the system.

  * •

Coverage probability is defined as the probability that the SINR is larger
than a predetermined acceptable threshold. It represents the probability that
the system can provide reliable connections.

  * •

Average achievable rate is the ergodic capacity from the Shannon–Hartley
theorem over a fading communication link.

After the distance to the associated satellite is determined by the contact
distance, the sum of the power of other LoS satellites is the interference.
SINR can be obtained then the coverage probability and average achievable rate
are calculated. When LEO satellite adopts multiple orthogonal bands, the
average achievable rate should be divided by the number of orthogonal bands.

Coverage probabilities in different application scenarios are expressed
differently. When LEO satellites rely on a ground GW as a relay, a user is
covered by a satellite when both satellite-GW and GW-user links achieve the
SINR requirements. Therefore, the system coverage probability should be the
product of the coverage probability on both links. In the terrestrial base
station and LEO satellite hybrid network scenario, satellites and BSs will
independently establish the network. The user is considered in outage when the
SINR of either the user-GW link or the GW-satellite link is below the
threshold.

Coverage probabilities in different system types is highlighted in the
literature. In an ideal System, the signal will not be affected by any
interference and noise in transmission. As long as there is a satellite in
service, the user can always acquire a successful transmission. Such systems
are used to study the upper bounds of performance and satellite availability.
A noise-limited system occurs in an environment where the noise power is much
stronger than the interference power. For example, when the frequency bands
used by satellites are orthogonal or the interference comes from a diffraction
signal, the interference power is often regarded as zero or a small constant.
In addition, although using orthogonal channels will significantly reduce the
Co-Channel Interference and improve the SINR and coverage probability, each
satellite can use only a portion of the system’s bandwidth, which reduces
capacity and data rates. In an interference-limited system, interference power
is dominant. Massive satellite constellations can cause huge interference. A
generic system takes both interference and noise power into account.

Coverage probabilities and average achievable rate under different types of
Small scale fading have completely different mathematical expressions.
Coverage probability and average achievable rate under Non-fading and Rayleigh
in the generic system are provided in [6]. In Rician fading and SR fading, the
associated satellite signal is transmitted in the direct path, while the
interfering satellite signals are transmitted through the reflection path. In
this case, the system is noise-limited. In [8], the coverage probability under
SR fading in a noise-limited system is obtained.

Finally, no matter what kind of Small scale fading is adopted, the Laplace
transform of interference is involved in the calculation of coverage
probability and achievable rate [6].

##  III SG-based Simulation Setup

Suppose two GWs at opposite ends of the earth, with zenith angles $0$ and
$\pi$ (not necessarily north and south poles), need to communicate. Satellites
are modeled as a BPP on a sphere concentric with the earth. The signals of the
inter-satellite link and the satellite-ground link travel at the speed of
light. As is shown in Fig. 2, GWs can send messages to the associated
satellites, while satellites can communicate with the GWs in the range of
line-of-sight. We designed two scenarios [15]: (i) in no inter-satellite links
scenario, the satellite has to rely on the GW as a relay. In other words, the
satellite uses GWs to forward data to other satellites until the message is
delivered to its destination, (ii) instead of requiring a GW relay, the
satellite sends a message directly to the next hop in inter-satellite links
scenario. Examples of both scenarios are shown in Fig. 2.

Figure 3: Average latency of satellite systems at different altitudes and
scales.

In Fig. 3, it can be seen that when the number of satellites is insufficient,
the time delay decreases with the increase of altitude. That’s because when a
satellite is looking for the next hop, it can only try to find satellites in
line of sight that might transmit signals. As it rises in altitude, the
satellite’s range of view expands, allowing it to choose which of several
satellites is in a better position. But for enough satellites, with the
increase of height, both the distances of inter-satellite link and satellite-
ground link will increase.

Overall, when S-S communication is not available, the delay increases by about
30 percent compared to the case when S-S communication is available. When the
number of satellites is increased from 100 to 300, the system delay can be
significantly reduced. From 300 to 1,000 satellites, performance no longer
improved much. Especially when S-S communication is not available, the curves
of 300 satellites and 1000 satellites overlap, and adding additional
satellites will not bring additional benefits in terms of delay. As the
altitude increases, the two considered scenarios tend to have different
latency values. The trends of the three curves when S-S communication is
available tend to be consistent faster. It suggests that increasing altitude,
increasing constellation size, and allowing interaction between satellites
lead to more satellites within a satellite’s reach. When there are enough
satellites in the range, the performance improvement of these methods is
limited, and the trend of the curves tends to be consistent.

In Fig. 4 and Fig. 5, we analyze the influence of satellite number, GW density
and constellation altitude on system coverage probability. Satellites follow
BPP, and GWs follow PPP. We assume that the transmitting power of the
satellite is 15dBw. The total coverage probability of the system is the
product of the coverage probability of two links: (i) satellite-GW link and
(ii) GW-user link. The satellite-GW link and GW-user link follow the same
channel fading, Large scale fading is 0dB in LoS range, Small scale fading
follows $\mathcal{SR}(1.29,0.158,19.4)$ [8], and the pass-loss exponent is 2.
When calculating the coverage probability, the threshold of SINR is -10dB.

Figure 4: Influence of satellite number and GW density on coverage
probability.

In Fig. 4, the satellites have an altitude of 1000km, that is, they are
distributed on a spherical surfaces with a radius of 7371km (radius of the
earth is 6371km). When satellites are insufficient, the system lacks available
satellites to provide coverage, or the associated satellites are too far away.
Conversely, too many satellites can cause significant interference. Therefore,
with the increase of the number of satellites, the coverage probability shows
a trend of first increasing and then decreasing. In addition, at an altitude
of 1000km, the optimal number of satellites for coverage probability is about
30. The optimal number of satellites is independent of the density of GWs.
Unlike satellites, increasing the density of GWs in a reasonable range will
improve the coverage performance of the system, and the GW subsystem shows the
characteristics of a noise-limited system. y

Figure 5: Influence of satellite number and constellation height on coverage
probability.

In Fig. 5, the density of the GW is fixed at 3/km2. The increase of
constellation height will reduce the coverage probability. There exists an
optimal number of satellites corresponding to altitude. With the increase of
constellation height, the number of satellites corresponding to the optimal
coverage decreases.

##  IV Challenges and Future scenarios

Although SG has many advantages in system performance analysis, the
traditional SG analysis framework has some limitations, and there are many
challenges in applying SG to more complex and sophisticated satellite systems
in the future. Therefore, in the right part of Fig. 1, we envision several
future scenarios. These scenarios can also be modeled, analyzed, and evaluated
using SG-based models as described in this article.

S-S Communication: In the previous analysis, we assumed that satellites were
independent of each other and analyzed the lower bound of the system
performance. But in future satellite systems, there will be more S-S
Communications, that will significantly improve the system’s performance.
Through interaction, satellites can learn more about the environment, such as
the number of nearby satellites and the frequency bands they use.
communicating with satellites in front of a predetermined trajectory allows
them to predict future communication environments. The problem of insufficient
capacity and achievable data rate can be effectively alleviated by smartly
designing association and hand over techniques between satellites. Conversely,
with the application of orthogonal channels, it is possible to have a high
SINR even if the user/GW is associated with a non-optimal satellite after the
handover.

Self-adapting to the environment: The second challenge faced by traditional SG
analysis frameworks is that, in reality, satellites do not remain inactive
when communications fail in harsh environments. While facing different
astronomical and geographical conditions, actively perceiving the environment
and selecting proper adaptive strategies are important. For example, the
satellite can always point the main lobe to the area with the greatest demand
for services, greatly improving local coverage. In the face of heavy rain
attenuation, satellites are expected to increase the transmission power to
maintain an acceptable SINR.

Virtual Satellite: In addition to S-S communications, there may be direct
cooperation between satellites. Similar to the idea of virtual hosting,
Multiple satellites with the same mission can be integrated into a virtual
satellite. The virtual satellite can decide the number of orthogonal channels
and allocate specific frequencies internally. What’s more, LEO satellites can
be used as the auxiliary of GEO satellites to complement the coverage of
regions that a single GEO cannot meet independently.

##  V Conclusions

This paper investigates the application of SG to simulating and analyzing LEO
satellite communications systems, introduces several satellite distributions
and analyzes the similarities and differences of channel models in the
literature. From a non-technical perspective, we describe coverage probability
in different scenarios, average achievable rate, and contact distance. This
paper analyzes several factors affecting latency and coverage probability
using SG-based models and gives the following simple conclusions. Increasing
the number of satellites and the height of the constellation within a specific
range can effectively reduce the latency. The ability of the satellites to
communicate with each other has a huge effect on reducing the value of latency
by 30 percent. This indicates that the recent technological advances in
enabling S-S communications will significantly widen the set of applications
that can rely on LEO satellite communications. For LEO satellites, there is an
optimal number of satellites to achieve the best coverage performance. The
optimal number of satellites decreases with constellation height and is
independent of the GW density. Increasing the density of the GW can
effectively improve the system coverage while increasing the height will
reduce the system coverage when the number of satellites remains unchanged.

## References

  * [1] I. Del Portillo, B. G. Cameron, and E. F. Crawley, “A technical comparison of three low earth orbit satellite constellation systems to provide global broadband,” _Acta Astronautica_ , vol. 159, pp. 123–135, 2019.
  * [2] S. Liu, Z. Gao, Y. Wu, D. W. K. Ng, X. Gao, K.-K. Wong, S. Chatzinotas, and B. Ottersten, “LEO satellite constellations for 5G and beyond: How will they reshape vertical domains?” 2021, available online: https://arxiv.org/abs/2106.0989.
  * [3] H. M. Mourad, A. A. M. Al-Bassiouni, S. S. Emam, and E. K. Al-Hussaini, “Generalized performance evaluation of low earth orbit satellite systems,” _IEEE communications letters_ , vol. 5, no. 10, pp. 405–407, 2001.
  * [4] M. Haenggi, _Stochastic Geometry for Wireless Networks_. Cambridge University Press, 2012.
  * [5] J. G. Andrews, F. Baccelli, and R. K. Ganti, “A tractable approach to coverage and rate in cellular networks,” _IEEE Transactions on communications_ , vol. 59, no. 11, pp. 3122–3134, 2011.
  * [6] N. Okati, T. Riihonen, D. Korpi, I. Angervuori, and R. Wichman, “Downlink coverage and rate analysis of low Earth orbit satellite constellations using stochastic geometry,” _IEEE Transactions on Communications_ , vol. 68, no. 8, pp. 5120–5134, 2020.
  * [7] A. Al-Hourani, “An analytic approach for modeling the coverage performance of dense satellite networks,” _IEEE Wireless Communications Letters_ , vol. 10, no. 4, pp. 897–901, 2021.
  * [8] A. Talgat, M. A. Kishk, and M.-S. Alouini, “Stochastic geometry-based analysis of LEO satellite communication systems,” _IEEE Communications Letters_ , vol. 25, no. 8, pp. 2458–2462, 2021.
  * [9] N. Okati and T. Riihonen, “Modeling and analysis of LEO mega-constellations as nonhomogeneous Poisson point processes,” in _IEEE 93rd Vehicular Technology Conference (VTC2021-Spring)_ , 2021, pp. 1–5.
  * [10] A. Talgat, M. A. Kishk, and M.-S. Alouini, “Nearest neighbor and contact distance distribution for binomial point process on spherical surfaces,” _IEEE Communications Letters_ , vol. 24, no. 12, pp. 2659–2663, 2020.
  * [11] B. A. Homssi and A. Al-Hourani, “Modeling uplink coverage performance in hybrid satellite-terrestrial networks,” 2021, available online: https://arxiv.org/abs/2106.04293.
  * [12] G. Alfano and A. De Maio, “Sum of squared shadowed-rice random variables and its application to communication systems performance prediction,” _IEEE Transactions on Wireless Communications_ , vol. 6, no. 10, pp. 3540–3545, 2007\.
  * [13] N. Okati and T. Riihonen, “Stochastic analysis of satellite broadband by mega-constellations with inclined LEOs,” in _31st Annual International Symposium on Personal, Indoor and Mobile Radio Communications_. IEEE, 2020, pp. 1–6.
  * [14] O. M. Bushnaq, M. A. Kishk, A. Celik, M.-S. Alouini, and T. Y. Al-Naffouri, “Optimal deployment of tethered drones for maximum cellular coverage in user clusters,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 3, pp. 2092–2108, 2020.
  * [15] F. Chiariotti, O. Vikhrova, B. Soret, and P. Popovski, “Information freshness of updates sent over LEO satellite multi-hop networks,” 2020, available online: https://arxiv.org/abs/2007.05449.

## biographies

Ruibo Wang is a Master candidate at King Abdullah University of Science and
Technology (KAUST). He received his B.S. degree in Communication Engineering
from University of Electronic Science and Technology of China in 2020.
---
Mustafa A. Kishk [S’16, M’18] is a postdoctoral research fellow at KAUST. He
received his Ph.D. degree in electrical engineering from Virginia Tech in
2018.
---
Mohamed-Slim Alouini [S’94, M’98, SM’03, F’09] is a professor at KAUST. He
received his Ph.D. degree in electrical engineering from California Institute
of Technology in 1998.
---
  *[SIR]: signal-to-interference ratio
  *[NLoS]: Non-Line-of-Sight
  *[LoS]: Line-of-Sight
  *[CCDF]: complementary cumulative distribution function
  *[SINR]: signal-to-interference plus noise ratio