Generalization of the Poisson distribution for the case of changing
probabilities of consecutive events

E.A.Kushnirenko†

###### Abstract

In this paper the generalization of the Poisson distribution is derived for
the case when each consecutive event changes event rate. A simple formula for
the probability of observing of a given number of events for the selected
period of time is derived for a given set of rates. Application of this
distribution in high-energy physics calculations is discussed.

Preface: This distribution was used to perform calculations of the backgrounds
in the electron-photon beam interactions [1],[2]. Mathematical properties of
the distribution are discussed here in more detail. This distribution can be
used when interaction causes change of the detector properties or of the
incident beam. Also it can be of certain interest as an example of yet another
generalization of the Poisson distribution.

##  1 Introduction

The Poisson distribution is used to describe a variety of different processes.
It is used very often in the analysis of data from the experimental sets in
accelerator physics, in cosmic rays physics, in radioactive decay studies, in
fluctuations of energy losses of a particle moving in the matter and in many
other cases. The only parameter that describes the distribution is the rate of
events $\mu$ or the average number of events per unit time. The assumption
that this parameter is constant gives the Poisson distribution unique features
and a certain beauty, however limits its applications.

Consider a several hundred $GeV$ electron beam going through a monochromatic
photon beam going in the opposite direction. An electron going through the
photon beam, can have several interactions with photons, and the energy of the
electron will drop down by an order of magnitude after 2 interactions. The
probability of the interaction of the electron with the monochromatic photon
beam of a particular energy is defined by the energy of the electron. In this
situation the rate $\mu$ of the Poisson process, which describes the
probability of the electron-photon interaction changes dramatically after each
interaction, for this reason even for an estimate of the number of the
electron-photon interactions the Poisson distribution can not be used.

Next we consider a process of shooting down a military aircraft with shells.
The probability of hitting the aircraft grows after each consecutive hit, as
its flying characteristics as well as defense systems go down and it becomes
an easier target to hit. There could be different scenarios which describe
variations of $\mu_{n}$. Clearly the classical Poison distribution will not
describe the process correctly and this problem requires generalization of the
Poisson distribution.

Also it is quite natural, that after the interaction the properties of the
detector change, which leads to the change of the rate of the detection of
events. For example the detection efficiency can degrade after the interaction
of the particle with the detector.

In these examples, the rate of observed events during the observation period
varies considerably, while during the period of constant rate $\mu$ the
observed process is the Poisson process. These examples show that it can be
useful to find distribution of events in processes where average rate changes.
In this paper we consider a process where the rate $\mu$ changes after each
consecutive event is happening, and the process for any length of time between
events is Poisson, with different rates $\mu_{i}$.

##  2 The generalization of the Poisson process

In general the problem can be formulated as follows: find a probability that
for a given period of time ${(0,t)}$ there would be $n$ consecutive events,
when the average rate of events $\mu$ is not constant, but changes
instantaneously after each consequent event. Before the first event the rate
is $\mu_{1}$, immediately after the first event and before the second event it
is $\mu_{2}$, immediately after the second event and before third event
$\mu_{3}$, immediately after event ${n-1}$ and before event $n$ the rate is
$\mu_{n}$. In each interval between consecutive events the rate is constant.

Generally speaking $\mu_{1},\mu_{2},\mu_{3},...,\mu_{n},....$ are independent
finite positive numbers. We will consider that all of these numbers are
different unless specifically stated otherwise.

The solution of the problem is the generalization of the Poisson distribution
which we look for. Obviously when ${\mu_{1}=\mu_{2}=...=\mu_{n}=...=\mu}$ the
distribution should transform into the Poisson distribution.

We will make the following assumptions:

  1. 1.

The probability ${P_{n}(t)}$ that during the period of time ${(0,t)}$ there
would be ${n\geq 1}$ events is defined by the rates
${\mu_{1},\mu_{2},...,\mu_{n+1}}$ which characterize the rates before
$1,2,\ldots,n$-th events and rate before event $n+1$.

  2. 2.

The probability of an event during the small period of time $\Delta t$ is
proportional to the duration $\Delta t$ and the expected rate $\mu_{i}$ of
events.

  3. 3.

The probability of 2 or more events during the small period of time $\Delta t$
is vanishingly small.

  4. 4.

For any $n\geq 1$ the event with number $n+1$ takes place later than the event
with number $n$.

Using these assumptions the probability $P_{0}(t+\Delta t)$ that during the
period of time $(0,t+\Delta t)$ there are no events means that there are no
events neither during the period $(0,t)$, nor during the period $(t,t+\Delta
t)$. For short periods of time $\Delta t$ this probability in accordance with
assumptions (2) and (3) is equal to:

$\displaystyle P_{0}(t+\Delta t)=P_{0}(t)\cdot(1-\mu_{1}\cdot\Delta t)$ (1)

where $P_{0}(t)$ and $(1-\mu_{1}\cdot\Delta t)$ are probabilities that there
are no events neither during the period of time $(0,t)$, nor during the period
of time $(t+\Delta t)$. With $\Delta t\rightarrow 0$ we have a differential
equation:

$\displaystyle\frac{dP_{0}(t)}{dt}=-\mu_{1}P_{0}(t)$ (2)

with the initial condition:

$\displaystyle P_{0}(0)=1$ (3)

which has the solution

$\displaystyle P_{0}(t)=e^{-\mu_{1}\cdot t}$ (4)

Let us now find out the probability $P_{1}(t+\Delta t)$ of single event
happening during the period of time $(t+\Delta t)$. In accordance with
assumptions (1),(2),(3) there could be 2 possibilities: either the first event
happens during the period $(0,t)$ and no events during the period $(t,t+\Delta
t)$, or there are no events during the period of time $(0,t)$ and single event
during the period $(t+\Delta t)$. So in accordance with assumptions (2),(3):

$\displaystyle P_{1}(t+\Delta t)=P_{1}(t)\cdot(1-\mu_{2}\cdot\Delta
t)+P_{0}(t)\cdot\mu_{1}\cdot\Delta t$ (5)

Considering the difference $P_{1}(t+\Delta t)-P_{1}(t)$ and taking $\Delta
t\rightarrow 0$ we obtain differential equation:

$\displaystyle\frac{dP_{1}(t)}{dt}=-\mu_{2}P_{1}(t)+\mu_{1}P_{0}(t)$ (6)

The solution of this equation using (Eq. 4) and initial condition

$\displaystyle P_{1}(0)=0$ (7)

is given by the formula:

$\displaystyle
P_{1}(t)=\mu_{1}\left[\frac{e^{-\mu_{1}t}}{\mu_{2}-\mu_{1}}+\frac{e^{-\mu_{2}t}}{\mu_{1}-\mu_{2}}\right]$
(8)

Analogous to the above one can solve the equation for 2 events with the
initial conditions

$P_{2}(0)=P_{3}(0)=0$ (9)

we find that:

$\displaystyle
P_{2}(t)=\mu_{1}\mu_{2}\left[\frac{e^{-\mu_{1}t}}{(\mu_{3}-\mu_{1})(\mu_{2}-\mu_{1})}+\right.$
$\displaystyle\left.+\frac{e^{-\mu_{2}t}}{(\mu_{3}-\mu_{2})(\mu_{1}-\mu_{2})}+\right.$
$\displaystyle\left.+\frac{e^{-\mu_{3}t}}{(\mu_{2}-\mu_{3})(\mu_{1}-\mu_{3})}\right]$
(10) $\displaystyle
P_{3}(t)=\mu_{1}\mu_{2}\mu_{3}\left[\frac{e^{-\mu_{1}t}}{(\mu_{4}-\mu_{1})(\mu_{3}-\mu_{1})(\mu_{2}-\mu_{1})}+\right.$
$\displaystyle\left.+\frac{e^{-\mu_{2}t}}{(\mu_{4}-\mu_{2})(\mu_{3}-\mu_{2})(\mu_{1}-\mu_{2})}+\right.$
$\displaystyle\left.+\frac{e^{-\mu_{3}t}}{(\mu_{4}-\mu_{3})(\mu_{2}-\mu_{3})(\mu_{1}-\mu_{3})}+\right.$
$\displaystyle\left.+\frac{e^{-\mu_{4}t}}{(\mu_{3}-\mu_{4})(\mu_{2}-\mu_{4})(\mu_{1}-\mu_{4})}\right]$
(11)

In general case the differential equation is

$\displaystyle\frac{dP_{n}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\mu_{n+1}P_{n}(t)+\mu_{n}P_{n-1}(t),{\rm when\,\,}n>0$
$\displaystyle\frac{dP_{n}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\mu_{n+1}P_{n}(t){\hskip 55.97205pt\rm when\,\,}n=0$ (12)

with initial condition

$\displaystyle P_{n}(0)$ $\displaystyle=$ $\displaystyle 0,{\rm when\,\,}n>0$
$\displaystyle P_{n}(0)$ $\displaystyle=$ $\displaystyle 1,{\rm when\,\,}n=0$
(13)

Comparing expressions for $P_{0}(t),P_{1}(t),P_{2}(t),P_{3}(t)$ and using
method of mathematical induction it is straightforward to show that the
solution for the probability $P_{n}(t)$ of $n$ events during the period
$(0,t)$

$\displaystyle
P_{n}(t)=\mu_{1}\mu_{2}...\mu_{n}\left[\frac{e^{-\mu_{1}t}}{(\mu_{n+1}-\mu_{1})(\mu_{n}-\mu_{1})...(\mu_{2}-\mu_{1})}+\right.$
$\displaystyle\left.+\frac{e^{-\mu_{2}t}}{(\mu_{n+1}-\mu_{2})(\mu_{n}-\mu_{2})...(\mu_{1}-\mu_{2})}+\right.\cdot\cdot\cdot$
$\displaystyle\left.\cdot\cdot\cdot+\frac{e^{-\mu_{n+1}t}}{(\mu_{n}-\mu_{n+1})(\mu_{n-1}-\mu_{n+1})...(\mu_{1}-\mu_{n+1})}\right]$
(14)

The structure of the obtained formula is fairly straightforward: $P_{n}(t)$
contains $n+1$ terms, where every $j$-th term $(1\leq j\leq n+1)$ contains a
fraction with $e^{-\mu_{j}t}$ in the numerator and the product
$(\mu_{n+1}-\mu_{j})(\mu_{n}-\mu_{j})\ldots(\mu_{1}-\mu_{j})$ with $n$ factors
but without $(\mu_{j}-\mu_{j})$ factor in the denominator.

Figure 1: The distributions $P_{n}(t)$ for $\mu_{1}=1$, $\mu_{2}=2$,
$\mu_{3}=3$, $\mu_{4}=0.5$, $\mu_{5}=1.50$, and $\mu_{6}=0$.

Expression $P_{n}(t)$ can be written in more compact form:

$\displaystyle
P_{n}(t)=\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot\sum\limits_{i=1}^{n+1}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}$ (15)

Remark: Clearly $P_{n}(t)$ is a function of not only time $t$, but also of
parameters $\mu_{1},\mu_{2},\ldots,\mu_{n},\mu_{n+1}$, so it would be more
appropriate to use notation $P_{n}(t,\mu_{1},\mu_{2},\ldots,\mu_{n+1})$.
Nevertheless we will use the notation $P_{n}(t)$ for the sake of shortness,
and will use the notation $P_{n}(t,\mu_{1},\mu_{2},\ldots,\mu_{n+1})$ in case
when it is absolutely necessary.

##  3 Analysis of the solution

###  3.1 Expression for the $m$-th derivative of $P_{n}(t)$
and some other useful formulas

Writing down differential equations (Eq. 12) for $n=1,2,3$ and taking into
account (Eq. 2) we obtain quite helpful for the further calculations system of
differential equations:

$\displaystyle\frac{dP_{0}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\mu_{1}P_{0}(t)$ $\displaystyle\frac{dP_{1}(t)}{dt}$
$\displaystyle=$ $\displaystyle-\mu_{2}P_{1}(t)+\mu_{1}P_{0}(t)$
$\displaystyle\frac{dP_{2}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\mu_{3}P_{2}(t)+\mu_{2}P_{1}(t)$ $\displaystyle\ldots$
$\displaystyle\frac{dP_{n}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\mu_{n+1}P_{n}(t)+\mu_{n}P_{n-1}(t)$ (16)

Then

$\displaystyle\frac{dP_{n}(t)}{dt}$ $\displaystyle=$
$\displaystyle(-1)^{1}\cdot\mu_{n+1}\cdot P_{n}(t)$ $\displaystyle+$
$\displaystyle(-1)^{0}\cdot\mu_{n}\cdot P_{n-1}(t)$
$\displaystyle\frac{d^{2}P_{n}(t)}{dt^{2}}$ $\displaystyle=$
$\displaystyle(-1)^{2}\cdot\mu_{n+1}^{2}\cdot P_{n}(t)$ $\displaystyle+$
$\displaystyle(-1)^{1}\cdot\mu_{n}\left(\mu_{n+1}+\mu_{n}\right)\cdot
P_{n-1}(t)$ $\displaystyle+$ $\displaystyle(-1)^{0}\cdot\mu_{n}\mu_{n-1}\cdot
P_{n-2}(t)$ $\displaystyle\frac{d^{3}P_{n}(t)}{dt^{3}}$ $\displaystyle=$
$\displaystyle(-1)^{3}\cdot\mu_{n+1}^{3}\cdot P_{n}(t)$ (17) $\displaystyle+$
$\displaystyle(-1)^{2}\cdot\mu_{n}\left(\mu_{n+1}^{2}+\mu_{n+1}\mu_{n}+\mu_{n}^{2}\right)\cdot
P_{n-1}(t)$ $\displaystyle+$
$\displaystyle(-1)^{1}\cdot\mu_{n}\mu_{n-1}\left(\mu_{n+1}+\mu_{n}+\mu_{n-1}\right)\cdot
P_{n-2}(t)$ $\displaystyle+$
$\displaystyle(-1)^{0}\cdot\mu_{n}\mu_{n-1}\mu_{n-2}\cdot P_{n-3}(t)$

General formula for the $m$-th derivative of $P_{n}(t)$, where $(n\geq 0)$,
which includes system of equations (Eq. 16) is given by:

$\displaystyle\frac{d^{m}P_{n}(t)}{dt^{m}}$ $\displaystyle=$
$\displaystyle(-1)^{m}\mu_{n+1}^{m}P_{n}(t)+$ (18) $\displaystyle+$
$\displaystyle\sum\limits_{s=1}^{m}(-1)^{m-s}\cdot\prod\limits_{j=n-s+1}^{n}\mu_{j}\cdot\left(\sum\limits_{i=n-s+1}^{n+1}\mu_{i}\right)^{\otimes(m-s)}\cdot
P_{n-s}(t)$

where the following notation is used:

$\displaystyle\left(\sum\limits_{i=1}^{r}\mu_{i}\right)^{\otimes
n}=\sum\limits_{k_{1}+k_{2}+\ldots+k_{r}=n}\mu_{1}^{k_{1}}\cdot\mu_{2}^{k_{2}}\cdot\ldots\cdot\mu_{r}^{k_{r}},$
(19)

For example:

$\displaystyle\left(\mu_{1}+\mu_{2}\right)^{\otimes 3}$ $\displaystyle=$
$\displaystyle\mu_{1}^{3}+\mu_{1}^{2}\cdot\mu_{2}+\mu_{1}\cdot\mu_{2}^{2}+\mu_{2}^{3}$
$\displaystyle\left(\mu_{1}+\mu_{2}+\mu_{3}\right)^{\otimes 2}$
$\displaystyle=$
$\displaystyle\mu_{1}^{2}+\mu_{1}\mu_{2}+\mu_{2}^{2}+\mu_{2}\mu_{3}+\mu_{1}\mu_{3}+\mu_{3}^{2}$

Also in usual notation we have:

$\displaystyle\left(\sum\limits_{i=1}^{3}\mu_{i}\right)^{\otimes
n}=\sum\limits_{s=0}^{n}\sum\limits_{m=0}^{s}\mu_{1}^{s-m}\mu_{2}^{m}\mu_{3}^{n-s}=\sum\limits_{s=0}^{n}\sum\limits_{m=0}^{n-s}\mu_{1}^{n-s-m}\mu_{2}^{m}\mu_{3}^{s}$
(20)

Expression for the derivative (Eq.18) was derived for $n\geq 0$ and $m<n$, but
it can be expanded for any $n$ or $m$ if we assume that:

$\displaystyle\left(\sum\limits_{i=1}^{r}\mu_{i}\right)^{\otimes n}$
$\displaystyle=$ $\displaystyle 1,\mbox{~{}when~{}}n=0,$
$\displaystyle\left(\sum\limits_{i=1}^{r}\mu_{i}\right)^{\otimes n}$
$\displaystyle=$ $\displaystyle 0,\mbox{~{}when~{}}n<0.$ (21) $\displaystyle
P_{n}(t)$ $\displaystyle=$ $\displaystyle 0,\mbox{~{}when~{}}n<0$

It is interesting to note the behaviour of (Eq. 18) at $t=0$. First of all at
$t=0$

$\displaystyle P_{n}(0)=0\mbox{ for }n\geq 1$ (22) $\displaystyle
P_{n}(0)=1\mbox{ for }n=0$ (23)

Then the first term in (Eq. 18) is equal to 0 for $n\geq 1$. In the second
term all the terms of the sum vanish except for the term with $s=n$. According
to (Eq. 21) we have non-zero terms only when $m\geq n$. So, given that
$P_{0}(0)=1$, the formula (Eq. 18) at $t=0$ transforms into:

$\displaystyle\left.\frac{d^{m}P_{n}(t)}{dt^{m}}\right|_{t=0}=(-1)^{(m-n)}\cdot\mathop{\displaystyle\coprod}\limits_{j=1}^{n}\mu_{j}\cdot\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes(m-n)}\mbox{~{}when~{}}m\geq
n.$ (24)

Let us use $m=n+r$ for convenience, then the above formula transforms into

$\displaystyle\left.\frac{d^{n+r}P_{n}(t)}{dt^{n+r}}\right|_{t=0}=(-1)^{r}\cdot\mathop{\displaystyle\coprod}\limits_{j=1}^{n}\mu_{j}\cdot\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes(r)}\mbox{~{}when~{}}r\geq
0.$ (25)

From expression (Eq. 18) given the remark (Eq. 21) we conclude that at $t=0$:
the derivatives :

$\displaystyle\left.\frac{d^{m}P_{n}(t)}{dt^{m}}\right|_{t=0}=0,\mbox{~{}when~{}}m<n$
(26)

According to (Eq. 25) the first derivative which does not vanish at $t=0$ of
the function $P_{n}(t)$ has the rank $m=n$ and it is equal to

$\displaystyle\left.\frac{d^{n}P_{n}(t)}{dt^{n}}\right|_{t=0}=\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}$
(27)

So we can combine the results in the following form:

$\displaystyle\left.\frac{d^{m}P_{n}(t)}{dt^{m}}\right|_{t=0}$
$\displaystyle=$
$\displaystyle(-1)^{(m-n)}\cdot\mathop{\displaystyle\coprod}\limits_{j=1}^{n}\mu_{j}\cdot\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes(m-n)}\mbox{~{}when~{}}m>n.$
(28) $\displaystyle\left.\frac{d^{m}P_{n}(t)}{dt^{m}}\right|_{t=0}$
$\displaystyle=$
$\displaystyle\mathop{\displaystyle\coprod}\limits_{j=1}^{n}\mu_{j},\mbox{~{}when~{}}m=n$
(29) $\displaystyle\left.\frac{d^{m}P_{n}(t)}{dt^{m}}\right|_{t=0}$
$\displaystyle=$ $\displaystyle 0,\mbox{~{}when~{}}m<n$ (30)

taking derivative $d^{m}/dt^{m}$ in (Eq. 15) we obtain that

$\displaystyle
P_{n}(t)=\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot\sum\limits_{i=1}^{n+1}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}$
$\displaystyle\frac{d^{m}P_{n}(t)}{dt^{m}}=\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot\sum\limits_{i=1}^{n+1}\frac{(-1)^{m}\cdot\mu_{i}^{m}\cdot
e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}$ (31)

Then for any given $m<n$ the sum of the set vanish:

$\displaystyle\sum\limits_{i=1}^{n+1}\frac{(-1)^{m}\cdot\mu_{i}^{m}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}=0$ (32)

Given the condition $P_{n}(0)=0$ for $n>0$ the following set vanishes:

$\displaystyle\sum\limits_{i=1}^{n+1}\frac{1}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}=0,\mbox{~{}when~{}}n>0,$ (33)

which is a special case of (Eq. 32) when $m=0$

Taking the derivative $d^{n}/dt^{n}$ of the (Eq. 15) we conclude that the sum
of the following set is equal to 1:

$\displaystyle\sum\limits_{i=1}^{n+1}\frac{(-1)^{n}\cdot\mu_{i}^{n}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}=1$ (34)

In the general case, which applies to any $m$ and $n$ we have:

$\displaystyle\sum\limits_{i=1}^{n+1}\frac{(-1)^{m}\cdot\mu_{i}^{m}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}=(-1)^{m-n}\cdot\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes(m-n)},$
(35)

assuming (Eq. 21) for ($m\leq n$)

Taking into account (Eq. 24) and (Eq. 27) the Taylor series expansion in the
vicinity of $t=0$ can be written as:

$\displaystyle P_{n}(t)$ $\displaystyle=$ $\displaystyle
P_{n}(0)+\sum\limits_{m=1}^{\infty}\frac{t^{m}}{m!}\cdot\left.\frac{d^{m}P_{n}(t)}{dt^{m}}\right|_{t=0}=P_{n}(0)+\sum\limits_{m=n}^{\infty}\frac{t^{m}}{m!}\cdot\left.\frac{d^{m}P_{n}(t)}{dt^{m}}\right|_{t=0}$
$\displaystyle=$ $\displaystyle
P_{n}(0)+\sum\limits_{r=0}^{\infty}\frac{t^{n+r}}{(n+r)!}\cdot\left.\frac{d^{n+r}P_{n}(t)}{dt^{n+r}}\right|_{t=0}$
$\displaystyle=$ $\displaystyle
P_{n}(0)+\left(\mathop{\displaystyle\coprod}\limits_{j=1}^{n}\mu_{j}\right)\cdot\sum\limits_{r=0}^{\infty}(-1)^{r}\cdot\frac{t^{n+r}}{(n+r)!}\cdot\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes(r)}$
$\displaystyle=$ $\displaystyle
P_{n}(0)+\left(\mathop{\displaystyle\coprod}\limits_{j=1}^{n}\mu_{j}\right)\cdot\frac{t^{n}}{n!}\sum\limits_{r=0}^{\infty}(-1)^{r}\cdot\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes(r)}\cdot\frac{t^{r}\cdot
n!}{(n+r)!}$

###  3.2 The limiting transition of
$P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{n+1},t)$ for the case when
$(\mu_{1},\mu_{2},\ldots,\mu_{n+1})\rightarrow\mu$

For $n>0$ let us consider $P_{n}(t)$ in the form

$\displaystyle
P_{n}(t)=\left(\mathop{\displaystyle\prod}\limits_{j=1}^{n}\mu_{j}\right)\cdot\frac{t^{n}}{n!}\cdot\left[\sum\limits_{r=0}^{\infty}(-1)^{r}\frac{n!}{(n+r)!}\cdot
t^{r}\cdot\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes r}\right]$ (37)

Let all $\mu_{i}=\mu$. According to the notation (Eq. 19):

$\displaystyle\left(\sum\limits_{i=1}^{n+1}\mu_{i}\right)^{\otimes
r}=\sum\limits_{k_{1}+k_{2}+\ldots+k_{n+1}=r}\mu_{1}^{k_{1}}\cdot\mu_{2}^{k_{2}}\cdot\ldots\cdot\mu_{n+1}^{k_{n+1}}$
(38)

The number $N$ of terms in the expression (Eq. 38) is given by the formula:

$\displaystyle N=C_{n+r}^{r}=\frac{(n+r)!}{n!\cdot r!}.$ (39)

Hence using formula (Eq. 37)

$\displaystyle\lim_{\mu_{i}\rightarrow\mu}P_{n}(t)$ $\displaystyle=$
$\displaystyle\frac{\mu^{n}t^{n}}{n!}\sum\limits_{r=0}^{\infty}(-1)^{r}\cdot\frac{t^{r}}{r!}\cdot\mu^{r}$
(40) $\displaystyle=$ $\displaystyle\frac{\left(\mu\cdot
t\right)^{n}}{n!}\cdot e^{-\mu\cdot t},\mbox{~{}when~{}}n>1$
$\displaystyle\lim_{\mu_{1}\rightarrow\mu}P_{0}(t)$ $\displaystyle=$
$\displaystyle e^{-\mu t},\mbox{~{}when~{}}n=0$ (41)

So for any $n\geq 0$ generalized Poisson distribution transforms into Poisson
distribution when $(\mu_{1},\mu_{2},\ldots,\mu_{n+1})\rightarrow\mu$

###  3.3 Case of $\mu_{n+1}=0$

It is quite clear that if $\mu_{n+1}=0$, which basically means that there are
only $n$ events, and event $n+1$ would never happen, then $P_{n}(t)$ is
assimptotically growing to 1, and $P_{n+1}(t)=P_{n+2}(t)=\ldots=0$. Let us
consider specific cases:

$\mu_{1}>0,\mu_{2}=0$, then

$P_{1}(t)=\mu_{1}\left[\frac{e^{-\mu_{1}t}}{\mu_{2}-\mu_{1}}+\frac{e^{-\mu_{2}t}}{\mu_{1}-\mu_{2}}\right]=1-e^{-\mu_{1}t}=1-P_{0}(t)$
(42)

$\mu_{1,2}>0,\mu_{3}=0$, then

$\displaystyle P_{2}(t)$ $\displaystyle=$
$\displaystyle\mu_{1}\mu_{2}\left[\frac{e^{-\mu_{1}t}}{(-\mu_{1})(\mu_{2}-\mu_{1})}+\frac{e^{-\mu_{2}t}}{(-\mu_{2})(\mu_{1}-\mu_{2})}+\frac{1}{\mu_{2}\cdot\mu_{1}}\right]$
(43) $\displaystyle=$ $\displaystyle
1-\mu_{1}\frac{e^{-\mu_{2}t}}{\mu_{1}-\mu_{2}}-\mu_{2}\frac{e^{-\mu_{1}t}}{\mu_{2}-\mu_{1}}$
$\displaystyle=$ $\displaystyle
1-P_{1}(t)+\mu_{1}\frac{e^{-\mu_{1}t}}{\mu_{2}-\mu_{1}}-\mu_{2}\frac{e^{-\mu_{1}t}}{\mu_{2}-\mu_{1}}$
$\displaystyle=$ $\displaystyle 1-P_{1}(t)-P_{0}(t)$

It is straightforward to show than when $\mu_{1,\ldots,n-1}>0,\mu_{n}=0$ then

$\displaystyle P_{n-1}(t)$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\prod}\limits_{i=1}^{n-1}\mu_{i}\cdot\sum\limits_{i=1}^{n}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n}(\mu_{j}-\mu_{i})}$ (44) $\displaystyle=$ $\displaystyle
1+\mathop{\displaystyle\prod}\limits_{i=1}^{n-1}\mu_{i}\cdot\sum\limits_{i=1}^{n-1}\frac{e^{-\mu_{i}t}}{(-\mu_{i})\cdot\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n-1}(\mu_{j}-\mu_{i})}$ $\displaystyle=$ $\displaystyle
1-P_{n-2}(t)-\ldots-P_{1}(t)-P_{0}(t)$

Indeed assuming (Eq. 44) being true, we can derive that, when
$\mu_{1,\ldots,n}>0$, and $\mu_{n+1}=0$ then

$\displaystyle P_{n}(t)$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot\sum\limits_{i=1}^{n+1}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}$ (45) $\displaystyle=$ $\displaystyle
1+\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot\sum\limits_{i=1}^{n}\frac{e^{-\mu_{i}t}}{(-\mu_{i})\cdot\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n}(\mu_{j}-\mu_{i})}$ $\displaystyle=$ $\displaystyle
1-\mathop{\displaystyle\prod}\limits_{i=1}^{n-1}\mu_{i}\cdot\frac{e^{-\mu_{n}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n}(\mu_{j}-\mu_{i})}+\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot\sum\limits_{i=1}^{n-1}\frac{e^{-\mu_{i}t}}{(-\mu_{i})\cdot\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n}(\mu_{j}-\mu_{i})}$ $\displaystyle=$ $\displaystyle
1-P_{n-1}(t)+\mathop{\displaystyle\prod}\limits_{i=1}^{n-1}\mu_{i}\cdot\sum\limits_{i=1}^{n-1}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n}(\mu_{j}-\mu_{i})}+\mathop{\displaystyle\prod}\limits_{i=1}^{n-1}\mu_{i}\cdot\sum\limits_{i=1}^{n-1}\frac{\mu_{n}\cdot
e^{-\mu_{i}t}}{(-\mu_{i})\cdot\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n}(\mu_{j}-\mu_{i})}$ $\displaystyle=$ $\displaystyle
1-P_{n-1}(t)+\mathop{\displaystyle\prod}\limits_{i=1}^{n-1}\mu_{i}\cdot\sum\limits_{i=1}^{n-1}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n}(\mu_{j}-\mu_{i})}\left[1+\frac{\mu_{n}}{-\mu_{i}}\right]$
$\displaystyle=$ $\displaystyle
1-P_{n-1}(t)+\mathop{\displaystyle\prod}\limits_{i=1}^{n-1}\mu_{i}\cdot\sum\limits_{i=1}^{n-1}\frac{e^{-\mu_{i}t}}{(-\mu_{i})\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n-1}(\mu_{j}-\mu_{i})}{\mathrm{\;\;hence\,\,using\,(Eq.\ref{eq:p_n-1_when0})}}$
$\displaystyle=$ $\displaystyle 1-P_{n-1}(t)-P_{n-2}(t)-\ldots-
P_{1}(t)-P_{0}(t)$

###  3.4 The Integral of $P_{n}(t)$

There are several ways to prove that:

$\displaystyle\mu_{n+1}\cdot\int\limits_{0}^{\infty}P_{n}(t)\cdot dt=1$ (46)

Indeed using direct integration and (Eq. 45)

$\displaystyle\int\limits_{0}^{\infty}P_{n}(t)\cdot dt$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot\sum\limits_{i=1}^{n+1}\frac{1}{\mu_{i}\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}=\frac{1}{\mu_{n+1}}\mathop{\displaystyle\prod}\limits_{i=1}^{n+1}\mu_{i}\cdot\sum\limits_{i=1}^{n+1}\frac{1}{\mu_{i}\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}$ (47) $\displaystyle=$
$\displaystyle\frac{1}{\mu_{n+1}}\left[P_{n}(0)+P_{n-1}(0)+\ldots+P_{1}(0)+P_{0}(0)\right]=\frac{1}{\mu_{n+1}}$

Similarly as $\mu_{1}\int\limits_{0}^{\infty}P_{0}(t)\cdot dt=1$ then using
mathematical induction applied to (Eq. 12) one can prove that

$\displaystyle\frac{dP_{n}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\mu_{n+1}P_{n}(t)+\mu_{n}P_{n-1}(t)$
$\displaystyle\mu_{n+1}\int\limits_{0}^{\infty}P_{n}(t)\cdot dt$
$\displaystyle=$ $\displaystyle
P_{n}(0)-P_{n}(\infty)+\mu_{n}\cdot\int\limits_{0}^{\infty}P_{n-1}(t)\cdot
dt=1$ (48)

###  3.5 Normalisation of $P_{n}(t)$

Let us prove that the limit of the sum
$S(t)=\sum\limits_{n=0}^{\infty}P_{n}(t)$ is equal to 1, at any given time
$t$.

In the special case when one of the $\mu_{i}=0$ the formula (Eq. 45) can be
applied:

$\displaystyle\sum\limits_{i=0}^{n-1}P_{i}(t)$ $\displaystyle=$ $\displaystyle
1-P_{n}(t)$ $\displaystyle P_{n+1}$ $\displaystyle=$ $\displaystyle
P_{n+2}=\ldots=0,$

then clearly $\sum\limits_{n=0}^{\infty}P_{n}(t)=1$

For the case when every $\mu_{i}>0$ let us consider the sum:

$\displaystyle S_{n}(t)=\sum\limits_{i=0}^{n}P_{i}(t)$ (49)

Let us find the limit, to which the first derivative $S_{n}(t)$ is
approximating when $n$ is growing.

$\displaystyle\lim_{n\rightarrow\infty}\frac{dS_{n}(t)}{dt}=\sum\limits_{i=0}^{\infty}\frac{dP_{i}(t)}{dt}$
(50)

Considering the system of equations (Eq. 16) and (Eq. 46)

$\displaystyle\sum\limits_{i=0}^{n}\frac{dP_{i}(t)}{dt}=-\mu_{n+1}\cdot
P_{n}(t)=-\frac{\int\limits_{0}^{t}P_{n}(t)\cdot
dt}{\int\limits_{0}^{\infty}P_{n}(t)\cdot dt}$ (51)

Clearly if we fix $t$ and keep increasing $n$ the fraction in (Eq. 51)
vanishes. So for any given $t\geq 0$ and for any given set $\mu_{i}>0$:

$\displaystyle\frac{dS(t)}{dt}=\lim_{n\rightarrow\infty}\frac{dS_{n}(t)}{dt}=0$
(52)

and hence $S(t)=S(0)=\mbox{const}$

So, taking into account (Eq. 13) we derive that $S_{n}(0)=1$, for $n>0$. That
means that for any $t\geq 0$ and for any given set $\mu_{i}$ the normalization
of $P_{n}(t)$ is:

$\displaystyle\sum\limits_{i=0}^{\infty}P_{i}(t)=S(t)=S(0)=1$ (53)

###  3.6 The symmetry of $P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{n+1},t)$ under the
permutation of parameters $\mu_{i},\mu_{j}$

One of very obvious properties of functions $P_{n}(t)$ is its independence
from mutual interchange of parameters $\mu_{1},\mu_{2},\ldots,\mu_{n}$,
because of terms simply interchange in (Eq. 15).

$\displaystyle
P_{n}(t)=\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot\sum\limits_{i=1}^{n+1}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{k=1;k\neq
i}^{n+1}(\mu_{k}-\mu_{i})}$

This gives us the possibility to put parameters $\mu_{i}$ in any order, for
example in the order of growing $\mu_{i}$.

Let us remark, that this does not apply to parameter $\mu_{n+1}$. In this case
when we interchange $\mu_{i}$ with $\mu_{n+1}$, the distribution changes
according to:
$P_{n}(\mu_{1},\ldots,\mu_{i},\ldots,\mu_{n},\mu_{n+1},t)\cdot\mu_{n+1}=P_{n}(\mu_{1},\ldots,\mu_{n+1},\ldots,\mu_{n},\mu_{i},t)\cdot\mu_{i}$

So the distribution changes under the permutation of parameters
$\mu_{i},\mu_{j}$ according to:

$\displaystyle
P_{n}(\mu_{1},\ldots,\mu_{i},\mu_{j},\ldots,\mu_{n},\mu_{n+1},t)$
$\displaystyle=$ $\displaystyle
P_{n}(\mu_{1},\ldots,\mu_{j},\mu_{i},\ldots,\mu_{n},\mu_{n+1},t)$ (54)
$\displaystyle P_{n}(\mu_{1},\ldots,\mu_{i},\ldots,\mu_{n},\mu_{n+1},t)$
$\displaystyle=$ $\displaystyle
P_{n}(\mu_{1},\ldots,\mu_{n+1},\ldots,\mu_{n},\mu_{i},t)\frac{\mu_{i}}{\mu_{n+1}}$
(55)

Given this property when studying properties of function $P_{n}(t)$ it is
usually enough to study dependency
$P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{n},\mu_{n+1},t)$ with respect to
$\mu_{n+1}$ and with respect to any of the parameters
$\mu_{1},\ldots,\mu_{n}$, and we can use $\mu_{1}$ for example.

###  3.7 The limiting transition in
$P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{n+1},t)$ when several parameters $\mu_{i}$
are equal to each other

.

Let us assume that $k$ different parameters $\mu_{i}$ are equal to each other.
Due to the symmetry of $P_{n}$ under the permutation of parameters $\mu_{i}$
we can calculate the case of
$P_{n}(\mu_{1}=\mu_{2}=\ldots=\mu_{k}=\mu,\mu_{k+1},\ldots,\mu_{n+1})$

Let us consider $k=2$, then $\mu_{1}=\mu$ and we calculate the limit when
$\mu_{2}\rightarrow\mu$:

$\displaystyle\frac{P_{n}}{\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}}$
$\displaystyle=$ $\displaystyle\frac{e^{-\mu\cdot
t}}{(\mu_{2}-\mu)\ldots(\mu_{n+1}-\mu)}+\frac{e^{-\mu_{2}\cdot
t}}{(\mu-\mu_{2})\ldots(\mu_{n+1}-\mu_{2})}+\sum\limits_{i=3}^{n+1}\frac{e^{-\mu_{i}\cdot
t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq i}^{n+1}(\mu_{j}-\mu_{i})}$

Let us define the third term as $R(3)$ and define $\Delta\mu=\mu_{2}-\mu$,
then expanding in terms power of $\Delta\mu$

$\displaystyle\frac{P_{n}}{\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}}$
$\displaystyle=$ $\displaystyle\frac{e^{-\mu\cdot
t}}{\Delta\mu(\mu_{3}-\mu)\ldots(\mu_{n+1}-\mu)}\left[1+\frac{e^{-\Delta\mu\cdot
t}}{(-1)(1-\frac{\Delta\mu}{\mu_{3}-\mu})\ldots(1-\frac{\Delta\mu}{\mu_{n+1}-\mu})}\right]+R(3)$
(57) $\displaystyle=$ $\displaystyle\frac{e^{-\mu\cdot
t}}{\Delta\mu(\mu_{3}-\mu)\ldots(\mu_{n+1}-\mu)}\left[1-\left(1+(-\Delta\mu\cdot
t)+\frac{(-\Delta\mu\cdot t)^{2}}{2!}+\ldots\right)\times\right.$
$\displaystyle\times$
$\displaystyle\left.\mathop{\displaystyle\prod}\limits_{i=3}^{n+1}\left(1+\frac{\Delta\mu}{\mu_{i}-\mu}+\left(\frac{\Delta\mu}{\mu_{i}-\mu}\right)^{2}+\ldots\right)\right]+R(3)$
$\displaystyle=$ $\displaystyle\frac{e^{-\mu\cdot
t}}{\Delta\mu(\mu_{3}-\mu)\ldots(\mu_{n+1}-\mu)}\left[1-\left(1+\left[\sum_{i=3}^{n+1}\frac{\Delta\mu}{\mu_{i}-\mu}-\Delta\mu\cdot
t\right]+\right.\right.$ $\displaystyle+$
$\displaystyle\left.\left.\left[\frac{(-\Delta\mu\cdot
t)^{2}}{2!}+(-\Delta\mu\cdot
t)\left(\sum_{i=3}^{n+1}\frac{\Delta\mu}{\mu_{i}-\mu}\right)^{\otimes
1}+\left(\sum_{i=3}^{n+1}\frac{\Delta\mu}{\mu_{i}-\mu}\right)^{\otimes
2}\right]+\ldots\right.\right.$ $\displaystyle+$
$\displaystyle\left.\left.\sum_{k=0}^{n}\left[\frac{(-\Delta\mu\cdot
t)^{k}}{k!}\left(\sum_{i=3}^{n+1}\frac{\Delta\mu}{\mu_{i}-\mu}\right)^{\otimes{n-k}}\right]+\ldots\right)\right]+R(3)$

where notation (Eq. 19) was used. In the limit $\Delta\mu\rightarrow 0$

$\displaystyle\frac{P_{n}(\mu_{1},\mu_{2}=\mu)}{\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}}$
$\displaystyle=$ $\displaystyle\frac{e^{-\mu\cdot
t}}{\mathop{\displaystyle\prod}\limits_{i=3}^{n+1}(\mu_{i}-\mu)}\left[t+\sum_{i=3}^{n+1}\frac{-1}{\mu_{i}-\mu}\right]+R(3)$
$\displaystyle\frac{P_{n}(\mu_{1},\mu_{2},\mu_{3}=\mu)}{\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}}$
$\displaystyle=$ $\displaystyle\frac{e^{-\mu\cdot
t}}{\mathop{\displaystyle\prod}\limits_{i=4}^{n+1}(\mu_{i}-\mu)}\left[\frac{t^{2}}{2!}+t\cdot\sum_{i=4}^{n+1}\frac{-1}{\mu_{i}-\mu}+\left(\cdot\sum_{i=4}^{n+1}\frac{-1}{\mu_{i}-\mu}\right)^{\otimes
2}\right]+R(4)$
$\displaystyle\frac{P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{k}=\mu)}{\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}}$
$\displaystyle=$ $\displaystyle\frac{e^{-\mu\cdot
t}}{\mathop{\displaystyle\prod}\limits_{i=k+1}^{n+1}(\mu_{i}-\mu)}\left[\sum_{j=0}^{k-1}\frac{t^{k-j}}{(k-j)!}\times\left(\sum_{i=k+1}^{n+1}\frac{-1}{\mu_{i}-\mu}\right)^{\otimes
j}\right]+R(k+1)$

where we use the notation

$R(k)=\sum\limits_{i=k}^{n+1}\frac{e^{-\mu_{i}\cdot
t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq i}^{n+1}(\mu_{j}-\mu_{i})}$
(59)

In a particular case when only $\mu_{n+1}$ differs from the rest of
$\mu_{i}=\mu$ we obtain:

$\displaystyle P_{2}(\mu_{1},\mu_{2}=\mu,\mu_{3})$ $\displaystyle=$
$\displaystyle\mu^{2}\left(\frac{e^{-\mu_{3}\cdot
t}}{(\mu-\mu_{3})^{2}}-e^{-\mu\cdot
t}\left[\frac{1}{(\mu-\mu_{3})^{2}}+\frac{t}{(\mu-\mu_{3})}\right]\right)$
$\displaystyle P_{3}(\mu_{1},\mu_{2},\mu_{3}=\mu,\mu_{4})$ $\displaystyle=$
$\displaystyle\mu^{3}\left(\frac{e^{-\mu_{4}\cdot
t}}{(\mu-\mu_{4})^{3}}-e^{-\mu\cdot
t}\left[\frac{1}{(\mu-\mu_{4})^{3}}+\frac{t}{(\mu-\mu_{4})^{2}}+\frac{t^{2}}{2!(\mu-\mu_{4})}\right]\right)$
$\displaystyle P_{n}(\mu_{1},\ldots,\mu_{n}=\mu,\mu_{n+1})$ $\displaystyle=$
$\displaystyle\mu^{n}\left(\frac{e^{-\mu_{n+1}\cdot
t}}{(\mu-\mu_{n+1})^{n}}-e^{-\mu\cdot
t}\left[\sum\limits_{j=0}^{n-1}\frac{t^{j}}{j!\left(\mu-\mu_{n+1}\right)^{(n-j)}}\right]\right)$
(60)

###  3.8 The limiting transition in
$P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{n+1},t)$ for the case when
$\mu_{i}\rightarrow\infty$

It is easy to see from compact form of the distribution (Eq. 15):

$\displaystyle
P_{n}(t)=\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot\sum\limits_{i=1}^{n+1}\frac{e^{-\mu_{i}t}}{\mathop{\displaystyle\prod}\limits_{j=1;j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}$

that

$\displaystyle\lim_{\mu_{1}\rightarrow\infty}P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{n+1},t)$
$\displaystyle=$
$\displaystyle\left(\mathop{\displaystyle\prod}\limits_{i=2}^{n}\mu_{i}\right)\cdot\sum\limits_{i=2}^{n+1}\frac{e^{-\mu_{i}\cdot
t}}{\mathop{\displaystyle\prod}\limits_{j=2,j\neq
i}^{n+1}(\mu_{j}-\mu_{i})}\,\rm{,or}$
$\displaystyle\lim_{\mu_{1}\rightarrow\infty}P_{n}(\mu_{1},\mu_{2},\ldots,\mu_{n+1},t)$
$\displaystyle=$ $\displaystyle P_{n-1}(\mu_{2},\ldots,\mu_{n+1},t)$

which certainly makes sense, as in the limit $\mu_{1}\rightarrow\infty$ 1-st
event happens right away and following events obey $P_{n-1}(t)$ distribution.
Due to the symmetry under the permutations (Eq. 54) same formula applies to
any $\mu_{i}$ where $i\leq n$.

When the last parameter $\mu_{n+1}\rightarrow\infty$, then using (Eq. 55)

$\displaystyle\lim_{\mu_{n+1}\rightarrow\infty}P_{n}(\mu_{1},\ldots,\mu_{n},\mu_{n+1},t)$
$\displaystyle=$
$\displaystyle\lim_{\mu_{n+1}\rightarrow\infty}P_{n}(\mu_{1},\ldots,\mu_{n+1},\mu_{n},t)\frac{\mu_{n}}{\mu_{n+1}}$
$\displaystyle=$ $\displaystyle
P_{n-1}(\mu_{1},\ldots,\mu_{n},t)\frac{\mu_{n}}{\mu_{n+1}}$

we can unify the result in the following way:

$\displaystyle\lim_{\mu_{i}\rightarrow\infty}P_{n}(t)$ $\displaystyle=$
$\displaystyle
P_{n-1}(\mu_{1},\ldots,\mu_{i-1},\mu_{i+1},\ldots,\mu_{n+1},t)\,\,{\rm
when}\,i\leq n$ (61) $\displaystyle\lim_{\mu_{i}\rightarrow\infty}P_{n}(t)$
$\displaystyle=$ $\displaystyle
P_{n-1}(\mu_{1},\ldots,\mu_{n},t)\frac{\mu_{n}}{\mu_{n+1}}\,\,{\rm
when}\,i=n+1$ (62)

###  3.9 Partial derivative $\partial P_{n}/\partial\mu_{i}$

Let us first calculate $\partial P_{n}/\partial\mu_{n+1}$ and for any other
$\mu_{i}$ the partial derivative can be calculated using the symmetry (Eq. 54)
and (Eq. 55) under the permutation of the parameters $\mu_{i}$.

$\displaystyle\frac{\partial P_{1}(\mu_{1},\mu_{2},t)}{\partial\mu_{2}}$
$\displaystyle=$
$\displaystyle-\mu_{1}\exp(-\mu_{2}t)\cdot\left[\frac{e^{-(\mu_{1}-\mu_{2})t}+(\mu_{1}-\mu_{2})t-1}{(\mu_{2}-\mu_{1})^{2}}\right]$
$\displaystyle\frac{\partial
P_{2}(\mu_{1},\mu_{2},\mu_{3},t)}{\partial\mu_{3}}$ $\displaystyle=$
$\displaystyle-\mu_{1}\mu_{2}e^{-\mu_{3}t}\cdot\left[\frac{e^{-(\mu_{1}-\mu_{3})t}+(\mu_{1}-\mu_{3})t-1}{(\mu_{3}-\mu_{1})^{2}(\mu_{2}-\mu_{1})}\right.$
$\displaystyle+$
$\displaystyle\left.\frac{e^{-(\mu_{2}-\mu_{3})t}+(\mu_{2}-\mu_{3})t-1}{(\mu_{3}-\mu_{2})^{2}(\mu_{1}-\mu_{2})}\right]$
$\displaystyle\frac{\partial P_{3}(\mu_{1},...,\mu_{4},t)}{\partial\mu_{4}}$
$\displaystyle=$
$\displaystyle-\left(\mathop{\displaystyle\prod}\limits_{i=1}^{3}\mu_{i}\right)e^{-\mu_{4}t}\left[\frac{e^{-(\mu_{1}-\mu_{4})t}+(\mu_{1}-\mu_{4})t-1}{(\mu_{4}-\mu_{1})^{2}(\mu_{3}-\mu_{1})(\mu_{2}-\mu_{1})}\right.$
$\displaystyle+$
$\displaystyle\left.\frac{e^{-(\mu_{2}-\mu_{4})t}+(\mu_{2}-\mu_{4})t-1}{(\mu_{4}-\mu_{2})^{2}(\mu_{3}-\mu_{2})(\mu_{1}-\mu_{2})}+\right.$
$\displaystyle+$
$\displaystyle\left.\frac{e^{-(\mu_{3}-\mu_{4})t}+(\mu_{3}-\mu_{4})t-1}{(\mu_{4}-\mu_{3})^{2}(\mu_{2}-\mu_{3})(\mu_{1}-\mu_{3})}\right]$
$\displaystyle\frac{\partial
P_{n}(\mu_{1},\ldots,\mu_{n+1},t)}{\partial\mu_{n+1}}$ $\displaystyle=$
$\displaystyle-\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot
e^{-\mu_{n+1}t}\sum\limits_{i=1}^{n}\frac{e^{-(\mu_{i}-\mu_{n+1})t}+(\mu_{i}-\mu_{n+1})t-1}{(\mu_{n+1}-\mu_{i})^{2}\cdot\mathop{\displaystyle\prod}\limits_{j=1,j\neq
i}^{n}(\mu_{j}-\mu_{i})}$ (63)

It is straightforward to show that $\partial P_{n}/\partial\mu_{n+1}<0$.
Indeed expanding $e^{-(\mu_{i}-\mu_{n+1})t}$ we obtain:

$\displaystyle\frac{\partial
P_{n}(\mu_{1},\ldots,\mu_{n+1},t)}{\partial\mu_{n+1}}$ $\displaystyle=$
$\displaystyle-\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot
e^{-\mu_{n+1}t}\sum\limits_{i=1}^{n}\frac{\sum\limits_{k=2}^{\infty}(-1)^{k}(\mu_{i}-\mu_{n+1})^{k}\frac{\cdot
t^{k}}{k!}}{(\mu_{n+1}-\mu_{i})^{2}\cdot\mathop{\displaystyle\prod}\limits_{j=1,j\neq
i}^{n}(\mu_{j}-\mu_{i})}$ $\displaystyle=$
$\displaystyle-\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot
e^{-\mu_{n+1}t}\sum\limits_{k=0}^{\infty}\frac{t^{k+2}}{(k+2)!}\sum\limits_{i=1}^{n}\frac{(-1)^{k}(\mu_{i}-\mu_{n+1})^{k}}{\mathop{\displaystyle\prod}\limits_{j=1,j\neq
i}^{n}(\mu_{j}-\mu_{i})}$

Then using (Eq. 35)

$\displaystyle\frac{\partial
P_{n}(\mu_{1},\ldots,\mu_{n+1},t)}{\partial\mu_{n+1}}$ $\displaystyle=$
$\displaystyle-\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot
e^{-\mu_{n+1}t}\sum\limits_{k=0}^{\infty}\frac{t^{k+2}}{(k+2)!}(-1)^{k-n+1}\left(\sum\limits_{i=1}^{n}\mu_{i}-\mu_{n+1}\right)^{\otimes{k-n+1}}$
$\displaystyle=$
$\displaystyle-\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\cdot
e^{-\mu_{n+1}t}\sum\limits_{r=0}^{\infty}\frac{t^{n+r+1}}{(n+r+1)!}(-1)^{r}\left(\sum\limits_{i=1}^{n}\mu_{i}-\mu_{n+1}\right)^{\otimes{r}}$

The last expression is similar to the Taylor series expansion of
$P_{n+1}(\mu_{i}-\mu_{n+1},t)$ (Eq. LABEL:eq:Pn_2) which is always positive.
Hence $\partial P_{n}/\partial\mu_{n+1}<0$

###  3.10 Analysis of functions $P_{n}(t)$

Expression for $P_{n}(t)$ where $(n\geq 1)$ (Eq. 15) is the sum of
monotonically decreasing functions or monotonically increasing functions,
depending on the sign of coefficeint before exponent. It is clear that
$P_{n}(t)$ is unique, smooth, differentiable function with no discontinuities,
which is equal to 0 when $t=0$ and when $t\rightarrow\infty$. From the Taylor
series expansion (Eq. LABEL:eq:Pn_2) the function $P_{n}(t)$ obeys the
following inequality:

$\displaystyle\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot\frac{t^{n}}{n!}\cdot\exp{\left[-\sum\limits_{i=1}^{n+1}\mu_{i}\cdot
t\right]}<P_{n}(t)<\left(\mathop{\displaystyle\prod}\limits_{i=1}^{n}\mu_{i}\right)\cdot\frac{t^{n}}{n!}\cdot
e^{-\mu_{\min}\cdot t}$ (66)

where $\mu_{\min}=\min(\mu_{1},\mu_{2},\ldots,\mu_{n+1}$). And it quite clear
that $P_{n}(t)$ is positive, bounded function.

##  4 Conclusion

So in this paper the generalization of the Poisson distribution is found for
the case of changing rates caused by consecutive events. In case of constant
event rate the distribution naturally transforms in the classical Poisson
distribution. The derived generalization can have different applications,
especially in simulation of cascade processes. It is possible that this
problem was already solved and published, but the author did not manage to
find such a publication. The author is grateful to S.M.Sergeev for the usefull
discussion.

## References

  * [1] E.Kushnirenko, ”Generalization of Poisson Distribution for the Case of Changing Probabilities of Consequential Events”, IX Workshop on High Energy Physics and Quantum Field Theory, NPI MSU, Zvenigorod, September 16-22, 1994. Workshop materials were published in Moscow State University, Moscow, 1995, pp.362-368. http://inspirehep.net/record/387538
  * [2] E.Kushnirenko, ”Estimates of the Electron-Gamma Conversion Background using Generalized Poisson Distribution”, VII International Workshop on Linear Colliders, Sept.29-Oct.3, 1997, Zvenigorod, Russia, Volume III, pp.1473-1476.
  * [3] Arsenin, V. Ya., Methods of Mathematical Physics, and Special Functions, ”Nauka”, Moscow, 1974. (Russian)
  * [4] G.Korn, T.Korn., Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review, 1970 (Russian)
  * [5] I.N. Bronshtein and K.A. Semendyayev. A Guide-Book to Mathematics for scientist and engineers, ”Nauka”, Moscow 1981. (Russian)
  * [6] Prudnikov, A.P. and Brychkov, I.U.A. and Marichev, O.I., Integrals and Series: Special functions, ”Nauka”, Moscow 1981. (Russian)


# Multiple UAV-Assisted Cooperative DF Relaying in Multi-User Massive MIMO IoT
Systems

Mobeen Mahmood, Yicheng Yuan, Tho Le-Ngoc
Department of Electrical and Computer Engineering, McGill University,
Montreal, QC, Canada
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>tho.le-
<EMAIL_ADDRESS>
###### Abstract

This work considers a multi-user massive multiple-input multiple-output (MU-
mMIMO) Internet-of-Things (IoT) system, where multiple unmanned aerial
vehicles (UAVs) operating as decode-and-forward (DF) relays connect the base
station (BS) to a large number of IoT devices. To maximize the total
achievable rate, we propose a novel joint optimization problem of hybrid
beamforming (HBF), multiple UAV relay positioning, and power allocation (PA)
to multiple IoT users. The study adopts a geometry-based millimeter-wave
(mmWave) channel model for both links and utilizes sequential optimization
based on K-means UAV-user association. The radio frequency (RF) stages are
designed based on the slow time-varying angular information, while the
baseband (BB) stages are designed utilizing the reduced-dimension effective
channel matrices. The illustrative results show that multiple UAV-assisted
cooperative relaying systems outperform a single UAV system in practical user
distributions. Moreover, compared to fixed positions and equal PA of UAVs and
BS, the joint optimization of UAV location and PA substantially enhances the
total achievable rate.

††This work was supported in part by Huawei Technologies Canada and in part by
the Natural Sciences and Engineering Research Council of Canada.

##  I Introduction

Relaying has long been considered as a potent means to mitigate signal
attenuation and fading effects, thereby extending the coverage and reliability
of wireless systems. The integration of relaying techniques in multi-user
massive multiple-input multiple-output (MU-mMIMO) Internet-of-Things (IoT)
systems is a practical consideration to further optimize spatial diversity and
multiplexing gains, thereby elevating the performance of wireless networks
[1]. In this regard, unmanned aerial vehicles (UAVs) as a moving relay can
offer the following advantages: 1) high mobility, which allows for the dynamic
adjustment of their locations to optimize communication conditions [2]; 2)
ease of deployment and low energy consumption; and 3) capability to reach
inaccessible locations. Thus, UAV-aided communications is a viable option for
future MU-mMIMO IoT applications, addressing coverage issues and minimizing
communications system overhead [3].

The high altitude deployment of UAVs increases the likelihood of line-of-sight
(LoS) dominated air-to-ground communication channels. While placing UAVs at
optimal locations is essential to improve channel quality, efficient
interference mitigation schemes are also required to enhance the network
capacity. The three-dimensional (3D) beamforming capability of MU-mMIMO can
effectively suppress interference among IoT nodes [4]. In this regard, two
common approaches for alleviating interference are fully-digital beamforming
(FDBF) and hybrid beamforming (HBF). FDBF requires radio frequency (RF) chains
equal to the number of antenna elements, which hinders its implementation in
UAV-assisted IoT systems due to the prohibitive cost, complexity, and limited
power supply of UAVs. Conversely, HBF involves the design of both the RF-stage
and baseband (BB)-stage, and can approach the performance of FDBF by reducing
the number of power-hungry RF chains, thereby improving energy efficiency [5,
6, 7]. Recent research studies have focused on optimizing UAV location with
particular emphasis on HBF solutions to maximize throughput or minimize
transmit power [8, 9, 10, 11]. In particular, [8, 9] investigate the joint
optimization of UAV deployment, while considering HBF at BS and UAV for
maximum sum-rate. An amply-and-forward UAV relay with analog beamforming
architecture is considered in [10] to maximize the capacity in a dual-hop
mMIMO IoT system. Similarly, [11] optimizes the source/relay power allocation
and UAV trajectory to enhance end-to-end system throughput.

Most existing studies consider a single UAV (for instance [8, 9, 10, 11]),
which can provide only limited user coverage and access. On the other hand, a
network of multiple UAVs can efficiently enlarge the coverage region and
increase the number of served users[12, 13, 14]. The authors in [12] propose a
multi-UAV relaying system, and compared the performance of a single multi-hop
link and multiple dual-hop links. In [13], the position of multiple UAV-
mounted BS is optimized to enhance the coverage area, while satisfying
quality-of-service (QoS) requirements. Similarly, the placement and the number
of UAVs are optimized in [14], while adhering to network capacity and coverage
constraints. Research works such as [8, 9, 10, 11, 12, 13, 14] have ignored
the direct link between the BS and IoT users. However, it is well known that
when direct link is non-negligible or not too weak, the spatial diversity can
be enhanced via direct and cooperative multiple path (through UAV relays). In
practice, assuming no direct link simplifies the design of the joint source-
relay beamforming. When the direct link is involved, the received signals via
source-to-destination link and multiple UAV relay-to-destination links are
combined at the destination (IoT user) to enhance the overall signal strength.
As such, the beamforming vector design at the BS needs to take into account
both of these links. Furthermore, the optimization problem becomes more
complex.

To solve this challenging issue, we consider a more practical cooperative
transmission approach, integrating both the direct link from BS to IoT devices
and the indirect links via multiple UAV relays to IoT users in MU-mMIMO IoT
systems to enhance the capacity and overcome the coverage issues. Our
objectives here are twofold: first, to show that the use of multiple UAVs as
relays can significantly increase the sum-rate capacity; and second, the
optimization of UAV location, user association, and power allocation (PA) at
BS and UAVs combined with HBF design can provide better performance than fixed
UAV locations with equal PA. To the best of our knowledge, the joint
optimization of multiple UAV placements, PA at the BS, and the design of HBF
for both the BS and UAVs is an unaddressed problem, presenting a significant
opportunity to advance the field of UAV-assisted MU-mMIMO IoT communications
networks. The joint optimization problem is highly non-convex, therefore, we
utilize structured sequential optimization to address the multi-faceted
optimization problem by splitting it into two subproblems. First, K-means-
based user clustering is used for UAV-users association based on the 3D
geometry-based millimeter-wave (mmWave) channel model. Then, the location of
each UAV is optimized jointly with PA using swarm intelligence. The RF
beamforming stages for BS and UAVs are designed based on the slow time-varying
angle-of-departure (AoD)/angle-of-arrival (AoA) information, and BB stages are
formulated using the reduced-dimensional effective channel matrices. The
illustrative results show that the proposed joint optimization scheme for
multiple UAV-assisted relaying significantly enhances the performance in MU-
mMIMO IoT systems.

##  II System and Channel Model

In this section, we introduce the system and channel models of the proposed
multi-UAV-assisted relaying and the HBF design for a multiple dual-hop MU-
mMIMO IoT system.

###  II-A System Model

We consider a downlink MU-mMIMO IoT network, where a large number of non-
overlapping IoT devices are connected to an IoT gateway via wire or wireless
links. Due to severe shadowing and blocking effect, many ground IoT devices
experience low signal-quality from the BS/eNodeB that is equipped with a large
array having $N_{b}$ antenna elements. To address this challenging scenario,
we consider a cooperative relaying system model as shown in Fig. 1(a), to
serve $\mathbb{K}=\\{1,\cdots,K\\}$ single-antenna IoT nodes, which are
clustered in $G$ groups, where $g^{th}$ group has $K_{g}$ IoT nodes such that
$K=\sum_{g=1}^{G}K_{g}$. Then, $M$ different UAVs, indexed by the set
$\mathbb{U}=\\{u_{1},u_{2},\cdots,u_{M}\\}$, are deployed to serve $K=MK_{m}$
IoT users, where $K_{m}$ is number of users served by $m^{th}$ UAV, which
operates as DF relay between BS/eNodeB and IoT node111We assume equal UAV-user
clustering for simplicity. However, it can be applicable for unequal user
clustering, which is left as our future work.. Let
$\left(x_{b},y_{b},z_{b}\right)$, $(x_{u}^{(m)},y_{u}^{(m)},z_{u}^{(m)})$ and
$\left(x_{k},y_{k},z_{k}\right)$ denote the locations of BS, $m^{th}$ UAV
relay, and $k^{th}$ IoT user, respectively. Then, we define the 3D distances
for multiple UAV-assisted MU-mMIMO IoT system as follows: $\tau_{1}^{(m)}$ =
$\sqrt{(x_{u}^{(m)}-x_{b})^{2}+(y_{u}^{(m)}-y_{b})^{2}+(z_{u}^{(m)}-z_{b})^{2}}$,
$\tau_{2,k}^{(m)}$ =
$\sqrt{(x_{u}^{(m)}-x_{k})^{2}+(y_{u}^{(m)}-y_{k})^{2}+(z_{u}^{(m)}-z_{k})^{2}}$,
$\tau_{k}$ = $\sqrt{(x_{b}-x_{k})^{2}+(y_{b}-y_{k})^{2}+(z_{b}-z_{k})^{2}}$,
where $\tau_{1}^{(m)}$, $\tau_{2_{k}}^{(m)}$ and $\tau_{k}$ are the 3D
distance between $m^{th}$ UAV & BS, between $m^{th}$ UAV and $k^{th}$ IoT
node, and between BS and $k^{th}$ IoT node, respectively. Each UAV is equipped
with $N_{r}(N_{t})$ antennas for receiving (transmitting) signals from (to) BS
(IoT users). For simplicity, we assume a homogeneous fleet of UAVs with
consistent specifications and functionality. Unlike traditional static
relaying, which uses fixed relay locations, we presumptively use multiple UAVs
as a movable relays.

(a)

(b)

Figure 1: Multiple UAV-assisted MU-mMIMO IoT communications. (a) network
model. (b) UAV as DF relay transmission phases. Figure 2: Multiple UAV-
assisted MU-mMIMO HBF system model.

For the downlink transmission of $N_{S}=K$ data streams, we consider HBF for
BS and all UAV relays as shown in Fig. 2. The BS consists of RF beamforming
stage $\mathbf{F}_{b}$ $\in$ $\mathbb{C}^{N_{b}\times N_{{RF}_{b}}}$, BB stage
$\mathbf{B}_{b}$ $\in$ $\mathbb{C}^{N_{{RF}_{b}}\times K}$, and MU PA matrix
$\mathbf{P}_{b}=\text{diag}(\sqrt{p_{b_{1}}},\cdots,\sqrt{p_{b_{K}}})$ $\in$
$\mathbb{C}^{K\times K}$. Here, $N_{{RF}_{b}}$ is the RF chains such that
$N_{S}\leq N_{{RF}_{b}}\leq N_{b}$ to guarantee multi-stream transmission. We
consider half-duplex (HD) DF relaying for each UAV. Therefore, each round of
information transmission from BS to IoT nodes can be divided into two phases:
1) source phase (SP); and 2) relay phase (RP) as illustrated in Fig. 1(b). In
SP, BS transmits $K$ data streams to the following: 1) $K$ IoT users through
channel $\mathbf{H}_{D}$ $\in$ $\mathbb{C}^{K\times N_{b}}$; and 2) each
$m^{th}$ UAV via channel $\mathbf{H}_{1}^{(m)}$ $\in$ $\mathbb{C}^{N_{r}\times
N_{b}}$. Using $N_{r}$ antennas, each UAV receives signals with RF stage
$\mathbf{F}_{u,r}^{(m)}$ $\in$ $\mathbb{C}^{N_{{RF}_{u}}\times N_{r}}$ and BB
stage $\mathbf{B}_{u,r}^{(m)}$ $\in$ $\mathbb{C}^{K\times N_{{RF}_{u}}}$. In
RP, each UAV decodes the received information and then forwards the decoded
information to $K$ IoT nodes using RF beamformer $\mathbf{F}_{u,t}^{(m)}\in$
$\mathbb{C}^{N_{t}\times N_{{RF}_{u}}}$, BB stage $\mathbf{B}_{u,t}^{(m)}\in$
$\mathbb{C}^{N_{{RF}_{u}}\times K}$, and MU PA matrix
$\mathbf{P}_{u}^{(m)}=\text{diag}(\sqrt{p_{u_{1}}^{(m)}},\cdots,\sqrt{p_{u_{K}}^{(m)}})$
$\in$ $\mathbb{C}^{K\times K}$ through channel $\mathbf{H}_{2}^{(m)}$ $\in$
$\mathbb{C}^{K\times N_{t}}$, where $p_{u_{k}}^{(m)}$ reflects the allocated
power to $k^{th}$ user from $m^{th}$ UAV. The implementation of all RF
beamforming/combining stages involves the use of phase-shifters (PSs) and
thus, impose a constant-modulus (CM) constraint, i.e.,
$|\mathbf{F}_{b}(i,j)|=\frac{1}{\sqrt{N_{b}}}$,
$|\mathbf{F}_{u,r}^{(m)}(i,j)|=\frac{1}{\sqrt{N_{r}}}$,
$|\mathbf{F}_{u,t}^{(m)}(i,j)|=\frac{1}{\sqrt{N_{t}}}\hskip 4.30554pt\forall
i,j,m$. For the data signal $\mathbf{d}=[d_{1},d_{2},\cdots,d_{K}]^{T}$ with
$\operatorname{\mathbb{E}}\\{\mathbf{d}\mathbf{d}^{H}\\}=\mathbf{I}_{K}$ $\in$
$\mathbb{C}^{K\times K}$, the signal received at $m^{th}$ UAV (after BB
processing) during SP is given as follows:

$\tilde{\mathbf{y}}_{u}^{(m)}=\mathbf{B}_{u,r}^{(m)}\mathbf{F}_{u,r}^{(m)}\mathbf{H}_{1}^{(m)}\mathbf{F}_{b}\mathbf{B}_{b}\mathbf{P}_{b}\mathbf{d}+\mathbf{B}_{u,r}^{(m)}\mathbf{F}_{u,r}^{(m)}\mathbf{n}_{u}^{(m)},\vspace{-0.5ex}$
(1)

where $n_{u}^{(m)}$ $\in$ $\mathbb{C}^{N_{u,r}}$ denotes the zero-mean complex
circularly symmetric Gaussian noise vector at $m^{th}$ UAV relay with
covariance matrix
$\operatorname{\mathbb{E}}\\{\mathbf{n}_{u}\mathbf{n}_{u}^{H}\\}=\sigma_{n_{u}}^{2}\mathbf{I}_{N_{u,r}}$
$\in$ $\mathbb{C}^{N_{u,r}\times N_{u,r}}$. Then, the signal transmitted by
$m^{th}$ UAV during RP is given as:

$\hat{\mathbf{s}}^{(m)}=\mathbf{F}_{u,t}^{(m)}\mathbf{B}_{u,t}^{(m)}\mathbf{P}_{u}^{(m)}\hat{\mathbf{d}}^{(m)},\vspace{-1ex}$
(2)

where $\mathbf{\hat{d}}^{(m)}$ is the re-encoded signal at $m^{th}$ UAV relay.
Each IoT node receives signal from BS and each $m^{th}$ UAV during SP and RP,
respectively. Then, the received signal at $k^{th}$ IoT node from BS and
$m^{th}$ UAV can be written as:

$\begin{split}y_{k}^{(m)}&=y_{k,SP}+y_{k,RP}^{(m)},\\\
&=\mathbf{h}_{D,k}^{T}\mathbf{F}_{b}\mathbf{B}_{b}\mathbf{P}_{b}\hat{d}_{k}+n_{D_{k}}+\mathbf{h}_{2,k}^{T(m)}\mathbf{F}_{u,t}^{(m)}\mathbf{B}_{u,t}^{(m)}\mathbf{P}_{u}^{(m)}\hat{d}_{k}^{(m)}+n_{2_{k}}^{(m)},\\\
\vspace{-4ex}&=\mbox{$\underbrace{\sqrt{p_{b_{k}}}\mathbf{h}_{D,k}^{T}\mathbf{F}_{b}\mathbf{B}_{b}\mathbf{P}_{b}\hat{d}_{k}+\sqrt{p_{u_{k}^{(m)}}}\mathbf{h}_{2,k}^{T(m)}\mathbf{F}_{u,t}^{(m)}\mathbf{B}_{u,t}^{(m)}\mathbf{P}_{u}^{(m)}\hat{d}_{k}^{(m)}}_{\text{Desired
Signal}}+$}\\\
&\mbox{$\underbrace{\sum_{\hat{k}=1}^{K}\sqrt{p_{b_{{\hat{k}}}}}\mathbf{h}_{D,k}^{T}\mathbf{F}_{b}\mathbf{B}_{b}\mathbf{P}_{b}\hat{d}_{{\hat{k}}}+\sum_{\hat{k}=1}^{K}\sqrt{p_{u_{k}^{(m)}}}\mathbf{h}_{2,k}^{T(m)}\mathbf{F}_{u,t}^{(m)}\mathbf{B}_{u,t}^{(m)}\mathbf{P}_{u}^{(m)}\hat{d}_{{\hat{k}}}^{(m)}}_{\text{Total
MU-interference}}$}\vspace{-1ex}\\\
&+\underbrace{n_{D_{k}}+n_{2_{k}}^{(m)}}_{\text{Total Noise}},\end{split}$ (3)

where $n_{D_{k}}$ $\sim$ $\mathcal{CN}(\mathbf{0},\sigma_{n_{D}}^{2})$ and
$n_{2_{k}}^{(m)}$ $\sim$ $\mathcal{CN}(\mathbf{0},\sigma_{n_{2}}^{2})$ are the
additive circular symmetric Gaussian noise at $k^{th}$ IoT node. The power
constraint of the beamforming matrices can be expressed as
$\lVert\mathbf{F}_{b}\mathbf{B}_{b}\mathbf{P}_{b}\rVert^{2}_{F}$ = $P_{T}$ and
$\lVert\mathbf{F}_{u,t}^{(m)}\mathbf{B}_{u,t}^{(m)}\mathbf{P}_{u}^{(m)}\rVert^{2}_{F}$
= $P_{u}^{(m)}$, where $P_{T}$ and $P_{u}^{(m)}$ denote the total transmit
power of BS and $m^{th}$ UAV, respectively. Then, the achievable rate of first
link (i.e., BS $\rightarrow$ $m^{th}$ UAV) is given as follows:

$\mbox{$\mathrm{R}_{1}^{(m)}=\log_{2}\left|\mathbf{I}_{K}+\mathbf{Q}_{1}^{\mathrm{-1}(m)}\mathbf{B}_{u,r}^{(m)}\mathbfcal{H}_{1}^{(m)}\mathbf{B}_{b}\mathbf{B}_{b}^{H}\mathbfcal{H}_{1}^{H(m)}\mathbf{B}_{u,r}^{H(m)}\right|$},$
(4)

where
$\mathbf{Q}_{1}^{\mathrm{-1}(m)}=(\sigma_{n_{u}}^{2}\mathbf{B}_{u,r}^{(m)}\mathbf{F}_{u,r}^{(m)})^{\mathrm{-1}}\mathbf{F}_{u,r}^{H(m)}\mathbf{B}_{u,r}^{H(m)}$
and
$\mathbfcal{H}_{1}^{(m)}=\mathbf{F}_{u,r}^{(m)}\mathbf{H}_{1}^{(m)}\mathbf{F}_{b}$.
The signal-to-interference-plus-noise ratio (SINR) of $k^{th}$ IoT node via
$m^{th}$ UAV is given as [15]:

$\gamma_{k}^{(m)}=\frac{p_{b_{k}}|\mathbf{h}_{D_{k}}^{H}\mathbf{F}_{b}\mathbf{b}_{b_{k}}|^{2}}{\sum\limits_{\hat{k}\neq
k}p_{b_{{\hat{k}}}}|\mathbf{h}_{D_{k}}^{H}\mathbf{F}_{b}\mathbf{b}_{b_{{\hat{k}}}}|^{2}+\sigma_{n_{D}}^{2}}+\frac{p_{u_{k}}^{(m)}|\mathbf{h}_{2_{k}}^{H(m)}\mathbf{F}_{u,t}^{(m)}\mathbf{b}_{u,t_{k}}^{(m)}|^{2}}{\sum\limits_{\hat{k}\neq
k}p_{u_{{\hat{k}}}}^{(m)}|\mathbf{h}_{2_{k}}^{H(m)}\mathbf{F}_{u,t}^{(m)}\mathbf{b}_{u,t_{{\hat{k}}}}^{(m)}|^{2}+\sigma_{n_{2}}^{2(m)}}.$
(5)

The achievable rate for the $m^{th}$ second link (i.e., from $m^{th}$ UAV to
$K_{m}$ users) can be written as:

$\mathrm{R}_{2}^{(m)}=\operatorname{\mathbb{E}}\bigl{\\{}\sum_{i=1}^{K_{m}}\log_{2}(1+\gamma_{i}^{(m)})\bigr{\\}}.$
(6)

For multiple dual-hop MU-mMIMO IoT system, where each $m^{th}$ UAV is deployed
at a fixed height $z_{u}^{m}$, and relaying data to $K_{m}$ IoT nodes, the
total achievable rate can be maximized by the joint optimization of
$\mathbf{F}_{b}$, $\mathbf{B}_{b}$, $\mathbf{F}_{u,t}^{(m)}$,
$\mathbf{F}_{u,r}^{(m)}$, $\mathbf{B}_{u,t}^{(m)}$, $\mathbf{B}_{u,r}^{(m)}$,
$\mathbf{P}_{b}$, $\mathbf{P}_{u}^{(m)}$ and the UAV location
$\mathbf{x}^{(m)}=[x_{o}^{(m)},y_{o}^{(m)}]^{T}\in\mathbb{R}^{2}$, which is to
be optimized within the given deployment area. Then, we can formulate the
optimization problem as:

$\begin{split}&\max_{\left\\{\mathbf{F}_{b},\mathbf{B}_{b},\mathbf{F}_{u,t}^{(m)},\mathbf{B}_{u,t}^{(m)},\mathbf{F}_{u,r}^{(m)},\mathbf{B}_{u,r}^{(m)},\mathbf{P}_{b},\mathbf{P}_{u}^{(m)},\mathbf{x}^{(m)}\right\\}}\quad\mathrm{R}_{T}\\\
&\hskip 17.22217pt\textrm{s.t.}\hskip
6.99997ptC_{1}:|\mathbf{F}_{u,t}^{(m)}(i,j)|=\frac{1}{\sqrt{N_{t}}},|\mathbf{F}_{u,r}^{(m)}(i,j)|=\frac{1}{\sqrt{N_{r}}},\hskip
2.15277pt\forall i,j,m,\\\ &\quad\hskip
25.83325ptC_{2}:|\mathbf{F}_{b}(i,j)|=\frac{1}{\sqrt{N_{b}}},\forall i,j,\\\
&\quad\hskip 25.83325ptC_{3}:\bigcup_{m\in M}K_{m}=K,\hskip 8.61108pt\forall
m,\\\ &\quad\hskip
25.83325ptC_{4}:\operatorname{\mathbb{E}}\\{\left\lVert\mathbf{s}_{1}\right\rVert_{2}^{2}\\}\leq
P_{T},\operatorname{\mathbb{E}}\\{||\hat{\mathbf{s}}^{(m)}||_{2}^{2}\\}\leq
P_{T},\forall m,\\\ &\quad\hskip 25.83325ptC_{5}:p_{b_{k}}\geq
0,p_{u_{k}}^{(m)}\geq 0,\hskip 8.61108pt\forall k,m,\\\ &\quad\hskip
25.83325ptC_{6}:\mathbf{x}_{\mathrm{min}}\leq\mathbf{x}_{o}^{(m)}\leq\mathbf{x}_{\mathrm{max}},\hskip
8.61108pt\forall m,\end{split}$ (7)

where
$\mathrm{R}_{T}=\sum\nolimits_{m=1}^{M}(1/2)\min(\mathrm{R}_{1}^{(m)},\mathrm{R}_{2}^{(m)})$
is the total achievable rate, $C_{1}$ and $C_{2}$ refers to the CM constraint
due to the use of PSs for UAV and BS, respectively, $C_{3}$ ensures that the
total user count is consistent with the counts for each UAV, $C_{4}$ indicates
the total transmit power constraint for UAV and BS, $C_{5}$ represents the
non-negative allocated power to each IoT node from BS and each UAV, and
$C_{6}$ implies UAV deployment within the given flying span. Here,
$[\mathbf{x}_{\mathrm{min}},\mathbf{x}_{\mathrm{max}}]=[(x_{\mathrm{min}},y_{\mathrm{min}}),(x_{\mathrm{max}},y_{\mathrm{max}})]$
represents the deployment range for each UAV in 2-D space. The optimization
problem defined in (7) is non-convex and intractable. Thus, we develop a sub-
optimal solution for (7) in Section III.

###  II-B Channel Model

We consider mmWave channels for both links. The channel between BS and
$m^{th}$ UAV is modeled based on the Saleh-Valenzuela channel model, and is
given as:

$\mathbf{H}_{1}^{(m)}=\sum\limits_{c=1}^{C}\sum\limits_{l=1}^{L}z_{1_{cl}}^{(m)}\tau_{1_{cl}}^{-\eta(m)}\mathbf{a}_{1,r}^{(m)}({{\theta^{(m)}_{r_{cl}}}},{{\phi^{(m)}_{r_{cl}}}})\mathbf{a}_{1,t}^{T(m)}({{\theta_{t_{cl}}^{(m)}}},{{\phi_{t_{cl}}^{(m)}}}),$
(8)

where $C$ is the total number of clusters, $L$ is the total number of paths,
$\eta$ is the path loss exponent, $z_{1_{cl}}^{(m)}$ is the complex gain of
$l^{th}$ path in $c^{th}$ cluster for $m^{th}$ UAV, and
${\bf{a}}_{1,j}^{(m)}(\cdot,\cdot)$ is the corresponding transmit or receive
array steering vector for uniform rectangular array (URA), which is given as
[16]:

$\begin{split}{\bf{a}}_{1,j}^{(m)}({{\theta^{(m)},\phi^{(m)}}})&=\big{[}{1,\cdots,{e^{-j2\pi
d\left({N_{x}-1}\right){{\sin(\theta^{(m)})\cos(\phi^{(m)})}}}}}\big{]}\\\
&\otimes\big{[}{1,\cdots,{e^{-j2\pi
d\left({N_{y}-1}\right){{\sin(\theta^{(m)})\sin(\phi^{(m)})}}}}}\big{]},\end{split}$
(9)

where $j=\\{t,r\\}$, $d$ is the inter-element spacing, and $N_{x}(N_{y})$ is
the horizontal (vertical) size of corresponding antenna array at BS and UAV.
Here, the angles
${\theta_{j_{cl}}^{(m)}}\in\big{[}{{\theta}_{j_{c}}^{(m)}-{\delta_{j_{c}}^{{}^{\theta(m)}}}},{{\theta}_{j_{c}}^{(m)}+{\delta_{j_{c}}^{{}^{\theta(m)}}}}\big{]}$
and
${\phi_{j_{cl}}^{(m)}}\in\big{[}{{\phi}_{j_{c}}^{(m)}-{\delta_{j_{c}}^{\phi(m)}}},{{\phi}_{j_{c}}^{(m)}+{\delta_{j_{c}}^{\phi(m)}}}\big{]}$
are the elevation and azimuth AoD($j=t$) or AoA ($j=r$) for $l^{th}$ path in
channel $\mathbf{H}_{1}^{(m)}$, respectively. Here, ${\theta}_{j_{c}}^{(m)}$
and ${\phi}_{j_{c}}^{(m)}$ are the mean elevation and azimuth angles,
respectively with ${\delta_{j_{c}}^{{}^{\theta(m)}}}$
(${\delta_{j_{c}}^{\phi(m)}}$) represents the elevation(azimuth) angle spread.
The channel vector between the UAV(or BS) and the $k^{th}$ IoT node can be
written as:

$\mathbf{h}_{i,k}^{T}=\sum\nolimits_{q=1}^{Q}z_{i,k_{q}}\tau_{i,k_{q}}^{-\eta}\mathbf{a}(\theta_{k_{q}},\phi_{k_{q}})=\mathbf{z}_{i,k}^{T}\mathbf{A}_{i,k}\in\mathbb{C}^{N},$
(10)

where $i=\\{D,2\\}$, $Q$ is the total number of downlink paths from UAV(or BS)
to IoT nodes, $z_{i,k_{q}}$ $\sim$ $\mathcal{CN}(0,\frac{1}{Q})$ is the
complex path gain of $q^{th}$ path, $\mathbf{a}(\cdot,\cdot)$ $\in$
$\mathbb{C}^{N}$ is the downlink array phase response vector. Then, the
complete channel matrix for $K$ IoT nodes can be written as:

$\mathbf{H}_{i}=[\mathbf{h}_{i,1},\cdots,\mathbf{h}_{i,K}]^{T}=\mathbf{Z}_{i}\mathbf{A}_{i}\in\mathbb{C}^{K\times
N_{t}},$ (11)

where
$\mathbf{Z}_{i}=[\mathbf{z}_{i,1},\cdots,\mathbf{z}_{i,K}]^{T}\in\mathbb{C}^{K\times
Q}$ is the complete path gain matrix for all downlink IoT nodes and
$\mathbf{A}_{i,k}\in\mathbb{C}^{Q\times N}$ is the slow time-varying array
phase response matrix.

##  III Proposed Joint User Association, Multiple UAV Positioning, PA & Hybrid
Beamforming

In this section, our objective is to optimize each UAV location jointly with
PA from BS and UAV, and sequentially design HBF stages for BS and each UAV to
reduce the channel state information (CSI) overhead size while maximizing the
throughput of a multiple UAV-assisted MU-mMIMO IoT system. First, we discuss
the UAV-user association using K-means-based user clustering.

###  III-A UAV-User Association

The proposed scheme leverages K-means-based user association to assign $K_{m}$
users to nearest $m^{th}$ UAV while maintaining exclusive user-UAV
associations. The objective is to minimize the sum of squared distances
between users and their assigned $m^{th}$ UAV (i.e.,
$\min\sum_{m=1}^{M}\sum_{k=1}^{K}z_{km}||\mathbf{x}_{k}-\mathbf{x}_{u}^{(m)}||^{2}$).
Here, $z_{km}$ is the assignment variable, which is defined as:

$z_{km}=\begin{cases}1,&\text{if user $k$ is assigned to UAV $m$},\\\
0,&\text{otherwise}.\end{cases}$ (12)

###  III-B Joint UAV Deployment, Optimal PA and HBF Design

In this section, our objective is to design the HBF stages for BS and $M$ UAVs
by using sequential optimization. Initially, both RF and BB stages are
constructed using some fixed UAV locations. Then, we employ swarm intelligence
to optimize each UAV location as well as PA from BS and each UAV for maximum
total achievable rate. Finally, the RF and BB stages are re-formulated for the
optimal UAV location as well as the allocated power in the MU PA blocks
$\mathbf{P}_{b}$ and $\mathbf{P}_{u}^{(m)}$.

####  III-B1 RF & BB Stage Design

The RF beamforming stage for BS and each UAV (both transmit and receive) are
designed as:

$\vspace{-0.5ex}{{\bf{F}}}=\big{[}{{\bf{e}}_{j}\big{(}{\lambda^{u_{1}}_{x},\lambda^{k_{1}}_{y}}\big{)},\cdots,{\bf{e}}_{j}\big{(}{\lambda^{u_{N_{{RF}}}}_{x},\lambda^{k_{N_{{RF}}}}_{y}}\big{)}}\big{]}\in\mathbb{C}^{N_{T}\times
N_{{RF}}},$ (13)

where $j$ $=$ $\\{t,r\\}$ and $\mathbf{e}(\cdot,\cdot)$ is the corresponding
transmit or receive steering vector, which is defined as
${\bf{e}}\left({{\theta,\phi}}\right)=\frac{1}{\mathcal{N}_{t}}\big{[}{1,{e^{j2\pi
d{{\sin(\theta)\cos(\phi)}}}},\cdots,{e^{j2\pi
d\left({\mathcal{N}_{x,t}-1}\right){{\sin(\theta)\cos(\phi)}}}}}\big{]}^{T}\otimes\big{[}{1,{e^{j2\pi
d{{\sin(\theta)\sin(\phi)}}}},\cdots,{e^{j2\pi
d\left({\mathcal{N}_{y,t}-1}\right){{\sin(\theta)\sin(\phi)}}}}}\big{]}^{T}$,
where $\mathcal{N}_{T}=\\{N_{b},N_{t},N_{r}\\}$. Here, the RF beamformers are
constructed via quantized angle-pairs, which are defined as
${{\lambda_{x}^{u}}=-1+\frac{2u-1}{{{\mathcal{N}_{x,t}}}}}$ for
$u=1,\cdots,{\mathcal{N}_{x,t}}$ and
${{\lambda_{y}^{k}}=-1+\frac{2k-1}{{{\mathcal{N}_{y,t}}}}}$ for
$k=1,\cdots,{\mathcal{N}_{y,t}}$. The quantized angle-pairs reduces the number
of RF chains at BS and each UAV while providing complete AoD/AoA supports,
which are defined as:

$\displaystyle\textrm{AoD}=\big{\\{}{\sin\left(\theta\right){{\left[{\cos\left(\phi\right),\sin\left(\phi\right)}\right]}}}\big{|}{\theta}\in\bm{\theta}_{t},{\phi}\in\bm{\phi}_{t}\big{\\}},$
(14)
$\displaystyle\textrm{AoA}=\big{\\{}{\sin\left(\theta\right){{\left[{\cos\left(\phi\right),\sin\left(\phi\right)}\right]}}}\big{|}\theta\in{\bm{\theta}}_{r},\phi\in{\bm{\phi}}_{r}\big{\\}},$
(15)

where
$\bm{\theta}_{i}=\big{[}{{\theta_{i}}-{\delta_{i}^{\theta}}},{{\theta_{i}}+{\delta_{i}^{\theta}}}\big{]}$
and
$\bm{\phi}_{i}=\big{[}{{\phi_{i}}-{\delta_{i}^{\phi}}},{{\phi_{i}}+{\delta_{i}^{\phi}}}\big{]}$
denote the azimuth and elevation angle supports, respectively. After designing
the transmit and receive RF beamformers for BS, and each UAV, the effective
channel matrices $\mathbfcal{H}_{1}^{(m)}$ and $\mathbfcal{H}_{2}^{(m)}$ as
seen from the BB-stages are given as follows:

$\mathbfcal{H}_{1}^{(m)}=\mathbf{F}_{u,r}^{(m)}\mathbf{H}_{1}^{(m)}\mathbf{F}_{b}=\mathbf{U}_{1}^{(m)}\mathbf{\Sigma}_{1}^{(m)}\mathbf{V}_{1}^{H(m)},$
(16)
$\displaystyle\mathbfcal{H}_{2}^{(m)}=\left[\begin{array}[]{ccc}\mathbf{H}_{2,1}^{(m)}\mathbf{F}_{u,t,1}^{(m)}&\ldots&\mathbf{H}_{2,1}^{(m)}\mathbf{F}_{u,t,G}^{(m)}\\\
\vdots&\ddots&\vdots\\\
\mathbf{H}_{2,G}^{(m)}\mathbf{F}_{u,t,1}^{(m)}&\ldots&\mathbf{H}_{2,G}^{(m)}\mathbf{F}_{u,t,G}^{(m)}\end{array}\right],\vspace{-0.5ex}$
(20)

where $\mathbf{U}_{i}^{(m)}$ $\in$ $\mathbb{C}^{N_{{RF}_{u}}\times
rank(\mathbfcal{H}_{1}^{(m)})}$ and $\mathbf{V}_{i}^{(m)}$ $\in$
$\mathbb{C}^{N_{{RF}_{b}}\times rank(\mathbfcal{H}_{1}^{(m)})}$ are tall
unitary matrices and $\mathbf{\Sigma}_{1}^{(m)}$ is the diagonal matrix with
singular values in the decreasing order for $m^{th}$ UAV. Then,
$\mathbf{B}_{u,r}^{(m)}$ for $m^{th}$ UAV is defined as:

$\mathbf{B}_{u,r}^{(m)}=\mathbf{U}_{1}^{H(m)}.$ (21)

The reduced-size effective CSI $\bm{\mathcal{H}}_{2}^{(m)}$ given in (20) is
employed for designing $\mathbf{B}_{u,t}^{(m)}$ by using regularized zero-
forcing (RZF) technique, and is defined as follows:

$\mathbf{B}_{u,t}^{(m)}=(\mathbfcal{H}_{2}^{H(m)}\mathbfcal{H}_{2}^{(m)}+\beta^{(m)}N_{{RF}_{u}}\mathbf{I}_{N_{{RF}_{u}}})^{-1}\mathbfcal{H}_{2}^{H(m)},\vspace{-1ex}$
(22)

where $\beta^{(m)}=\frac{\sigma^{2(m)}}{P_{T}^{(m)}}$ is the regularization
parameter and $\mathbf{I}_{N_{{RF}_{u}}}$ $\in$
$\mathbb{C}^{N_{{RF}_{u}}\times N_{{RF}_{u}}}$. The optimal design of
$\mathbf{B}_{b}$ is formulated using an effective channel
$\mathbfcal{H}_{1D}$, which constitutes all channel components from BS (i.e.,
$\mathbf{H}_{D}$ and $\mathbf{H}_{2}^{(m)},\forall m=1,\cdots,M$). Then, the
total effective channel can be written as:

$\mathbfcal{H}_{1D}=[\mathbf{H}_{D},\mathbf{H}_{1}^{(1)},\mathbf{H}_{1}^{(2)},\cdots,\mathbf{H}_{1}^{(M)}]\in\mathbb{C}^{(MN_{r}+K)\times
K}.$ (23)

By using SVD of effective channel $\mathbfcal{H}_{1D}$, we can design
$\mathbf{B}_{b}$ by using tall unitary matrix $\mathbf{V}_{1D}$ $\in$
$\mathbb{C}^{N_{{RF}_{b}}\times rank(\mathbfcal{H}_{1D})}$ as [16]:

$\mathbf{B}_{b}=\sqrt{\frac{P_{T}}{K}}\mathbf{V}_{1D}\in\mathbb{C}^{N_{{RF}_{b}}\times
K}.\vspace{-0.5ex}$ (24)

TABLE I: Simulation Parameters [17]

Number of antennas | $(N_{b},N_{r},N_{t})=64$
---|---
Number of paths | Path loss exponent | $L=10$ | $3.6$
Frequency | Channel Bandwidth | 28 GHz | 100 MHz
Noise PSD | Number of UAVs $M$ |  $-174$ dBm/Hz | 2 or 3
BS | UAV height | 10 m | 20 m
UAV x-axis range | UAV y-axis range | $[x_{min},x_{max}]=\left[0,100\right]m$ | $[y_{min},y_{max}]=\left[0,100\right]m$
Azimuth AoD/AoA ($1^{st}$ link) | Azimuth AoD/AoA ($2^{nd}$ link) | $120^{\circ}$ | $150^{\circ}$
Elevation AoD/AoA ($1^{st}$ link) | Elevation AoD/AoA ($2^{nd}$ link) | $60^{\circ}$ | $30^{\circ}$
Azimuth/Elevation Angle Spread | # of network realizations | $\pm 10^{\circ}$ | 2000

####  III-B2 Joint Multiple UAV Positioning and Optimal PA

After the design of RF and BB stages for BS and UAVs, the optimization problem
given in (7) can be reformulated as:

$\displaystyle\max_{\left\\{\mathbf{P}_{b},\mathbf{P}_{u}^{(m)},\mathbf{x}^{(m)}\right\\}}\quad\mathrm{R}_{T}$
(25) $\displaystyle\hskip 21.52771pt\textrm{s.t.}\hskip
21.52771ptC_{4}-C_{6}.$

This resulting problem in (25) is still non-convex due to the joint dependence
of both the allocated powers $p_{b_{k}}$, $p_{u_{k}}^{(m)}$ and the UAV
location $\mathbf{x}^{(m)}=[x_{o}^{(m)},y_{o}^{(m)}]^{T}$ on the SINR
expression in (5), which is used in the sum-rate $\mathrm{R}_{2}^{(m)}$
calculation as given in (6). To overcome this challenge, we propose sequential
optimization using swarm intelligence, which employs multiple agents, called
particles, to explore the search space of objective function given in (25).
Initially, $N_{p}$ particles are randomly placed in search space, where each
particle communicates with other particles to share their personal best and
update the current global best solution for the objective function. The
particles then move iteratively for $T$ iterations to reach the global optimum
solution. In particular, each UAV location
$\mathbf{x}_{u}^{(m)}=[x_{u}^{(m)},y_{u}^{(m)}]^{T}$ and $\mathbf{P}_{b}$,
$\mathbf{P}_{u}^{(m)}$ are optimized by using particle swarm optimization
(PSO)-based algorithmic solution while maximizing the total achievable rate.
Here, the $i^{th}$ particle at the $t^{th}$ iteration now represents an
instance of the each UAV location and multi-user PA matrices, which is given
as follows:

$\begin{split}&\mbox{$\mathbf{J}_{p_{i}}^{(t)}=[\mathbf{X}_{i}^{(t)},\hat{\mathbf{P}}_{b_{i}}^{(t)},\hat{\mathbf{P}}_{u_{i}}^{(t)}]^{T}=[x_{1_{i}}^{(t)},y_{1_{i}}^{(t)},\cdots,x_{M_{i}}^{(t)},y_{M_{i}}^{(t)},\sqrt{\hat{p}_{b,1_{i}}^{(t)}},\cdots,\sqrt{\hat{p}_{b,K_{i}}^{(t)}},\cdots,$}\\\
&\mbox{$\cdots,\sqrt{\hat{p}_{u_{1},1_{i}}^{(t)}},\cdots,\sqrt{\hat{p}_{u_{1},K_{i}}^{(t)}},\cdots,\sqrt{\hat{p}_{u_{M},1_{i}}^{(t)}},\cdots,\sqrt{\hat{p}_{u_{M},K_{i}}^{(t)}}]^{T}\in\mathbb{R}^{2K+2M}$},\end{split}$
(26)

where each particle $i$ represents the $M$ UAV positions and PA to $K$ IoT
users from BS and each $m^{th}$ UAV, and calculates the objective function as
$R_{T}(\mathbf{F}_{b},\mathbf{B}_{b},\mathbf{F}_{u,t}^{(m)},\mathbf{B}_{u,t}^{(m)},\mathbf{F}_{u,r}^{(m)},\mathbf{B}_{u,r}^{(m)},\kappa_{b_{i}}^{(t)}\hat{\mathbf{P}}_{b_{i}}^{(t)},\kappa_{u_{i}}^{(t)}\hat{\mathbf{P}}_{u_{i}}^{(t)},\mathbf{X}_{i}^{(t)})$.
We define
$\hat{\mathbf{P}}_{b_{i}}^{(t)}=\operatorname{diag}(\sqrt{\hat{p}_{b,1_{i}}^{(t)}},\cdots,\sqrt{\hat{p}_{b,K_{i}}^{(t)}})\in\mathbb{R}^{K\times
K}$ and
$\hat{\mathbf{P}}_{u_{i}}^{(t)}=[\hat{\mathbf{P}}_{u_{1_{i}}}^{(t)},\cdots,\hat{\mathbf{P}}_{u_{M_{i}}}^{(t)}]$,
where
$\hat{\mathbf{P}}_{u_{m_{i}}}^{(t)}=\operatorname{diag}(\sqrt{\hat{p}_{u_{m},1_{i}}^{(t)}},\cdots,\sqrt{\hat{p}_{u_{m},K_{i}}^{(t)}})\in\mathbb{R}^{K\times
K}$ as the normalized PA matrices with $\hat{p}_{b,k_{i}}^{(t)}$,
$\hat{p}_{u_{m},k_{i}}^{(t)}\in[0,1]$. Then, the transmit power constraints
for $\hat{\mathbf{P}}_{b_{i}}^{(t)}$ and $\hat{\mathbf{P}}_{u_{i}}^{(t)}$ are
satisfied by defining
$\mathbf{P}_{b_{i}}^{(t)}=\kappa_{b_{i}}^{(t)}\mathbf{\hat{P}}_{b_{i}}^{(t)}$
and
$\mathbf{P}_{u_{i}}^{(t)}=\kappa_{u_{i}}^{(t)}\mathbf{\hat{P}}_{u_{i}}^{(t)}$.
The position $\mathbf{J}_{p_{i}}^{(t)}$ and velocity
$\mathbf{J}_{v_{i}}^{(t)}$ for $i^{th}$ particle during $t^{th}$ iteration are
updated as follows:

$\mathbf{J}_{p_{i}}^{(t+1)}=\mathbf{J}_{p_{i}}^{(t)}+\mathbf{J}_{v_{i}}^{(t+1)},$
(27)
$\mbox{$\mathbf{J}_{v_{i}}^{(t+1)}=\gamma_{1}\mathbf{Y}_{1}^{(t)}(\mathbf{J}_{p_{\mathrm{best}}}^{(t)}-\mathbf{J}_{p_{i}}^{(t)})+\gamma_{2}\mathbf{Y}_{2}^{(t)}(\mathbf{J}_{p_{\mathrm{best}_{i}}}^{(t)}-\mathbf{J}_{p_{i}}^{(t)})+\gamma_{3}^{(t)}\mathbf{J}_{v_{i}}^{(t)}$}.\vspace{-1ex}$
(28)

Finally, the personal and global best solutions for $i^{th}$ particle during
$t^{th}$ iteration are obtained as follows:

$\begin{split}&\mathbf{J}_{p_{\mathrm{best}_{i}}}^{(t)}=\operatorname*{arg\,max}_{\mathbf{J}_{p_{i}}^{(t^{*})},\forall
t^{*}=0,1,\cdots,t}R_{T}(\mathbf{F}_{b},\mathbf{B}_{b},\mathbf{F}_{u,t}^{(m)},\mathbf{B}_{u,t}^{(m)},\mathbf{F}_{u,r}^{(m)},\cdots,\\\
&\hskip
65.00009pt\cdots,\mathbf{B}_{u,r}^{(m)},\kappa_{b_{i}}^{(t^{*})}\hat{\mathbf{P}}_{b_{i}}^{(t^{*})},\kappa_{u_{i}}^{(t^{*})}\hat{\mathbf{P}}_{u_{i}}^{(t^{*})},\mathbf{X}_{i}^{(t^{*})}),\vspace{-0.5ex}\end{split}$
(29)
$\begin{split}&\mathbf{J}_{p_{\mathrm{best}}}^{(t)}=\operatorname*{arg\,max}_{\mathbf{J}_{p_{\mathrm{best}_{i}}}^{(t)},\forall
i=0,1,\cdots,M_{p}}R_{T}(\mathbf{F}_{b},\mathbf{B}_{b},\mathbf{F}_{u,t}^{(m)},\mathbf{B}_{u,t}^{(m)},\mathbf{F}_{u,r}^{(m)},\cdots,\\\
\vspace{-2ex}&\vspace{-2ex}\hskip
40.00006pt\cdots,\mathbf{B}_{u,r}^{(m)},\kappa_{\mathrm{best},b_{i}}^{(t)}\hat{\mathbf{P}}_{\mathrm{best,b_{i}}}^{(t)},\kappa_{\mathrm{best},u_{i}}^{(t)}\hat{\mathbf{P}}_{\mathrm{best,u_{i}}}^{(t)},\mathbf{X}_{i}^{(t)}).\end{split}\vspace{-5ex}$
(30)

After $T$ iterations, we update
$\mathbf{x}^{(m)}=\mathbf{X}_{\mathrm{best}}^{(T)}$,
$\mathbf{P}_{b}=\kappa_{\mathrm{best},b}^{(T)}\hat{\mathbf{P}}_{\mathrm{best},b}^{(T)}$
and
$\mathbf{P}_{u}^{(m)}=\kappa_{\mathrm{best},u}^{(T)}\hat{\mathbf{P}}_{\mathrm{best},u}^{(T)}$.

(a)

(b)

Figure 3: Achievable rate $\mathrm{R}_{2}$ vs. $(x-y)$-coordinates at
$P_{T}^{(m)}$ = 20 dBm. (a) Single UAV deployment ($M=1$). (b) Multiple UAV
deployment ($M=2$).

##  IV Illustrative Results

In this section, the Monte-Carlo simulation results are presented based on the
proposed scheme. Table I outlines the simulation setup based on the 3D micro-
cell scenario [16] for the results discussed hereafter. The PSO parameters are
chosen as: $N_{p}$$=$$20$, $\gamma_{1}$$=$$\gamma_{2}$$=2$ and
$\gamma_{3}$$=1.1$. In Fig. 3(a), we compare the achievable rate
$\mathrm{R}_{2}$ versus transmit power for the following two cases: 1) a
single UAV ($M=1$) deployed at initial fixed location $(x_{u},y_{u})=(50,50)$;
and 2) multiple UAVs ($M=2$) deployed at initial fixed locations
$(x_{u}^{(1)},y_{u}^{(1)})=(50,50)$, $(x_{u}^{(2)},y_{u}^{(2)})=(100,50)$. We
consider a practical user distribution scenario where the users are located at
multiple locations (i.e., ($x_{k},y_{k}$) $\in$ [50,100]) from BS, which is
located at $(x_{b},y_{b})=(0,0)$. It can be seen from Fig. 3(a) that the
proposed scheme can optimize the UAV location, however, it can only achieve a
sub-optimal solution as a single UAV can not be positioned optimally to
support a large number of users. Moreover, each $k^{th}$ user experience
interference from $K-1$ users, which leads to low achievable rate. To improve
the performance, a multiple UAV-assisted system (i.e., $M=2$) is used in Fig.
3(b), where each UAV can support $K_{m}=K/M$ users based on the proposed
scheme and find the optimal deployments close to its associated $K_{m}$ users.
It can be seen that using the proposed PSO-based UAV location and PA (J-HBF-
PSOLPA) scheme, each UAV can optimally cluster its users, and then find the
optimal deployment while achieving the maximum achievable rate. Thus, a
multiple UAV-assisted MU-mMIMO system can overcome the coverage and capacity
limitation of a single UAV in more practical scenarios.

Figure 4: Total AR $\mathrm{R}_{T}$ vs. $P_{T}$ for single and multiple UAV
system.

Fig. 4 compares the achievable rate of a multiple UAV system versus a single
UAV system for the following cases: 1) J-HBF-PSOLPA with direct link from BS
(we consider optimal PA at each UAV while considering equal PA from BS); 2)
J-HBF-PSOLPA without direct link from BS; and 3) fixed UAV location and equal
PA (FL-EQPA). It can be seen that by using $M=2$ UAV can provide approximately
2.5 times the achievable rate when compared to a single UAV case ($M=1$). For
instance, the achievable rate at $P_{T}=50$ dBm is increased from 27 bps/Hz to
58 bps/Hz, which indicates around 250$\%$ increase in total achievable rate.
Moreover, compared to FL-EQPA, the proposed scheme can significantly enhance
the performance by optimizing the UAV locations and the allocated powers. The
inclusion of direct link from BS to users (i.e., cooperative relaying) can
provide an additional 4 to 5 bps/Hz rate improvement for all cases. Fig. 5
compares the achievable rate for different number of UAVs (i.e., for $M=2$ or
$M=3$) for the following four cases: 1) J-HBF-PSOLPA with direct link and
optimal PA at BS and each UAV; and 2) J-HBF-PSOLPA with direct link and
optimal PA at each UAV only; 3) J-HBF-PSOLPA without direct link; and 4)
J-HBF-FL-EQPA. The analysis can be summarized as follows: 1) the optimal PA
from BS and each UAV can provide an improved performance (e.g., an increase of
$\approx$ 10% rate) when compared to only optimal PA at each UAV; 2) by
increasing the number of UAVs, we can achieve a higher achievable rate (e.g.,
the rate can be increased by 5-7$\%$ when $M$ is increased from 2 to 3).
Moreover, by increasing the number of UAVs, we can further improve the
performance as each UAV cluster less number of users, which leads to reduced
MU interference, and thus, an increased achievable rate. However, it must be
noted that optimizing the number of UAVs is beyond the scope of this paper,
and is left as our future work.

##  V Conclusions

In this paper, we considered a MU-mMIMO IoT cooperative relaying system, where
multiple UAV DF relays connect the BS to a large number of users. We have
proposed a sequential optimization scheme that employed swarm intelligence to
assign users to UAVs through K-means clustering, and optimized the UAV
locations and power allocation from BS and each UAV, followed by the design of
RF and BB stages for maximum achievable rates. The RF stages are designed
using the angular information of UAVs and users, while BB stages are designed
using reduced-dimension effective channel matrices. Our findings show that
multiple UAV-assisted cooperative relaying system works better than single
UAV, especially when taking into account the practical user distributions.
Moreover, compared to fixed positions and equal power distribution of UAV, the
optimization of UAV locations and power allocation substantially improves the
achievable rate.

Figure 5: Total AR $\mathrm{R}_{T}$ vs. $P_{T}$ for 2 and 3 UAV system.

## References

  * [1] U. Uyoata _et al._ , “Relaying in the Internet of Things (IoT): A survey,” _IEEE Access_ , vol. 9, pp. 132 675–132 704, 2021.
  * [2] Q. Wu and R. Zhang, “Common throughput maximization in UAV-enabled OFDMA systems with delay consideration,” _IEEE Trans. Commun._ , vol. 66, no. 12, pp. 6614–6627, 2018.
  * [3] W. Feng _et al._ , “UAV-aided MIMO communications for 5G Internet of Things,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 1731–1740, 2018\.
  * [4] Y. Huang _et al._ , “Massive MIMO for cellular-connected UAV: Challenges and promising solutions,” _IEEE Commun. Mag._ , vol. 59, no. 2, pp. 84–90, 2021.
  * [5] M. Mahmood _et al._ , “Energy-efficient MU-massive-MIMO hybrid precoder design: Low-resolution phase shifters and digital-to-analog converters for 2D antenna array structures,” _IEEE Open J. Commun. Soc._ , vol. 2, pp. 1842–1861, 2021.
  * [6] M. Mahmood _et al._ , “3-D antenna array structures for millimeter wave multi-user massive MIMO hybrid precoder design: A performance comparison,” _IEEE Commun. Lett._ , vol. 26, no. 6, pp. 1393–1397, 2022.
  * [7] M. Mahmood _et al._ , “2D antenna array structures for hybrid massive MIMO precoding,” in _Proc. IEEE Global Commun. Conf. (GLOBECOM)_ , Taipei, Taiwan, Dec. 7-11, 2020.
  * [8] M. Mahmood _et al._ , “PSO-based joint UAV positioning and hybrid precoding in UAV-assisted massive MIMO systems,” in _Proc. IEEE 96th Veh. Technol. Conf. (VTC-Fall)_ , 2022, pp. 1–6.
  * [9] J. Du _et al._ , “Energy-saving UAV-assisted multiuser communications with massive MIMO hybrid beamforming,” _IEEE Commun. Lett._ , vol. 24, no. 5, pp. 1100–1104, 2020.
  * [10] M. Mahmood _et al._ , “Spherical array-based joint beamforming and UAV positioning in massive MIMO systems,” in _Proc. IEEE 97th Veh. Technol. Conf. (VTC-Spring)_ , 2023, pp. 1–5.
  * [11] X. Jiang _et al._ , “Power and trajectory optimization for UAV-enabled amplify-and-forward relay networks,” _IEEE Access_ , vol. 6, pp. 48 688–48 696, 2018.
  * [12] Y. Chen _et al._ , “Multiple UAVs as relays: Multi-hop single link versus multiple dual-hop links,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 9, pp. 6348–6359, 2018.
  * [13] C. Zhang _et al._ , “3D deployment of multiple UAV-mounted base stations for UAV communications,” _IEEE Trans. Commun._ , vol. 69, no. 4, pp. 2473–2488, 2021.
  * [14] N. Nouri _et al._ , “Multi-UAV placement and user association in uplink MIMO ultra-dense wireless networks,” _IEEE Trans. Mobile Comput._ , vol. 22, no. 3, pp. 1615–1632, 2023.
  * [15] K. Xiong _et al._ , “Optimal cooperative beamforming design for MIMO decode-and-forward relay channels,” _IEEE Trans. Signal Process._ , vol. 62, no. 6, pp. 1476–1489, 2014.
  * [16] M. Mahmood _et al._ , “Deep learning meets swarm intelligence for UAV-assisted IoT coverage in massive MIMO,” _IEEE Internet Things J._ , vol. 11, no. 5, pp. 7679–7696, 2024.
  * [17] A. Koc _et al._ , “3D angular-based hybrid precoding and user grouping for uniform rectangular arrays in massive MU-MIMO systems,” _IEEE Access_ , vol. 8, pp. 84 689–84 712, 2020.