# Elevation Angle-Dependent 3D Trajectory Design for Aerial RIS-aided
Communication

Yifan Liu, , Bin Duo, , Qingqing Wu, , Xiaojun Yuan, , Jun Li, , and Yonghui
Li,  Yifan Liu is with the College of Mechanical and Electrical Engineering,
Chengdu University of Technology, Chengdu 610059, China. Bin Duo is with the
College of Computer Science and Cyber Security, Chengdu University of
Technology, Chengdu 610059, China. Qingqing Wu is with the State Key
Laboratory of Internet of Things for Smart City, University of Macau, Macao
999078, China. Xiaojun Yuan is with the National Key Laboratory of Science and
Technology on Communications, University of Electronic Science and Technology
of China, Chengdu 611731, China. Jun Li is with the School of Electronic and
Optical Engineering, Nanjing University of Science and Technology, Nanjing
210094, China. Yonghui Li is with the School of Electrical and Information
Engineering, The University of Sydney, Sydney, NSW 2006, Australia.

###### Abstract

This paper investigates an aerial reconfigurable intelligent surface
(RIS)-aided communication system under the probabilistic line-of-sight (LoS)
channel, where an unmanned aerial vehicle (UAV) equipped with an RIS is
deployed to assist two ground nodes in their information exchange. An
optimization problem with the objective of maximizing the minimum average
achievable rate is formulated to jointly design the communication scheduling,
the RIS’s phase shift, and the three-dimensional (3D) UAV trajectory. To solve
such a non-convex problem, we propose an efficient iterative algorithm to
obtain its suboptimal solution. Simulation results show that our proposed
design significantly outperforms the existing schemes and provides new
insights into the elevation angle and distance trade-off for the UAV-borne RIS
communication system.

###### Index Terms:

UAV communication, reconfigurable intelligent surface, probabilistic LoS
channel, trajectory design.

##  I Introduction

Besides terrestrial deployment, wireless networks are gradually evolving into
air-ground integrated networks to achieve ubiquitous wireless connection and
network capacity upgrades [1, 2]. Recently, unmanned aerial vehicles (UAVs)
have obtained substantial attention in wireless communication. Thanks to their
high mobility, low cost, and line-of-sight (LoS) transmission, UAVs can
further improve communication coverage, throughput, and average secrecy rates
[3, 4, 5].

Recently, the reconfigurable intelligent surface (RIS) has attracted
considerable attention due to its low profile, low energy consumption, and
ability to overcome the non-LoS (NLoS) transmission [6]. Typically, the RIS
contains many reflective elements, each of which is manipulated to induce
changes in the amplitude and phase shift of incident signals to create
favorable propagation environment. The RIS becomes a promising technology for
the future mobile communications. It can solve the pain points of fifth-
generation (5G), such as high energy consumption and coverage voids. By
optimizing the phase shift of all elements of the RIS, the signals from
different transmission paths can be precisely aligned at the desired receiver
to increase the signal power [7].

For typical application scenarios of UAV communications in urban areas, such
as cargo delivery, traffic monitoring, and so on, their communication links
are often blocked by tall building, which leads to severe degradation of
channel quality. Fortunately, with its low power consumption and lightweight,
the RIS can be installed at an appropriate location to reconfigure the
propagation environment of air-ground links, thereby improving communication
performance. Several works have studied various RIS-assisted UAV communication
systems. In general, these studies mainly fall into two categories, one for
terrestrial RIS[8, 9, 10, 11] and the other for aerial RIS[12, 13, 14, 15,
16]. In particular, for the first category, the UAV trajectory and the phase
shift of the RIS mounted on building facades are jointly designed to optimize
different objectives such as communication coverage [8], energy efficiency
[9], confidentiality [10], and communication rates[11].

Taking advantage of the UAV’s ability to fly freely in the three-dimensional
(3D) space, the RIS can be mounted on the UAV. This allows the RIS to fly
along with the UAV, which is more flexible in adjusting its 3D location than
the conventional terrestrial RIS, thus enhancing communication services [12,
13]. An aerial RIS was deployed in [14] to expand the coverage of
communication services, where the worst-case signal-to-noise ratio (SNR) was
improved by jointly optimizing the transmit beamforming, the RIS’s placement,
and the 3D passive beamforming. In [15], multiple users were served by a base
station (BS), and their received powers were significantly enhanced with the
aid of an aerial RIS. To provide communication services for blocked users that
locate far apart while preventing information leakage, the UAV equipped with
one RIS was deployed in [16] to improve the security and energy efficiency.
Above works only took the ideal deterministic LoS channel (DLC) into account,
which has two limitations in practice. $i)$ The DLC model cannot fully capture
the critical effects of the UAV location-dependent LoS and NLoS states in
urban areas with typically high and dense buildings/trees[17]; $ii)$ the DLC
model cannot accurately describe the elevation angle and distance trade-off,
since the elevation angles between the UAV and ground nodes (GNs) are closely
related to the 3D UAV trajectory [18]. Therefore, it is intuitive that the UAV
trajectory designed under the simplified DLC model will undoubtedly cause
significant performance loss in practical urban environments.

Motivated by the above, this paper considers an aerial RIS-aided communication
system, where the UAV-borne RIS assists in information exchange with two GNs.
In particular, we adopt a more accurate probabilistic LoS channel (PLC) model
to characterize the complex channel states of LoS and NLoS in an urban
environment. To maximize the minimum average achievable rate, we jointly
optimize the communication scheduling, the RIS’s phase shift, and the 3D UAV
trajectory. The formulated problem is non-convex and difficult to solve, since
it contains intractable non-convex constraints, binary scheduling variables,
and the complicated achievable rate expression concerning UAV trajectory
variables. To tackle such challenges, we propose an efficient iterative
algorithm to obtain a high-quality solution. Simulation results show that our
proposed joint design for the aerial RIS-aided communication system under the
PLC model can significantly improve the max-min rate compared to that under
the conventional DLC model. This is because the optimized UAV trajectory
better balances the elevation angle and distance trade-off, resulting in the
enhanced gain of the cascaded channel between the UAV and the GNs.

Note that although the authors in [19] and [20] also considered a UAV-borne
RIS communication system under the PLC model, in this paper, both the system
model and the proposed algorithm are different, which results in the
fundamental trade-off between the elevation angle and distance for rate
enhancement. Furthermore, in [17], the 3D UAV trajectory was designed for data
collection under the PLC model; however, this paper considers a new UAV-borne
RIS communication system and derives a new expression for the expected
achievable rate function.

##  II System Model and Problem Formulation

In this paper, we consider an aerial RIS-aided communication system, where two
GNs exchange information via a UAV-borne RIS111Note that the movement of the
UAV changes the orientation of the RIS, which will cause the difference
between the phase shift designed in this paper and the required phase shift in
practice, resulting in performance loss. However, it is reasonable to assume
that the RIS can always remain stable in this paper, since we can employ a
three-axis gimbal to ensure that the orientation of the RIS keeps unchangeable
for the duration of the UAV flight [21]. due to the blockage of dense
buildings. We characterize the position of the UAV and the two GNs via the 3D
Cartesian coordinate system. It is assumed that the UAV flies over a given
duration $T$ to assist in reflecting signals via the RIS between the GNs,
whose locations are denoted by $\mathbf{w}_{k}=[x_{k},y_{k}]$,
$k\in\mathcal{K}=\left\\{1,2\right\\}$. To ease the 3D UAV trajectory design,
$T$ is divided into $N$ time slots that are equal in length, i.e.,
$T=N\delta_{t}$, where $\delta_{t}$ is the length of each time slot. Thus, the
trajectory of the UAV can be approximated by a 3D sequence
$\left\\{\left(\mathbf{q}[n],h[n]\right)\right\\},n\in\mathcal{N}=\left\\{1,2,...,N\right\\}$,
where the discrete way-points $\mathbf{q}[n]=[x[n],y[n]]$ and $h[n]$ represent
the horizontal and vertical locations, which satisfy the following
constraints:

$\|\mathbf{q}[n+1]-\mathbf{q}[n]\|^{2}\leq\hat{\Omega}^{2},\forall n,$ (1)
$\mathbf{q}[N+1]=\mathbf{q}_{F},\mathbf{q}[1]=\mathbf{q}_{0},$ (2)
$|{h}[n+1]-{h}[n]|^{2}\leq\tilde{\Omega}^{2},H_{\rm{min}}\leq h[n]\leq
H_{\rm{max}},\forall n,$ (3)

where $\mathbf{q}_{0}$ and $\mathbf{q}_{F}$ denote the initial and final
horizontal positions of the UAV, respectively,
$\hat{\Omega}=\hat{v}_{\rm{max}}\delta_{t}$ and
$\tilde{\Omega}=\tilde{v}_{\rm{max}}\delta_{t}$ are the maximum horizontal and
vertical distance that the UAV can reach in each time slot, respectively, and
$\hat{v}_{\rm{max}}$ and $\tilde{v}_{\rm{max}}$ are the corresponding maximum
horizontal and vertical flying speed of the UAV, respectively. Furthermore,
$H_{\rm{min}}$ and $H_{{\rm{max}}}$ indicate the minimum and maximum altitude
that the UAV can reach at any given time.

We assume that both GN 1 and GN 2 are equipped with an omni-directional single
antenna. The RIS is equipped with a uniform planar array (UPA) of
$M=M_{x}\times M_{y}$ ($M_{x}$ rows and $M_{y}$ columns) reflective elements,
each of which can be manipulated by an embedded development board, such as
Raspberry Pi $4$, mounted on the UAV. Since the UAV and the RIS are assembled
compactly, we assume that their 3D coordinates are identical, causing
negligible performance loss due to high flying altitude of the UAV [16]. We
denote
$\Theta[n]=\operatorname{diag}\left\\{e^{j\theta_{1,1}[n]},e^{j\theta_{1,2}[n]},\cdots,e^{j\theta_{M_{x},M_{y}}[n]}\right\\}$
as the diagonal phase-shift matrix of the RIS in the $n$th time slot, where
$\theta_{m_{x},m_{y}}[n]\in[0,2\pi)$ is the phase shift of $m_{x}$-th row and
$m_{y}$-th column in the $n$th time slot, where
$m_{x}\in\mathcal{M}_{x}=\\{1,\cdots,M_{x}\\}$ and
$m_{y}\in\mathcal{M}_{y}=\\{1,\cdots,M_{y}\\}.$ Assume that the phase shift of
each element is continuously controllable, which should satisfy
$|e^{j\theta_{m_{x},m_{y}}[n]}|=1$.

To accurately represent the channel states in urban environments, we adopt the
PLC model[22], for which the ground-UAV channel can be represented by LoS or
NLoS state. Thus, for GN $k$, the LoS probability in the $n$th time slot is
given by

$P_{k}^{\mathrm{L}}[n]=\frac{1}{1+ae^{\left(-b\left[\psi_{k}[n]-a\right]\right)}},$
(4)

where _a_ > 0 and _b_ > 0 are constants specified by the actual environment,
and

$\psi_{k}[n]=\frac{180}{\pi}\arctan\left(\frac{h[n]}{\left\|\mathbf{q}[n]-\mathbf{w}_{k}\right\|}\right)$
(5)

is the elevation angle from GN $k$ to the UAV in the $n$th time slot. The
relevant NLoS probability can then be acquired as
$P_{k}^{\mathrm{N}}[n]=1-P_{k}^{\mathrm{L}}[n]$. The channel gain between GN
$k$ and the UAV conditioned on the LoS state in the $n$th time slot can be
expressed as

$\displaystyle\mathbf{h}_{k}^{L}[n]=\tau[n]$
$\displaystyle\left[1,\ldots,e^{-j\frac{2\pi
d}{\lambda}(M_{x}-1){\rm{sin}}\varphi_{k}[n]{\rm
cos}\omega_{k}[n]}\right]^{T}\otimes$
$\displaystyle\left[1,\ldots,e^{-j\frac{2\pi
d}{\lambda}(M_{y}-1){\rm{sin}}\varphi_{k}[n]{\rm
sin}\omega_{k}[n]}\right]^{T},$ (6)

where $\tau[n]=\sqrt{\beta_{0}d_{k}^{-\alpha_{\mathrm{L}}}[n]}$, $\beta_{0}$
is the path loss at the reference distance of $D_{0}=1$ meter (m),
$d_{k}[n]=\sqrt{(\mathbf{q}[n]-\mathbf{w}_{k})^{2}+h[n]^{2}}$ is the distance
from GN $k$ to the UAV in the $n$th time slot, $\alpha_{\mathrm{L}}$ denotes
the path loss exponent for the LoS state. Furthermore, $d$ is the antenna
separation, $\lambda$ is the carrier wavelength, $\varphi_{k}[n]$ and
$\omega_{k}[n]$ represent the elevation and azimuth angles in the $n$th time
slot, respectively. Furthermore, ${\rm{sin}}\varphi_{k}[n]{\rm
cos}\omega_{k}[n]=\frac{x[n]-x_{k}}{d_{k}[n]}$ and
${\rm{sin}}\varphi_{k}[n]{\rm sin}\omega_{k}[n]=\frac{y[n]-x_{k}}{y_{k}[n]}$.
The channel gain between GN $k$ and the UAV conditioned on the NLoS state in
the $n$th time slot is given by

$\mathbf{h}_{k}^{N}[n]=\zeta[n]\mathbf{\widetilde{h}{}}_{k},$ (7)

where $\zeta[n]=\sqrt{\beta_{0}d_{k}^{-\alpha_{\mathrm{N}}}[n]}$,
$\alpha_{\mathrm{N}}$ denotes the path loss exponent for the NLoS state, and
$\mathbf{\widetilde{h}}_{k}\sim\mathcal{CN}(0,1)$ is the small-scale fading
component modeled by a circularly symmetric complex Gaussian (CSCG) random
variable.

Assume that GN $k$ operates in the half-duplex mode, i.e., it can only receive
or transmit in each time slot. Thus we define a binary variable that indicates
whether GN $k$ is scheduled to receive reflected signals from the UAV in the
$n$th time slot or not, i.e., GN $k$ receives signals from the other GN via
the RIS if $\alpha_{k}[n]$ = 1, and transmits otherwise. Assume that only one
GN is allowed to transmit or receive signals to or from the UAV-borne RIS in
the $n$th time slot, so we have the following scheduling constraints:

$\displaystyle\sum_{k=1}^{2}\alpha_{k}[n]\leq 1,\forall n\in\mathcal{N},$ (8)
$\displaystyle\alpha_{k}[n]\in\\{0,1\\},\forall k,n.$ (9)

In a statistical sense, the expected achievable rate in bits/second/Hertz
(bps/Hz) from GN $\bar{k}$ through the UAV to GN ${k}$ in the $n$th time slot
can be given by

$\displaystyle\mathbb{E}\left[R_{k}[n]\right]=$ $\displaystyle
P_{k}^{\mathrm{L}}[n]P_{\bar{k}}^{\mathrm{L}}[n]R_{k}^{\mathrm{LL}}[n]+P_{k}^{\mathrm{L}}[n]P_{\bar{k}}^{\mathrm{N}}[n]R_{k}^{\mathrm{LN}}[n]$
$\displaystyle+$ $\displaystyle
P_{k}^{\mathrm{N}}[n]P_{\bar{k}}^{\mathrm{L}}[n]R_{k}^{\mathrm{NL}}[n]+P_{k}^{\mathrm{N}}[n]P_{\bar{k}}^{\mathrm{N}}[n]R_{k}^{\mathrm{NN}}[n],$
(10)

where
$R_{k}^{{fo}}[n]=\log_{2}\left(1+{\gamma_{\bar{k}}\left|\left(\mathbf{h}_{k}^{f}[n]\right)^{H}\Theta[n]\mathbf{h}_{\bar{k}}^{o}[n]\right|^{2}}\right)$,
$f,o\in\\{L,N\\}$ denote respectively the achievable rates at GN $k$
conditioned on the LoS and NLoS states of the air-ground channels.
Furthermore, $\gamma_{\bar{k}}=P_{\bar{k}}/\sigma^{2}$, where $P_{\bar{k}}$ is
maximum transmit power of GN ${\bar{k}}$ and $\sigma^{2}$ is the variance of
Gaussian noise at GN $k$.

We aim to maximize the minimum average achievable rate by jointly optimizing
the communication scheduling $\mathbf{A}$, the horizontal UAV trajectory
$\mathbf{Q}$, the vertical UAV trajectory $\mathbf{H}$, and the RIS’s phase
shift $\mathbf{\Theta}$ for the entire $N$ time slot. Thus, the optimization
problem can be formulated as

$\displaystyle\max_{\mathbf{A},\mathbf{Q},\mathbf{\Theta},\mathbf{H},\mathbf{\Psi_{k}}}$
$\displaystyle\quad\eta$ (11a) s.t.
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\alpha_{k}[n]\mathbb{E}\left[R_{k}[n]\right]\geq\eta,\forall
k\in\mathcal{K},$ (11b)
$\displaystyle|e^{j\theta_{m_{x},m_{y}}[n]}|=1,\quad\forall n,$ (11c)
$\displaystyle(1)-(3),(5),(8)-(9).$

Problem (11) is non-convex because the rate function $\bar{R}_{{k}}[n]$ in
(11b) is not jointly concave with respect to $\mathbf{A}$, $\mathbf{\Theta}$,
$\mathbf{Q}$, and $\mathbf{H}$, the binary scheduling constraints in (9) are
non-convex, the phase shift constraints in (11c) are non-convex, and the
elevation angle constraints in (5) are non-affine, which make it difficult to
obtain the optimal solution. In the following section, we propose an efficient
iteration algorithm to solve problem (11).

##  III Proposed Algorithm

In this section, an alternating optimization method222The optimization
algorithm can be executed on a ground control center. Specifically, the ground
control center first executes the proposed algorithm and then transmits the
optimized variables to the UAV through the control and non-payload
communication (CNPC) link so that the UAV can perform the mission according to
the designed trajectory [23]. is proposed to obtain a high-quality solution to
problem (11), where $\mathbf{A}$, $\mathbf{\Theta}$, $\mathbf{Q}$, and
$\mathbf{H}$ are iteratively optimized. Specifically, the original problem is
partitioned into three subproblems, each of which is solved in an iterative
manner until the algorithm converges.

###  III-1 GNs’ Scheduling Optimization

With any given feasible $\mathbf{Q}$, $\mathbf{\Theta}$, and $\mathbf{H}$,
this subproblem can be expressed as

$\displaystyle\max_{\mathbf{A},\eta}$ $\displaystyle\quad\eta$ (12a) s.t.
$\displaystyle 0\leq\alpha_{k}[n]\leq 1,\forall k,n,$ (12b)
$\displaystyle(8),(11b).$

Since problem (12) is a standard linear program, it can be solved efficiently
by CVX [24].

###  III-2 RIS’s Phase Shift Design

With any given feasible $\mathbf{A}$, $\mathbf{Q}$, and $\mathbf{H}$, problem
(11) can be rewritten as

$\displaystyle\max_{\Theta,\eta}\eta$ (13) $\displaystyle\text{ s.t.
}(11b)-(11c).$

Since the constraints (11c) are unit modulus, it is difficult to solve problem
(13). To overcome this difficulty, we employ the semidefinite relaxation
(SDR). To begin, let us perform the following transformation:

$\displaystyle\left(\mathbf{h}_{k}^{f}[n]\right)^{H}\Theta[n]\mathbf{h}_{\bar{k}}^{o}[n]=$
$\displaystyle\left(\mathbf{h}_{k}^{f}[n]\right)^{H}\operatorname{diag}(\mathbf{h}_{\bar{k}}^{o}[n])\boldsymbol{v}[n],$
$\displaystyle f,o\in\\{L,N\\},$ (14)

where
$\boldsymbol{v}[n]=[e^{j\theta_{1,1}[n]},e^{j\theta_{1,2}[n]},\cdots,e^{j\theta_{M_{x},M_{y}}[n]}]^{T}$.
According to (14), we can reinterpret the square term in the expected
achievable rate as

$\displaystyle\left|\left(\mathbf{h}_{k}^{f}[n]\right)^{H}\Theta[n]\mathbf{h}_{\bar{k}}^{o}[n]\right|^{2}=\left(\mathbf{h}_{k}^{f}[n]\right)^{H}\operatorname{diag}(\mathbf{h}_{\bar{k}}^{o}[n])\boldsymbol{v}[n]\boldsymbol{v}^{H}[n]$
$\displaystyle\operatorname{diag}(\mathbf{h}_{\bar{k}}^{o}[n])^{H}\left(\mathbf{h}_{k}^{f}[n]\right)=\operatorname{Tr}(\boldsymbol{V}[n]\mathbf{G}^{fo}[n]),\quad
f,o\in\\{L,N\\},$ (15)

where
$\mathbf{G}^{fo}[n]=\operatorname{diag}(\mathbf{h}_{\bar{k}}^{o}[n])^{H}\left(\mathbf{h}_{k}^{f}[n]\right)\left(\mathbf{h}_{k}^{f}[n]\right)^{H}\operatorname{diag}(\mathbf{h}_{\bar{k}}^{o}[n])$.
Then, the expected achievable rate can be re-expressed as

$\displaystyle\mathbb{E}\left[\widetilde{R}_{k}[n]\right]$
$\displaystyle=P_{k}^{\mathrm{L}}[n]P_{\bar{k}}^{\mathrm{L}}[n]\log_{2}\left(1+\operatorname{Tr}(\boldsymbol{V}[n]\mathbf{G}^{LL}[n])\right)$
$\displaystyle+P_{k}^{\mathrm{L}}[n]P_{\bar{k}}^{\mathrm{N}}[n]\log_{2}\left(1+\operatorname{Tr}(\boldsymbol{V}[n]\mathbf{G}^{LN}[n])\right)$
$\displaystyle+P_{k}^{\mathrm{N}}[n]P_{\bar{k}}^{\mathrm{L}}[n]\log_{2}\left(1+\operatorname{Tr}(\boldsymbol{V}[n]\mathbf{G}^{NL}[n])\right)$
$\displaystyle+P_{k}^{\mathrm{N}}[n]P_{\bar{k}}^{\mathrm{N}}[n]\log_{2}\left(1+\operatorname{Tr}(\boldsymbol{V}[n]\mathbf{G}^{NN}[n])\right).$
(16)

Since the expected achievable rate is a function of the newly introduced
variable $\boldsymbol{V}[n]$, we can alleviate the difficulties associated
with using $\boldsymbol{v}[n]$ as the optimization variable. In this respect,
$|\boldsymbol{v}_{m,m}[n]|=1$ can be rewritten as $\boldsymbol{V}[n]\succeq 0$
and $\boldsymbol{V}_{m,m}[n]=1$ for the new optimization variable
$\boldsymbol{V}[n]$. Thus, problem (13) can be rewritten as

$\displaystyle\max_{\mathbf{\Theta},\eta}$ $\displaystyle\quad\eta$ (17) s.t.
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\alpha_{k}[n]\mathbb{E}\left[\widetilde{R}_{k}[n]\right]\geq\eta,\forall
k\in\mathcal{K},$ (17a) $\displaystyle\boldsymbol{V}[n]\succeq 0,$ (17b)
$\displaystyle\boldsymbol{V}_{m,m}[n]=1,m=1,2,...,M.$ (17c)

As problem (17) is a semidefinite programming, it can be solved efficiently by
CVX. However, a rank-one solution may not be obtained. Hence, we recover
$\boldsymbol{v}[n]$ from $\boldsymbol{V}[n]$ using the Gaussian randomization
method, which is similar to that in [10] and thus omitted here for brevity.

###  III-3 UAV Horizontal Trajectory Optimization

With any feasible $\mathbf{H}$ and $\mathbf{A}$ and $\mathbf{\Theta}$ obtained
by solving problem (12) and (17), respectively, the UAV horizontal trajectory
optimization problem can be written as

$\displaystyle\max_{\mathbf{Q},\mathbf{\Psi_{k}},\eta}\eta$ (18)
$\displaystyle\quad\text{ s.t. }(1)-(2),(5),(11b).$

To facilitate the solution to problem (18), we firstly rewrite Eq. (6) as
$\boldsymbol{h}^{L}_{k}[n]=\tau[n]\boldsymbol{h}^{L^{{}^{\prime}}}_{k}[n]$.
Note that not only $\tau[n]$ but also
$\boldsymbol{h}^{L^{{}^{\prime}}}_{k}[n]$ is relevant to the UAV trajectory.
It is observed that $\boldsymbol{h}^{L^{{}^{\prime}}}_{k}[n]$ is complex and
non-linear with respect to the UAV trajectory variables, which makes the UAV
trajectory design intractable. To handle such difficulty, we use the UAV
trajectory of the $(l-1)$ th iteration to obtain an approximate
$\boldsymbol{h}^{L^{{}^{\prime}}}_{k}[n]$ in the $l$ th iteration [10].
Therefore, constraint (11b) can be re-expressed as

$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\alpha_{k}[n]\left(P_{k}^{\mathrm{L}}[n]P_{\bar{k}}^{\mathrm{L}}[n]\log_{2}\left(1+\xi_{k}^{LL}[n]d_{k}^{-\alpha_{\mathrm{L}}}[n]d_{\bar{k}}^{-\alpha_{\mathrm{L}}}[n]\right)\right.$
$\displaystyle\left.+P_{k}^{\mathrm{L}}[n]P_{\bar{k}}^{\mathrm{N}}[n]\log_{2}\left(1+\xi_{k}^{LN}[n]d_{k}^{-\alpha_{\mathrm{L}}}[n]d_{\bar{k}}^{-\alpha_{\mathrm{N}}}[n]\right)\right.$
$\displaystyle\left.+P_{k}^{\mathrm{N}}[n]P_{\bar{k}}^{\mathrm{L}}[n]\log_{2}\left(1+\xi_{k}^{NL}[n]d_{k}^{-\alpha_{\mathrm{N}}}[n]d_{\bar{k}}^{-\alpha_{\mathrm{L}}}[n]\right)\right.$
$\displaystyle\left.+P_{k}^{\mathrm{N}}[n]P_{\bar{k}}^{\mathrm{N}}[n]\log_{2}\left(1+\xi_{k}^{NN}[n]d_{k}^{-\alpha_{\mathrm{N}}}[n]d_{\bar{k}}^{-\alpha_{\mathrm{N}}}[n]\right)\right)\geq\eta,$
(19)

where
$\xi_{k}^{LL}[n]=\frac{p_{\bar{k}}|\beta_{0}|^{2}}{\sigma^{2}}\left|\left(\mathbf{h}_{k}^{L^{{}^{\prime}},(l-1)}[n]\right)^{H}\Theta[n]\mathbf{h}_{\bar{k}}^{L^{{}^{\prime}},(l-1)}[n]\right|^{2}$,
$\xi_{k}^{LN}[n]=\frac{p_{\bar{k}}|\beta_{0}|^{2}}{\sigma^{2}}\left|\left(\mathbf{h}_{k}^{L^{{}^{\prime}},(l-1)}[n]\right)^{H}\Theta[n]\mathbf{\widetilde{h}}_{\bar{k}}\right|^{2}$,
$\xi_{k}^{NL}[n]=\frac{p_{\bar{k}}|\beta_{0}|^{2}}{\sigma^{2}}\left|\left(\mathbf{\widetilde{h}}_{k}\right)^{H}\Theta[n]\mathbf{h}_{\bar{k}}^{L^{{}^{\prime}},(l-1)}[n]\right|^{2}$,
$\xi_{k}^{NN}[n]=\frac{p_{\bar{k}}|\beta_{0}|^{2}}{\sigma^{2}}\left|\left(\mathbf{\widetilde{h}}_{k}\right)^{H}\Theta[n]\mathbf{\widetilde{h}}_{\bar{k}}\right|^{2}$.

To deal with the non-convex constraint (19), we introduce the slack variables
$\mathbf{x}={\\{x_{k}[n],\forall k,n\\}}$, $\mathbf{y}={\\{y_{k}[n],\forall
k,n\\}}$, and $\mathbf{z}={\\{z_{k}[n],\forall k,n\\}}$ into the rate function
(11). Thus, $\mathbb{E}\left[R_{k}[n]\right]$ can be rewritten as

$\displaystyle\mathbb{E}\left[R_{k}[n]\right]$
$\displaystyle=\frac{1}{x_{k}[n]x_{\bar{k}}[n]}\log_{2}\left(1+\frac{\xi_{k}^{LL}[n]}{y_{k}^{\alpha_{\mathrm{L}}/2}[n]y_{\bar{k}}^{\alpha_{\mathrm{L}/2}}[n]}\right)$
$\displaystyle+\frac{1}{x_{k}[n]z_{\bar{k}}[n]}\log_{2}\left(1+\frac{\xi_{k}^{LN}[n]}{y_{k}^{\alpha_{\mathrm{L}}/2}[n]y_{\bar{k}}^{\alpha_{\mathrm{N}/2}}[n]}\right)$
$\displaystyle+\frac{1}{z_{k}[n]x_{\bar{k}}[n]}\log_{2}\left(1+\frac{\xi_{k}^{NL}[n]}{y_{k}^{\alpha_{\mathrm{N}}/2}[n]y_{\bar{k}}^{\alpha_{\mathrm{L}/2}}[n]}\right)$
$\displaystyle+\frac{1}{z_{k}[n]z_{\bar{k}}[n]}\log_{2}\left(1+\frac{\xi_{k}^{NN}[n]}{y_{k}^{\alpha_{\mathrm{N}}/2}[n]y_{\bar{k}}^{\alpha_{\mathrm{N}/2}}[n]}\right),$
(20)

where

$\displaystyle x_{k}[n]\geq 1+ae^{\left(-b\left[\psi_{k}[n]-a\right]\right)},$
(21) $\displaystyle
y_{k}[n]\geq\left\|\mathbf{q}[n]-\mathbf{w}_{k}\right\|^{2}+h[n]^{2},$ (22)
$\displaystyle z_{k}[n]\geq
1+\frac{1}{a}e^{\left(b\left[\phi_{k}[n]-a\right]\right)}.$ (23)

Furthermore,

$\displaystyle\psi_{k}[n]\leq\frac{180}{\pi}\arctan\left(\frac{h[n]}{\left\|\mathbf{q}[n]-\mathbf{w}_{k}\right\|}\right),$
(24)
$\displaystyle\phi_{k}[n]\geq\frac{180}{\pi}\arctan\left(\frac{h[n]}{\left\|\mathbf{q}[n]-\mathbf{w}_{k}\right\|}\right),$
(25)

are the relaxed constraints for the sake of handling the non-affine
constraints (5). We can prove by contradiction that constraints (21)-(25) must
hold with equalities to ensure that the objective value of problem (18) does
not decrease. Note that, after the variable replacement, $\bar{R}_{k}[n]$ in
(20) is jointly convex with respect with $x_{k}[n]$, $y_{k}[n]$, and
$z_{k}[n]$.

Although constraint (24) is non-convex, the right-hand-side (RHS) of (24) is
convex with respect to $||\mathbf{q}[n]-\mathbf{w}_{k}||$. Since the first-
order Taylor approximation of a convex function is a global underestimator, it
can be applied at any local points $x_{k}^{(l)}[n]$, $y_{k}^{(l)}[n]$,
$z_{k}^{(l)}[n]$, and $||\mathbf{q}^{(l)}[n]-\mathbf{w}_{k}||$ in the $l$th
iteration for (20) and (24), i.e., (26) which is shown at the top of the next
page, and

$\displaystyle\psi_{k}[n]\leq\frac{180}{\pi}\arctan\left(\frac{h[n]}{\left\|\mathbf{q}[n]-\mathbf{w}_{k}\right\|}\right)$
$\displaystyle=\frac{180}{\pi}\left(F_{k}^{(l)}[n]-G_{k}^{(l)}[n]\left(\left\|\mathbf{q}[n]-\mathbf{w}_{k}\right\|-\left\|\mathbf{q}^{(l)}[n]-\mathbf{w}_{k}\right\|\right)\right),$
(27)

$\displaystyle\frac{\log_{2}\left(B^{(l)}[n]\right)}{x^{(l)}_{k}[n]x^{(l)}_{\bar{k}}[n]}+\frac{\log_{2}\left(C^{(l)}[n]\right)}{x^{(l)}_{k}[n]z^{(l)}_{\bar{k}}[n]}+\frac{\log_{2}\left(D^{(l)}[n]\right)}{z^{(l)}_{k}[n]x^{(l)}_{\bar{k}}[n]}+\frac{\log_{2}\left(E^{(l)}[n]\right)}{z^{(l)}_{k}[n]z^{(l)}_{\bar{k}}[n]}-\left(\frac{\log_{2}\left(B^{(l)[n]}\right)}{x_{k}^{(l)}[n]\left(x_{\bar{k}}^{(l)}[n]\right)^{2}}+\frac{\log_{2}\left(D^{(l)[n]}\right)}{z_{k}^{(l)}[n]\left(x_{\bar{k}}^{(l)}[n]\right)^{2}}\right)(x_{\bar{k}}[n]-x_{\bar{k}}^{(l)}[n])-\left(\frac{\log_{2}\left(B^{(l)}[n]\right)}{\left(x_{k}^{(l)}[n]\right)^{2}x_{\bar{k}}^{(l)}[n]}\right.$
$\displaystyle\left.+\frac{\log_{2}\left(C^{(l)}[n]\right)}{\left(x_{k}^{(l)}[n]\right)^{2}z_{\bar{k}}^{(l)}[n]}\right)(x_{k}[n]-x_{k}^{(l)}[n])-\left(\frac{\log_{2}\left(C^{(l)}[n]\right)}{x_{k}^{(l)}[n]\left(z_{\bar{k}}^{(l)}[n]\right)^{2}}+\frac{\log_{2}\left(E^{(l)}[n]\right)}{z_{k}^{(l)}[n]\left(z_{\bar{k}}^{(l)}[n]\right)^{2}}\right)\left(z_{\bar{k}}[n]-z_{\bar{k}}^{(l)}[n]\right)-\left(\frac{\log_{2}\left(D^{(l)}[n]\right)}{x_{\bar{k}}^{(l)}[n]\left(z_{{k}}^{(l)}[n]\right)^{2}}+\frac{\log_{2}\left(E^{(l)}[n]\right)}{z_{\bar{k}}^{(l)}[n]\left(z_{{k}}^{(l)}[n]\right)^{2}}\right)\left(z_{{k}}[n]\right.$
$\displaystyle\left.-z_{{k}}^{(l)}[n]\right)-\left(\frac{\alpha_{L}\xi^{LL}_{k}[n]\left(y_{k}^{(l)}[n]\right)^{{\frac{-\alpha_{L}-1}{2}}}\log_{2}e}{2x_{k}^{(l)}[n]x_{\bar{k}}^{(l)}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{\alpha_{L}}{2}}B^{(l)}[n]}+\frac{\alpha_{L}\xi^{LN}_{k}[n]\left(y_{k}^{(l)}[n]\right)^{{\frac{-\alpha_{L}-1}{2}}}\log_{2}e}{2x_{k}^{(l)}[n]z_{\bar{k}}^{(l)}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{\alpha_{N}}{2}}C^{(l)}[n]}+\frac{\alpha_{N}\xi^{NL}_{k}[n]\left(y_{k}^{(l)}[n]\right)^{{\frac{-\alpha_{N}-1}{2}}}\log_{2}e}{2z_{k}^{(l)}[n]x_{\bar{k}}^{(l)}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{\alpha_{L}}{2}}D^{(l)}[n]}\right.$
$\displaystyle\left.+\frac{\alpha_{N}\xi^{NN}_{k}[n]\left(y_{k}^{(l)}[n]\right)^{{\frac{-\alpha_{N}-1}{2}}}\log_{2}e}{2z_{k}^{(l)}[n]z_{\bar{k}}^{(l)}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{\alpha_{N}}{2}}E^{(l)}[n]}\right)(y_{k}[n]-y_{k}^{(l)}[n])-\left(\frac{\alpha_{L}\xi^{LL}_{k}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{-\alpha_{L}-1}{2}}\log_{2}e}{2x_{k}^{(l)}[n]x_{\bar{k}}^{(l)}[n]\left(y_{k}^{(l)}[n]\right)^{\frac{\alpha_{L}}{2}}B^{(l)}[n]}+\frac{\alpha_{N}\xi^{LN}_{k}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{-\alpha_{N}-1}{2}}\log_{2}e}{2x_{k}^{(l)}[n]z_{\bar{k}}^{(l)}[n]\left(y_{k}^{(l)}[n]\right)^{\frac{\alpha_{L}}{2}}C^{(l)}[n]}\right.$
$\displaystyle\left.+\frac{\alpha_{L}\xi^{NL}_{k}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{-\alpha_{L}-1}{2}}\log_{2}e}{2z_{k}^{(l)}[n]x_{\bar{k}}^{(l)}[n]\left(y_{k}^{(l)}[n]\right)^{\frac{\alpha_{N}}{2}}D^{(l)}[n]}+\frac{\alpha_{N}\xi^{NN}_{k}[n]\left(y_{\bar{k}}^{(l)}[n]\right)^{\frac{-\alpha_{N}-1}{2}}\log_{2}e}{2z_{k}^{(l)}[n]z_{\bar{k}}^{(l)}[n]\left(y_{k}^{(l)}[n]\right)^{\frac{\alpha_{N}}{2}}E^{(l)}[n]}\right)(y_{\bar{k}}[n]-y^{(1)}_{\bar{k}}[n]),$
(26)

where
$B^{(l)}[n]=1+\frac{\xi^{LL}_{k}[n]}{\left(y_{k}^{(l)}[n]\right)^{\alpha_{L}/2}\left(y_{\bar{k}}^{(l)}[n]\right)^{\alpha_{L}/2}}$,
$C^{(l)}[n]=1+\frac{\xi^{LN}_{k}[n]}{\left(y_{k}^{(l)}[n]\right)^{\alpha_{L}/2}\left(y_{\bar{k}}^{(l)}[n]\right)^{\alpha_{N}/2}}$,
$D^{(l)}[n]=1+\frac{\xi^{NL}_{k}[n]}{\left(y_{k}^{(l)}[n]\right)^{\alpha_{N}/2}\left(y_{\bar{k}}^{(l)}[n]\right)^{\alpha_{L}/2}}$,
and
$E^{(l)}[n]=1+\frac{\xi^{NN}_{k}[n]}{\left(y_{k}^{(l)}[n]\right)^{\alpha_{N}/2}\left(y_{\bar{k}}^{(l)}[n]\right)^{\alpha_{N}/2}}$.
Furthermore,
$F_{k}^{(l)}[n]=\arctan\left(\frac{h[n]}{\left\|\mathbf{q}^{(l)}[n]-\mathrm{w}_{k}\right\|}\right)$
and
$G_{k}^{(l)}[n]=\frac{h[n]}{\left\|\mathbf{q}^{(l)}[n]-w_{k}\right\|^{2}+H^{2}}.$

With (26)-(27), problem (18) can be reformulated into the following convex
optimization problem,

$\displaystyle\max_{\mathbf{Q},\Psi_{k},\Phi_{k}\mathbf{x}_{k},\mathbf{z}_{k},\mathbf{y}_{k},\eta}\eta$
(28) $\displaystyle\text{ s.t. }(1)-(2),(21)-(23),(25)-(27),$

where $\Psi_{k}=\left\\{\psi_{k}[n],\forall k,n\right\\}$ and
$\Phi_{k}=\left\\{\phi_{k}[n],\forall k,n\right\\}$. As such, problem (28) can
be efficiently solved by $\mathrm{CVX}$ [24].

###  III-4 UAV Vertical Trajectory Optimization

With the optimal $\mathbf{A}$, $\mathbf{Q}$, and $\mathbf{\Theta}$ obtained by
problem (12), (17), and (28), the UAV vertical trajectory optimization problem
can be reformulated as

$\displaystyle\max_{\mathbf{H},\mathbf{\Psi_{k}},\eta}\eta$ (29)
$\displaystyle\text{ s.t. }(3),(5),(11b),$

Since problems (18) and (29) are similar in form and differ only slightly in
terms of optimization variables $\mathbf{H}$, the procedure for solving
problem (18) can be similarly applied to solve problem (29). We omit the
detailed derivation owing to the page limitation.

###  III-5 Overall Algorithm

By applying our proposed algorithm, problem (15) can be solved by alternately
optimizing variables $\mathbf{A}$, $\mathbf{\Theta}$, $\mathbf{Q}$, and
$\mathbf{H}$, while its solution converges to a preset accuracy $\epsilon$.
Note that the binary solution can be reconstructed with high precision from
the obtained continuous variables of GN’s transmission scheduling by applying
the proposed reconstruction method in [4]. Since the four subproblems are
solved by applying CVX via the standard interior point method, their
computational complexity can be obtained as
$O_{1}\left((KN)^{3.5}\text{log}(1/\epsilon)\right)$,
$O_{2}\left(\sqrt{M}\text{log}(1/\epsilon)\left((KN+1)M^{3}+(K\right.\right.\\\
\left.\left.N+1)^{2}M^{2}+(KN+1)^{3}\right)\right)$,
$O_{3}\left((2N+6KN)^{3.5}\text{log}(1/\epsilon)\right)$, and
$O_{4}\left((2N+6KN)^{3.5}\text{log}(1/\epsilon)\right)$, respectively. Thus,
the total computational complexity of our proposed algorithm is in the order
of $O_{1}+O_{2}+O_{3}+O_{4}$.

Figure 1: UAV’s trajectories.

Figure 2: Achieved expected max-min
rate versus $T$.

Figure 3: Achieved expected max-min
rate versus $M$.

##  IV Simulation Results

In this section, we provide simulation results to show the interesting
elevation angle and distance trade-off of the 3D UAV trajectory design under
the PLC model (denoted as PLC for brevity). The following schemes are used for
comparison: 1) UAV horizontal trajectory design under the PLC scheme with
fixed altitude being $H=200$ m (denoted as PLCFA); 2) 3D UAV trajectory design
under the LC model (denoted as DLC). We assume that GN1 and GN 2 are located
at (0, 0, 0) m and (800, 0, 0) m, respectively. Unless otherwise stated, the
simulation results are set as: $\hat{v}_{max}=40$ m/s, $\tilde{v}_{max}=20$
m/s $p_{1}=p_{2}=0.1\mathrm{~{}W}$, $\beta_{0}=-40$ $\mathrm{dB}$,
$H_{\rm{min}}=100$ m, $H_{\rm{max}}=500$ m, $d=\frac{\lambda}{2}$,
$\sigma^{2}=-169$ $\mathrm{dBm}$, $\alpha_{\mathrm{L}}=2.2$,
$\alpha_{\mathrm{N}}=3.2$, $a=11.95$, and $b=0.14$ [17, 18]. Furthermore, the
initialized trajectory of the UAV is set to be a straight-line trajectory with
a fix altitude $H=200$ m, i.e., the UAV flies from
$\mathbf{q}_{0}=[-200,-200]$ to $\mathbf{q}_{F}=[1000,200]$ at its maximum
flying speed.

Fig. 2 shows the different 3D UAV trajectories by three schemes with $T=150$
s. It can be seen that the UAV in the DLC first descends quickly to
$H_{\rm{min}}$, hovers over GN 1, then flies horizontally at its maximum
speed, and hovers over GN 2. Finally, the UAV rises and flies back to the
final location. This is because hovering above each GN at $H_{\rm{min}}$
suffers from the least path loss for the cascade channel between the UAV and
the GNs. By contrast, in the PLC scheme, the UAV first ascends rapidly to
increase the elevation angle between the UAV and each GN for a higher LoS
probability. Then, the UAV hovers above the midpoint of the two GNs, which
results in less path loss while maintaining a larger LoS probability of the
cascaded channel. Compared to the PLC scheme that the UAV hovers only above
the midpoints of the two GNs, the UAV in the PLCFA scheme hovers close to both
sides of the midpoint of the two GNs. This is because that although hovering
at the midpoint of the two GNs can maximize the LoS probability of the
cascaded channel, i.e., $P_{k}^{\mathrm{L}}[n]P_{\bar{k}}^{\mathrm{L}}[n]$,
the larger path loss of the communication links also degrades the rate
performance. Therefore, our proposed PLC scheme can take advantage of the
additional design brought by the UAV vertical trajectory to obtain a more
efficient angle-distance trade-off than the PLCFA scheme.

Fig. 3 illustrates the achieved expected max-min rates of different schemes
versus $T$ when $M=100$. It can be seen that the PLCFA scheme achieves larger
rates than the DLC scheme, which indicates the necessity of adopting the more
accurate PLC model to describe the LoS/NLoS channel states in the UAV-borne
RIS communication system. Furthermore, our proposed PLC scheme significantly
improves performance over the PLCFA scheme. The reason is that the
additionally designing the UAV vertical trajectory can further increase the
LoS probability of the cascaded channel. Fig. 4 shows the achieved expected
max-min rates for different schemes with $T=150$ s versus different $M$. As
expected, the rate performance is significantly increased when more elements
are equipped in the RIS due to the larger passive beamforming gain. In
practice, we cannot increase $M$ indefinitely to obtain higher rates due to
the limitation of the size of rotary-wing UAVs. An oversized $M$ will result
in a larger RIS size, which will increase the weight of the UAV and result in
greater energy consumption (i.e., shorter endurance). Therefore, an
interesting trade-off exists between $M$ (i.e., related to the RIS size) and
the UAV energy consumption for rate improvement.

##  V Conclusion

This paper investigated the potential of rate enhancement for the aerial RIS-
aided communication system under the accurate PLC model in the dense urban
environment. The objective was to maximize the minimum average achievable
rate. We proposed an efficient iterative algorithm to jointly optimize the
communication schedule, the RIS’s phase shift, and the 3D UAV trajectory.
Numerical results showed that the proposed scheme has a significant
improvement compared to the conventional DLC scheme. Furthermore, our proposed
scheme enjoys the additional gain of elevation angle-dependent 3D UAV
trajectory design and can effectively balance the elevation angle and distance
trade-off between the UAV and the GNs, whereas the DLC scheme cannot. This
validates the practical importance of considering the more accurate PLC model
to support UAV-borne RIS communications in urban environments.

Note that in the suburban/urban/dense urban environment, the error in the
approximated PLC model [22] is extremely small and thus can be ignored.
However, this PLC model suffers from errors in the high-rise urban
environment, which may cause performance loss. In this case, the hybrid
_offline-online_ 3D UAV trajectory design proposed in [17] may be an
alternative to better characterize the real-time location-dependent air-ground
channel states, which is interesting to study whether the performance loss can
be further compensated by employing additional online design.

## References

  * [1] Y. Zeng, Q. Wu, and R. Zhang, “Accessing from the sky: A tutorial on UAV communications for 5G and beyond,” _Proc. IEEE_ , vol. 107, no. 12, pp. 2327–2375, 2019.
  * [2] R. Zhou, X. Zhang, X. Wang, G. Yang, H.-N. Dai, and M. Liu, “Device-oriented keyword searchable encryption scheme for cloud-assisted industrial IoT,” _IEEE Internet Things J._ , pp. 1–1, 2021.
  * [3] G. Zhang, Q. Wu, M. Cui, and R. Zhang, “Securing UAV communications via joint trajectory and power control,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 2, pp. 1376–1389, 2019.
  * [4] Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for Multi-UAV enabled wireless networks,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 2109–2121, 2018.
  * [5] B. Duo, Q. Wu, X. Yuan, and R. Zhang, “Energy efficiency maximization for full-duplex UAV secrecy communication,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 4, pp. 4590–4595, 2020.
  * [6] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Trans. Commun._ , vol. 58, no. 1, pp. 106–112, 2020.
  * [7] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [8] L. Yang, F. Meng, J. Zhang, M. O. Hasna, and M. D. Renzo, “On the performance of RIS-assisted dual-hop UAV communication systems,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 9, pp. 10 385–10 390, 2020.
  * [9] X. Liu, Y. Liu, and Y. Chen, “Machine learning empowered trajectory and passive beamforming design in UAV-RIS wireless networks,” _IEEE J. Select. Areas Commun._ , vol. 39, no. 7, pp. 2042–2055, 2021.
  * [10] S. Li, B. Duo _et al._ , “Robust secure UAV communications with the aid of reconfigurable intelligent surfaces,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2021.
  * [11] S. Li _et al._ , “Reconfigurable intelligent surface assisted UAV communication: Joint trajectory design and passive beamforming,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 5, pp. 716–720, 2020.
  * [12] A. S. Abdalla _et al._ , “UAVs with reconfigurable intelligent surfaces: Applications, challenges, and opportunities,” 2020. [Online]. Available: https://arxiv.org/abs/2012.04775
  * [13] B. Shang, R. Shafin, and L. Liu, “UAV swarm-enabled aerial reconfigurable intelligent surface,” 2021. [Online]. Available: https://arxiv.org/abs/2103.06361
  * [14] H. Lu _et al._ , “Aerial intelligent reflecting surface: Joint placement and passive beamforming design with 3D beam flattening,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 7, pp. 4128–4143, 2021.
  * [15] S. Jiao _et al._ , “Joint beamforming and phase shift design in downlink UAV networks with IRS-assisted NOMA,” _J. Commun. Inf. Networks_ , vol. 5, no. 2, pp. 138–149, 2020.
  * [16] H. Long _et al._ , “Reflections in the sky: Joint trajectory and passive beamforming design for secure UAV networks with reconfigurable intelligent surface,” 2020. [Online]. Available: https://arxiv.org/abs/2005.10559
  * [17] C. You and R. Zhang, “Hybrid offline-online design for UAV-enabled data harvesting in probabilistic LoS channels,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 6, pp. 3753–3768, 2020.
  * [18] B. Duo _et al._ , “Anti-jamming 3D trajectory design for UAV-enabled wireless sensor networks under probabilistic LoS channel,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 12, pp. 16 288–16 293, 2020.
  * [19] T. Shafique _et al._ , “Optimization of wireless relaying with flexible UAV-borne reflecting surfaces,” _IEEE Trans. Commun._ , vol. 69, no. 1, pp. 309–325, 2021.
  * [20] A. Mahmoud _et al._ , “Intelligent reflecting surfaces assisted UAV communications for IoT networks: Performance analysis,” _IEEE Trans. Green Commun. Networking_ , vol. 5, no. 3, pp. 1029–1040, 2021.
  * [21] F. Zhang, “Simultaneous self-calibration of nonorthogonality and nonlinearity of cost-effective multiaxis inertially stabilized gimbal systems,” _IEEE Rob. Autom. Lett_ , vol. 3, no. 1, pp. 132–139, 2018.
  * [22] A. Al-Hourani _et al._ , “Optimal LAP altitude for maximum coverage,” _IEEE Wireless Commn. Lett._ , vol. 3, no. 6, pp. 569–572, 2014.
  * [23] B. Duo, H. Hu, Y. Li, Y. Hu, and X. Zhu, “Robust 3D trajectory and power design in probabilistic LoS channel for UAV-enabled cooperative jamming,” _Veh. Commun._ , vol. 32, p. 100387, 2021.
  * [24] B. S. Grant, M. CVX: Matlab Software for Disciplined Convex Programming, Version 2.2. (2020). [Online]. Available: http://cvxr.com/cvx


# Evolutionary behavior in a two-locus system

A. M. Diyorov, U. A. Rozikov  A. M. Diyorov
The Samarkand branch of TUIT, Samarkand, Uzbekistan<EMAIL_ADDRESS>U.A.
Rozikova,b,c a V.I.Romanovskiy Institute of Mathematics, 9, Universitet str.,
100174, Tashkent, Uzbekistan; b AKFA University, 264, Milliy Bog street,
Yangiobod QFY, Barkamol MFY, Kibray district, 111221, Tashkent region,
Uzbekistan; c National University of Uzbekistan, 4, Universitet str., 100174,
Tashkent, Uzbekistan<EMAIL_ADDRESS>
###### Abstract.

In this short note we study a dynamical system generated by a two-parametric
quadratic operator mapping 3-dimensional simplex to itself. This is an
evolution operator of the frequencies of gametes in a two-locus system. We
find the set of all (a continuum set) fixed points and show that each fixed
point is non-hyperbolic. We completely describe the set of all limit points of
the dynamical system. Namely, for any initial point (taken from the
3-dimensional simplex) we find an invariant set containing the initial point
and a unique fixed point of the operator, such that the trajectory of the
initial point converges to this fixed point.

Mathematics Subject Classifications (2010). 37N25, 92D10.

Key words. loci, gamete, dynamical system, fixed point, trajectory, limit
point.

##  1\. Introduction

In this paper following [2, page 68] we define an evolution operator of a
population assuming viability selection, random mating and discrete non-
overlapping generations. Consider two loci $A$ (with alleles $A_{1}$, $A_{2}$)
and $B$ (with alleles $B_{1}$, $B_{2}$). Then we have four gametes:
$A_{1}B_{1}$, $A_{1}B_{2}$, $A_{2}B_{1}$, and $A_{2}B_{2}$. Denote the
frequencies of these gametes by $x$, $y$, $u$, and $v$ respectively. Thus the
vector $(x,y,u,v)$ can be considered as a state of the system, and therefore,
one takes it as a probability distribution on the set of gametes, i.e. as an
element of 3-dimensional simplex, $S^{3}$. Recall that $(m-1)$-dimensional
simplex is defined as

$S^{m-1}=\\{x=(x_{1},...,x_{m})\in\mathbb{R}^{m}:x_{i}\geq
0,\sum^{m}_{i=1}x_{i}=1\\}.$

Following [2, Section 2.10] we define the frequencies
$(x^{\prime},y^{\prime},u^{\prime},v^{\prime})$ in the next generation as

$W:\begin{array}[]{llll}x^{\prime}=x+a\cdot(yu-xv)\\\\[5.69054pt]
y^{\prime}=y-a\cdot(yu-xv)\\\\[5.69054pt] u^{\prime}=u-b\cdot(yu-
xv)\\\\[5.69054pt] v^{\prime}=v+b\cdot(yu-xv),\end{array}$ (1.1)

where $a,b\in[0,1]$.

It is easy to see that this quadratic operator, $W$, maps $S^{3}$ to itself.
Indeed, we have $x^{\prime}+y^{\prime}+u^{\prime}+v^{\prime}=1$ and each
coordinate is non-negative, for example, we check it for $y^{\prime}$:

$y^{\prime}=y-a\cdot(yu-xv)=y(1-au)+axv\geq y(1-au)\geq 0,$

these inequalities follow from the conditions that $x,y,u,v,a\in[0,1]$, and
therefore, we have $0\leq au\leq 1$.

The operator (1.1), for any initial point (state)
$t_{0}=(x_{0},y_{0},u_{0},v_{0})\in S^{3}$, defines its trajectory:
$\\{t_{n}=(x_{n},y_{n},u_{n},v_{n})\\}_{n=0}^{\infty}$ as

$t_{n}=(x_{n},y_{n},u_{n},v_{n})=W^{n}(t_{0}),n=0,1,2,...$

Here $W^{n}$ is the $n$-fold composition of $W$ with itself:

$W^{n}(\cdot)=\underbrace{W(W(W\dots(W}_{n\,{\rm times}}(\cdot)))\dots).$

The main problem in theory of dynamical system (see [1]) is to study the
sequence $\\{t_{n}\\}_{n=0}^{\infty}$ for each initial point $t_{0}\in S^{3}$.

In general, if a dynamical system is generated by a nonlinear operator then
complete solution of the main problem may be very difficult. But in this short
note we will completely solve this main problem for nonlinear operator (1.1).

######  Remark 1.

Using $1=x+y+u+v$ (on $S^{3}$) one can rewrite operator (1.1) as

$W:\begin{array}[]{llll}x^{\prime}=x^{2}+xy+xu+(1-a)xv+ayu\\\\[5.69054pt]
y^{\prime}=xy+y^{2}+(1-a)yu+axv+yv\\\\[5.69054pt]
u^{\prime}=xu+(1-b)yu+u^{2}+bxv+uv\\\\[5.69054pt]
v^{\prime}=(1-b)xv+yv+byu+uv+v^{2}.\end{array}$ (1.2)

Note that the operator (1.2) is in the form of quadratic stochastic operator
(QSO), i.e., $V:S^{m-1}\to S^{m-1}$ defined by

$V:x_{k}^{\prime}=\sum_{i,j=1}^{m}P_{ij,k}x_{i}x_{j},$

where $P_{ij,k}\geq 0$, $\sum_{k}P_{ij,k}=1$.

The operator is not studied in general, but some large class of QSO’s are
studied (see for example [3], [4], [5], [6], [7], [8] and the references
therein). But the operator (1.1) was not studied yet.

##  2\. The set of limit points

######  Remark 2.

The case $a=b=0$ is very trivial, so we will not consider this case.

Recall that a point $t\in S^{3}$ is called a fixed point for $W:S^{3}\to
S^{3}$ if $W(t)=t$.

Denote the set of all fixed points by Fix$(W)$.

It is easy to see that for any $a,b\in[0,1]$, $a+b\neq 0$ the set of all fixed
points of (1.1) is

${\rm Fix}(W)=\\{t=(x,y,u,v)\in S^{3}:yu-xv=0\\}.$

This is a continuum set of fixed points.

The main problem is completely solved in the following result:

######  Theorem 1.

For any initial point $(x_{0},y_{0},u_{0},v_{0})\in S^{3}$ the following
assertions hold

  * 1.

If $(x_{0}+y_{0})(u_{0}+v_{0})=0$ then $(x_{0},y_{0},u_{0},v_{0})$ is fixed
point.

  * 2.

If $(x_{0}+y_{0})(u_{0}+v_{0})\neq 0$ then trajectory has the following limit:

$\lim_{n\to\infty}(x_{n},y_{n},u_{n},v_{n})=$
$\left(A(x_{0},u_{0})(x_{0}+y_{0}),A(y_{0},v_{0})(x_{0}+y_{0}),A(x_{0},u_{0})(u_{0}+v_{0}),A(y_{0},v_{0})(u_{0}+v_{0})\right)\in{\rm
Fix}(W),$

where

$A(x,u)={bx+au\over(u_{0}+v_{0})a+(x_{0}+y_{0})b}.$

###### Proof.

We note that for each $\alpha\in[0,1]$ the following set is invariant:

$X_{\alpha}=\\{t=(x,y,u,v)\in S^{3}:x+y=\alpha,\ \ u+v=1-\alpha\\},$

i.e., $W(X_{\alpha})\subset X_{\alpha}$.

Note also that

$S^{3}=\bigcup_{\alpha\in[0,1]}X_{\alpha}.$

The part 1 of theorem follows in the case $\alpha=0$ and $\alpha=1$. Indeed,
for $\alpha=0$ we have

$X_{0}=\\{t=(0,0,u,v)\in S^{3}:u+v=1\\},$

and in the case of $\alpha=1$ we get

$X_{1}=\\{t=(x,y,0,0)\in S^{3}:x+y=1\\}.$

Note that in both case the restriction of operator on the corresponding set is
an id-operator, i.e., all points of the set are fixed points.

Now to prove part 2 we consider the case $\alpha\in(0,1)$.

Since $X_{\alpha}$ is an invariant, it suffices to study limit points of the
operator on sets $X_{\alpha}$, for each $\alpha\in(0,1)$ separately. To do
this, we reduce operator $W$ on the invariant set $X_{\alpha}$ (i.e., replace
$y=\alpha-x$, $v=1-\alpha-u$):

$W_{\alpha}:\begin{array}[]{ll}x^{\prime}=(1-a+a\alpha)x+a\alpha
u\\\\[5.69054pt] u^{\prime}=(1-\alpha)bx+(1-b\alpha)u,\end{array}$ (2.1)

where $a,b\in[0,1]$, $\alpha\in(0,1)$ $x\in[0,\alpha]$, $u\in[0,1-\alpha]$.

It is easy to find the set of all fixed points:

${\rm
Fix}(W_{\alpha})=\\{(x,u)\in[0,\alpha]\times[0,1-\alpha]:(1-\alpha)x-\alpha
u=0\\}.$

The operator $W_{\alpha}$ is a linear operator given by the matrix

$M_{\alpha}=\left(\begin{array}[]{cc}1-a+a\alpha&a\alpha\\\\[5.69054pt]
(1-\alpha)b&1-b\alpha\end{array}\right).$ (2.2)

Eigenvalues of the linear operator are

$\lambda_{1}=1,\ \ \lambda_{2}=1-(1-\alpha)a-\alpha b.$ (2.3)

For any $a,b\in[0,1]$, $a+b\neq 0$, $\alpha\in(0,1)$ we have
$0<(1-\alpha)a+\alpha b<1$, therefore, $0<\lambda_{2}<1.$

By (2.1) we define trajectory of an initial point $(x_{0},u_{0})$ as

$(x_{n+1},u_{n+1})=M_{\alpha}(x_{n},u_{n})^{T},\ \ n\geq 0.$

Thus

$(x_{n},u_{n})=M_{\alpha}^{n}\,(x_{0},u_{0})^{T},\ \ n\geq 1.$ (2.4)

Therefore we need to find $M_{\alpha}^{n}$. To find it we use a little Cayley-
Hamilton Theorem111https://www.freemathhelp.com/forum/threads/formula-for-
matrix-raised-to-power-n.55028/ to obtain the following formula

$M_{\alpha}^{n}={\lambda_{2}\,\lambda^{n}_{1}-\lambda_{1}\,\lambda^{n}_{2}\over\lambda_{2}-\lambda_{1}}\cdot
I_{2}+{\lambda^{n}_{2}-\lambda^{n}_{1}\over\lambda_{2}-\lambda_{1}}\cdot
M_{\alpha},$

where $I_{2}$ is $2\times 2$ unit matrix and $\lambda_{1},\lambda_{2}$ are
eigenvalues (defined in (2.3)).

By explicit formula (2.3) we get the following limit

$\lim_{n\to\infty}M_{\alpha}^{n}={\lambda_{2}\over\lambda_{2}-\lambda_{1}}\cdot
I_{2}-{1\over\lambda_{2}-\lambda_{1}}\cdot
M_{\alpha}={1\over(1-\alpha)a+\alpha b}\cdot\left(\begin{array}[]{cc}\alpha
b&\alpha a\\\\[5.69054pt] (1-\alpha)b&(1-\alpha)a\end{array}\right).$

Using this limit, for any initial point
$(x_{0},u_{0})\in[0,\alpha]\times[0,1-\alpha]$ we get

$\lim_{n\to\infty}(x_{n},u_{n})=\lim_{n\to\infty}M_{\alpha}^{n}\,(x_{0},u_{0})^{T}={bx_{0}+au_{0}\over(1-\alpha)a+\alpha
b}\cdot(\alpha,1-\alpha)\in{\rm Fix}(W_{\alpha}).$ (2.5)

By (2.5) we obtain

######  Lemma 1.

For any initial point $(x_{0},y_{0},u_{0},v_{0})\in S^{3}\setminus(X_{0}\cup
X_{1})$ there exists $\alpha\in(0,1)$ such that $(x_{0},y_{0},u_{0},v_{0})\in
X_{\alpha}$ and the trajectory of this initial point (under operator $W$,
defined in (1.1)) has the following limit

$\lim_{n\to\infty}(x_{n},y_{n},u_{n},v_{n})=$
$\left(A(x_{0},u_{0})\alpha,A(y_{0},v_{0})\alpha,A(x_{0},u_{0})(1-\alpha),A(y_{0},v_{0})(1-\alpha)\right)\in{\rm
Fix}(W),$

where

$A(x,u)={bx+au\over(1-\alpha)a+\alpha b}.$

In this lemma we note that $\alpha=x_{0}+y_{0}$ and $1-\alpha=u_{0}+v_{0}$,
therefore, the part 2 of Theorem follows, where limit point of trajectory of
each initial point is given as function of the initial point only. Theorem is
proved. ∎

##  3\. Biological interpretations

The results of Theorem 1 have the following biological interpretations:

Let $t=(x_{0},y_{0},u_{0},v_{0})\in\mathcal{S}^{3}$ be an initial state (the
probability distribution on the set
$\\{A_{1}B_{1},A_{1}B_{2},A_{2}B_{1},A_{2}B_{2}\\}$ of gametes). Theorem 1
says that, as a rule, the population tends to an equilibrium state with the
passage of time.

Part 1 of Theorem 1 means that if at an initial time we had only two gametes
then the (initial) state remains unchanged.

Part 2 means that depending on the initial state future of the population is
stable: gametes survive with probability

$A(x_{0},u_{0})(x_{0}+y_{0}),A(y_{0},v_{0})(x_{0}+y_{0}),A(x_{0},u_{0})(u_{0}+v_{0}),A(y_{0},v_{0})(u_{0}+v_{0})$

respectively. From the existence of the limit point of any trajectory and from
the explicit form of ${\rm Fix}(W)$ it follows that

$\lim_{n\to\infty}(y_{n}u_{n}-x_{n}v_{n})=0.$

This property, biologically means ([2, page 69]), that the population
asymptotically goes to a state of linkage equilibrium with respect to two
loci.

##  Acknowledgements

Rozikov thanks Institut des Hautes Études Scientifiques (IHES), Bures-sur-
Yvette, France for support of his visit to IHES. The work was partially
supported by a grant from the IMU-CDC.

## References

  * [1] R.L. Devaney, An introduction to chaotic dynamical system, Westview Press, 2003.
  * [2] W.J. Ewens, Mathematical population genetics. Mathematical biology, Springer, 2004.
  * [3] R.N. Ganikhodzhaev, F.M. Mukhamedov, U.A. Rozikov, Quadratic stochastic operators and processes: results and open problems. Inf. Dim. Anal. Quant. Prob. Rel. Fields. 2011. V.14, No.2, p.279-335.
  * [4] Yu.I. Lyubich, Mathematical structures in population genetics. Biomathematics, 22, Springer-Verlag, 1992.
  * [5] U.A. Rozikov, Population dynamics: algebraic and probabilistic approach. World Sci. Publ. Singapore. 2020, 460 pp.
  * [6] U.A. Rozikov, N.B. Shamsiddinov, On non-Volterra quadratic stochastic operators generated by a product measure. Stoch. Anal. Appl. 2009, V.27, No.2, p.353-362.
  * [7] U.A. Rozikov, A.Zada, On $\ell$\- Volterra quadratic stochastic operators. Inter. Journal Biomath. 2010, V. 3, No. 2, p. 143–159.
  * [8] U.A. Rozikov, U.U. Zhamilov, On dynamics of strictly non-Volterra quadratic operators on two-dimensional simplex. Sbornik: Math. 2009, V. 200, No.9, p.1339-1351.


LU TP 97-11

Statistical Properties of Unrestricted
Crew Scheduling Problems

Martin<EMAIL_ADDRESS>Carsten<EMAIL_ADDRESS>and Bo<EMAIL_ADDRESS>
Complex Systems Group, Department of Theoretical Physics
University of Lund, Sölvegatan 14A, S-223 62 Lund, Sweden

Submitted to Operations Research

Abstract:

A statistical analysis is performed for a random unrestricted local crew
scheduling problem, expressed in terms of pairing arrivals with departures.
The analysis is aimed at understanding the structure of similar problems with
global restrictions, and estimating their difficulty. The methods developed
are of a general nature and can be of use in other problems with a similar
structure. For large random problems, the ground-state energy scales like
$\sqrt{N}$ and the average excitation like $N$, where $N$ is the number of
arrivals/departures. The average ground-state degeneracy is such that the
probability of hitting an optimal pairing by chance scales like $2N2^{-N}$ for
large $N$. By insisting on the local ground-state energy for a restricted
problem, airports can be split into smaller parts, and the state space reduced
by typically a factor $\sim 2^{N_{a}}$, with $N_{a}$ the total number of
airports.

##  1 Introduction

Airline crew scheduling represents an important class of optimization problems
where the topological structure is important. In ref. [1], a novel Potts
artificial neural network approach is developed for attacking semi-realistic
airline crew scheduling problems of the following type: A weekly flight
schedule is given in terms of a set of flights, each with a specified airport
and time for departure and arrival. The object is to assign a crew to each
flight, while minimizing a cost-function defined by the total required crew
time (including the waiting-time at airports). The solution is subject to a
set of global constraints: The crews are required to travel along closed
tours, starting and ending at a certain airport, the home-base. These tours
are subject to limitations as to duration and leg-count.

As shown in ref. [1], a great deal of simplification is gained by
reformulating the problem as that of mapping arrivals on departures at each
airport, implying an implicit representation of the crews. In fact, without
the global restrictions, the problem is reduced to a set of independent local
subproblems, one at each airport. Each local problem amounts to minimizing the
local waiting-time, and is simply solvable in polynomial time.

In this paper, we focus on this kind of unrestricted local problems, in
particular their statistical properties. These are quite interesting, and in
no way trivial, in spite of the triviality of the problems. In particular, we
consider the ensemble of random local problems of a fixed size $N$, as defined
by the number of arrivals/departures. In addition, we analyze the properties
of random solutions to such problems.

Such a statistical analysis of this type of problem does not exist in the
literature, and we feel it is interesting for the following reasons: The
results illuminate the structure of the corresponding restricted problems, and
as a by-product, useful tools are provided for probing their difficulty, and
for simplifying their solution. Some of the methods used are novel, and may be
used also in other contexts, where a similar structure occurs. In addition, a
lower bound to the waiting-time is provided by the solutions to the
unrestricted problem. This bound is often saturated [1].

The methodology we use contains the following steps. First, the analysis of a
local problem is simplified by considering its topology (defined by the
relative ordering in time between arrivals and departures) separately from its
geometry (defined by the lengths of the time intervals between consecutive
events). The ensemble of problems is thus factorized into the direct product
of the ensemble of topologies and the ensemble of geometries.

After introducing a notation for the topology, we consider problems with a
fixed topology, and evaluate averages over the geometry of entities related to
the waiting-time. These are simple, since the effect of the geometry on the
waiting-time spectrum is a mere shift.

The apparent difficulty of a problem is probed by analyzing the distribution
of waiting-times of random solutions. A nice feature is that the waiting-time
spectrum for each problem is quantized in steps of the basic period of the
schedule.

Subsequently, all variables of interest are averaged also over the topology.
This is a more difficult task, and requires the use of a subtle recursive
method.

An alternative measure of the difficulty of a problem is the ground-state
degeneracy, i.e. the number of solutions with minimal waiting-time, as
compared to the total number of solutions, given by $N!$. This is independent
of geometry. Due to the character of the dependence on topology, non-standard
methods are required to compute the average over topology.

This paper is organized as follows: In Section 2, the problem ensemble is
defined. In Section 3, a formalism is introduced for characterizing the
topology. The degeneracy structure as a function of topology is analyzed, and
various energy moments for fixed topology are computed, by averaging over
geometry and/or pairing. In Section 4, a statistical analysis of the detailed
degeneracy structure is considered. In particular, the average ground-state
degeneracy is computed. Section 5 contains our conclusions.

##  2 Unrestricted Crew Scheduling

###  2.1 The Local Problem

A local problem of size $N$ is defined by specifying the times for $N$
arrivals (arr’s) and $N$ departures (dep’s), denoted respectively by
$t^{a}_{i}$ and $t^{d}_{i},\;i=1\ldots N$. The object is simply to find a one-
to-one mapping (a pairing) between the arr’s and dep’s, such that the energy
(or objective function) $E$, given by the total waiting-time, is minimal. In
general this can be done in more than one way, implying a degeneracy of the
ground-state.

The pairing of an arr $A$ with a dep $D$ implies that the crew assigned to $A$
should next be assigned to $D$. The periodicity of the schedule implies that
any arr may be mapped on any dep: if the dep is earlier, it is taken as the
same dep in the next period.

In what follows, we will use the period as the unit of time (and thus energy).
Then the times for the arr’s and dep’s can be restricted to the unit interval.
If an arr $i$ is mapped on a dep $j$, their contribution to the total waiting-
time is given by

$t^{w}_{ij}=\left(t^{d}_{j}-t^{a}_{i}\right)\mbox{ mod }1\;\in[0,1].$ (1)

Thus, whatever the mapping, the energy $E$ is restricted to the interval
$[0,N]$, and between different pairings it can only change by an integer
amount. Pairings yielding the lowest possible energy $E_{0}$ are said to
belong to the ground-state, while a pairing with $E=E_{0}+k,\;k>0$ is said to
belong to $k$:th excited state.

Figure 1: An example of a local flight schedule.

In figure 1, an example of a flight schedule is depicted.

###  2.2 The Random Local Problem

A random local problem is defined by independently choosing the times for $N$
arr’s and $N$ dep’s randomly on the unit period. It is useful to divide the
characterization of a given problem in two parts: its topology and its
geometry. The topology is defined by the relative ordering in the combined
(cyclic) sequence of arr’s and dep’s. The probability is the same for every
distinct topology. The geometry is defined by the sizes of the $2N$ inter-
spaces $x_{i},i=1,\ldots,2N$, into which the period is divided.

For a given problem, the solution space is defined by the set of pairings,
i.e. the $N!$ possible mappings between arr’s and dep’s.

We will be interested mainly in the statistical properties of the following
entities: the ground-state energy (minimal waiting-time) $E_{0}$, the energy
$E$ of a random pairing, and the (integer) difference $D=E-E_{0}$ defining the
excitation energy. We will also be interested in the degeneracies of the
ground-state and the excited states. To that end, we will consider three kinds
of averages: over respectively the topology, the geometry, and the pairing.

The degeneracies of the ground-state and the excited states depend only on the
topology, i.e. the combined ordering of arr’s and dep’s.

For a fixed topology, the excitation energy $D$ is independent of the
geometry, and depends entirely on the pairing. Conversely, the ground-state
energy $E_{0}$ obviously is independent of the pairing, and depends only on
the geometry. Thus, for a fixed topology, $D$ and $E_{O}$ are completely
uncorrelated, in the combined ensemble of random geometries and random
pairings.

##  3 Analysis for Fixed Topology

###  3.1 Characterization of the Topology

A simple way to achieve a ground-state pairing (i.e. solve the problem) for a
given topology is as follows:

  1. 1.

An arr immediately followed by a dep is paired with that dep, and both are
removed from the sequence. Note that the dep could be in the next period.

  2. 2.

The process is continued until all arr’s and dep’s are used.

As an example, consider the sequence
$[A_{1}A_{2}A_{3}D_{1}D_{2}A_{4}D_{3}D_{4}]$, corresponding to the topology of
the schedule in figure 1.

  * •

Pairing $A_{3}$ with $D_{1}$ leaves $[A_{1}A_{2}D_{2}A_{4}D_{3}D_{4}]$.

  * •

Pairing $A_{2}$ with $D_{2}$ leaves $[A_{1}A_{4}D_{3}D_{4}]$.

  * •

Pairing $A_{4}$ with $D_{3}$ leaves $[A_{1}D_{4}]$.

  * •

Pairing $A_{1}$ with $D_{4}$ finishes the process.

A graphical representation of the topology is now defined as follows.

  * •

If necessary, rotate the sequence such that no pairing crosses the interval
border.

  * •

Represent the arr’s and dep’s by equally spaced points on the interval.

  * •

For each pairing in turn, draw a line from the arr to the dep on the lowest
level (one). All previously drawn lines that overlap with the new line are
lifted one level.

Thus a set of lines at different levels are obtained, each line starting at an
arr and ending on a dep. For the example above, the result is shown in fig. 2.
Each line represents a crew waiting for a dep.

Figure 2: Graphical representation of the topology $[AAADDADD]$

It turns out that for the properties of the total energy spectrum, full
knowledge of the topology is not necessary; it suffices to know how many lines
there are at the different levels. Thus let $P_{k}$ be the number of lines at
level $k$ ($k=1,2,3,\ldots$). For the example above, the $P$ sequence is
$[1,2,1]$, i.e. $P_{1}=1,P_{2}=2,P_{3}=1$, while $P_{k}=0$ for $k>3$.

###  3.2 Degeneracy Structure

The ground-state degeneracy is now given by the product of the line-levels,
i.e.

$g_{0}=\prod_{k\geq 1}k^{P_{k}}.$ (2)

This is because each dep must terminate a line, and the number of available
crews is equal to the number of lines alive at that point, which is given by
the line’s level. For the example we get $3\times 2\times 2\times 1=12$,
corresponding to half of the 24 possible pairings.

Naively, the degeneracy of the $m$th excited state can be obtained by adding
$m$ dummy lines, covering the entire interval; they represent extra crews.
Then each dep has $m$ additional crews to choose between, and we get for the
naive degeneracy444For the example above, this precisely corresponds to the
spectrum ${\bf O}|n\\!>=a_{n}|n\\!>$ of the quantum-mechanical operator ${\bf
O}={\bf a}{\bf a}{\bf a}{\bf a}^{\dagger}{\bf a}^{\dagger}{\bf a}{\bf
a}^{\dagger}{\bf a}^{\dagger}$, where ${\bf a}^{\dagger},{\bf a}$ are harmonic
oscillator creation and annihilation operators, satisfying $[{\bf a},{\bf
a}^{\dagger}]=1$, and $|n\\!>\propto\left({\bf a}^{\dagger}\right)^{n}|0\\!>$
is the $n$:th excited state.

$a_{m}=\prod_{k\geq 1}(k+m)^{P_{k}},$ (3)

where we have assumed indistinguishable dummy lines; otherwise we would have
an extra factor $m!$. This defines an infinite sequence.

However, some of the possible pairings will contain permanently grounded crews
(closed lines), or lines extending over more than one period. The contribution
to the naive multiplicity from solutions with $n$ grounded crews and/or
excessive periods depends on the proper degeneracy $n$ steps down. Denoting
the proper degeneracies by $g_{m}$, the relation is

$a_{m}=\sum_{n=0}^{m}g_{m-n}\left(\begin{array}[]{c}N+n\\\
n\end{array}\right),$ (4)

where the last factor is a binomial coefficient. This represents a kind of
renormalization, and can be inverted to yield the proper degeneracies

$g_{m}=\sum_{n=0}^{m}a_{m-n}(-)^{n}\left(\begin{array}[]{c}N+1\\\
n\end{array}\right).$ (5)

This must define a finite sequence, since $g_{m}\geq 0$ and the total number
of pairings, $\sum_{m}g_{m}=N!$, is finite.

From eq. (3), it is obvious that the naive degeneracy $a_{m}$ is an $N$th
degree polynomial in $m$; hence it can be written as

$a_{m}=\sum_{k=0}^{N}c_{k}\left(\begin{array}[]{c}m+k\\\ k\end{array}\right),$
(6)

with some coefficients $c_{k}$. Define the generating functions

$A(x)=\sum_{m}a_{m}x^{m}\;,\;\;G(x)=\sum_{m}g_{m}x^{m}.$ (7)

Due to eq. (5), these are related by $G(x)=(1-x)^{N+1}A(x)$. Then, in terms of
$c_{k}$, we have

$\displaystyle A(x)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{N}c_{k}(1-x)^{-k-1},$ (8) $\displaystyle G(x)$
$\displaystyle=$ $\displaystyle\sum_{k=0}^{N}c_{k}(1-x)^{N-k}.$ (9)

From this, we see that the $g$ sequence is indeed finite: $g(m)=0$ for all
$m>N$. The individual degeneracies $g_{m}$ can be obtained from $G$ and its
derivatives at $x=0$,

$\displaystyle g_{0}$ $\displaystyle=$ $\displaystyle G(0)\;=\;\sum_{k}c_{k},$
(10) $\displaystyle g_{1}$ $\displaystyle=$ $\displaystyle
G^{\prime}(0)\;=\;-\sum_{k=0}^{N}(N-k)c_{k},$

etc.

###  3.3 The Excitation Energy for a Random Pairing

Conversely, the moments over $g_{m}$ are readily obtained from $G$ and its
derivatives at $x=1$,

$\displaystyle\sum_{m}g_{m}$ $\displaystyle=$ $\displaystyle G(1)\;=\;c_{N},$
(11) $\displaystyle\sum_{m}mg_{m}$ $\displaystyle=$ $\displaystyle
G^{\prime}(1)\;=\;-c_{N-1},$ $\displaystyle\sum_{m}m(m-1)g_{m}$
$\displaystyle=$ $\displaystyle G^{\prime\prime}(1)\;=\;2c_{N-2},$

etc. In order to relate $c_{k}$ to $P_{k}$, we can express the coefficients of
the polynomial $a_{m}$ in two different ways. From eq. (3) we get for the
leading coefficients

$a_{m}=m^{N}+m^{N-1}\sum_{k}kP_{k}+\frac{m^{N-2}}{2}\left\\{\left(\sum_{k}kP_{k}\right)^{2}-\sum_{k}k^{2}P_{k}\right\\}+\ldots,$
(12)

while eq. (6) gives, upon expanding the binomial coefficients,

$\displaystyle a_{m}$ $\displaystyle=$
$\displaystyle\frac{m^{N}}{N!}c_{N}\;+\;\frac{m^{N-1}}{N!}\left\\{\frac{1}{2}N(N+1)c_{N}+Nc_{N-1}\right\\}$
$\displaystyle+$
$\displaystyle\frac{m^{N-2}}{N!}\left\\{\frac{1}{24}N(N-1)(N+1)(3N+2)c_{N}+\frac{1}{2}N^{2}(N-1)c_{N-1}+N(N-1)c_{N-2}\right\\}+\ldots$

From this we obtain the relations

$\displaystyle c_{N}$ $\displaystyle=$ $\displaystyle N!,$ (14) $\displaystyle
c_{N-1}$ $\displaystyle=$ $\displaystyle
N!\left\\{\frac{1}{N}\sum_{k}kP_{k}-\frac{N+1}{2}\right\\},$ $\displaystyle
c_{N-2}$ $\displaystyle=$ $\displaystyle
N!\left\\{\frac{1}{2N(N-1)}\left(\sum_{k}kP_{k}\right)^{2}-\frac{1}{2N(N-1)}\sum_{k}k^{2}P_{k}-\frac{1}{2}\sum_{k}kP_{k}+\frac{(N+1)(3N-2)}{24}\right\\},$

etc. This gives (as it should) $\sum_{m}g_{m}=N!$, the number of possible
pairings. In a given topology, the fraction of pairings having an energy $D=m$
steps above $E_{0}$ is given by $g_{m}/N!$. Thus, the first few moments of the
excitation energy $D$ for a random pairing are:

$\displaystyle<D>_{p}$ $\displaystyle=$
$\displaystyle\frac{1}{N!}\sum_{m}mg_{m}\;=\;\frac{N+1}{2}-\frac{1}{N}\sum_{k}kP_{k},$
$\displaystyle<D^{2}>_{p}$ $\displaystyle=$
$\displaystyle\frac{1}{N!}\sum_{m}m^{2}g_{m}$ $\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\left(\sum_{k}kP_{k}\right)^{2}-\frac{1}{N(N-1)}\sum_{k}k^{2}P_{k}-\frac{N+1}{N}\sum_{k}kP_{k}+\frac{(N+1)(3N+4)}{12},$

where $<>_{p}$ denotes the average over pairings for a fixed topology.

###  3.4 The Ground-State Energy in a Random Geometry

For a given topology, the ground-state energy depends on the geometry, and is
given by

$E_{0}=\sum_{i=1}^{2N}k_{i}x_{i},$ (16)

where $x_{i}$ is the length of the $i$th sub-interval, and $k_{i}$ the number
of lines in that interval.

To perform averages over a random geometry, we need to analyze the
distribution of the intervals $x_{i}$. Independently of the topology, they
obey the distribution

$dP=(2N-1)!\;\delta\left(\sum_{i}x_{i}-1\right)\prod_{i=1}^{2N}\left(\Theta(x_{i})dx_{i}\right).$
(17)

For a single $x_{i}$, this implies the distribution

$f(x_{i})=(2N-1)\left(1-x_{i}\right)^{2N-2}\;,\;0\leq x_{i}\leq 1.$ (18)

Using the identity $\sum_{i}x_{i}=1$, and the permutation symmetry between
different $x_{i}$, we have

$\displaystyle<x_{i}>$ $\displaystyle=$ $\displaystyle\frac{1}{2N},$ (19)
$\displaystyle<x_{i}x_{j}>$ $\displaystyle=$
$\displaystyle\frac{1+\delta_{ij}}{2N(2N+1)}.$ (20)

In what follows, we will also need the number of intervals with $k$ lines, to
be denoted by $Q_{k}$; in terms of $P_{k}$ (defining $P_{0}=0$) it is simply

$Q_{k}=P_{k}+P_{k+1},\;k\geq 0,$ (21)

since adding a line on level $k$ implies adding two intervals, with
respectively $k$ and $k-1$ lines. Obviously, we have $\sum_{k}Q_{k}=2N$.

Now we are ready to compute the first few moments of $E_{0}$, yielding

$\displaystyle<E_{0}>_{g}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{2N}k_{i}<x_{i}>_{g}=\frac{1}{2N}\sum_{i}k_{i}=\frac{1}{2N}\sum_{k}kQ_{k}=\frac{1}{N}\sum_{k}kP_{k}-\frac{1}{2},$
$\displaystyle<E_{0}^{2}>_{g}$ $\displaystyle=$
$\displaystyle\sum_{ij}k_{i}k_{j}<x_{i}x_{j}>=\frac{1}{2N(2N+1)}\sum_{ij}k_{i}k_{j}(1+\delta_{ij})$
$\displaystyle=$
$\displaystyle\frac{1}{2N(2N+1)}\left(\sum_{kl}klQ_{k}Q_{l}+\sum_{k}k^{2}Q_{k}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2N(2N+1)}\left(4\left(\sum_{k}kP_{k}\right)^{2}+2\sum_{k}k^{2}P_{k}-2(2N+1)\sum_{k}kP_{k}+N(N+1)\right),$

where $<>_{g}$ denotes average over the geometry for fixed topology.

###  3.5 The Total Energy

The total energy $E$ is given by the ground-state and excitation energies,
$E=E_{0}+D$. Averaging over both geometry and pairing, we obtain the
surprisingly simple result

$<E>_{gp}\,=\,<E_{0}>_{g}+<D>_{p}\,=\,\frac{N}{2},$ (23)

independently of topology.

Since we already have the second moments of $E_{0}$ and $D$, we only need
$<E_{0}D>$ to compute the second moment of $E$. Using the fact that $E_{0}$
and $D$ are statistically independent for a fixed topology, we get

$<E_{0}D>_{gp}\,=\,<E_{0}>_{g}<D>_{p}\,=\,-\frac{1}{N^{2}}\left(\sum_{k}kP_{k}\right)^{2}+\frac{N+2}{2N}\sum_{k}kP_{k}-\frac{N+1}{4}.$
(24)

This implies

$<E^{2}>_{gp}=\frac{2(N^{2}+N+1)}{N^{2}(N-1)(2N+1)}\left(\sum_{k}kP_{k}\right)^{2}-\frac{3(N+1)}{N(N-1)(2N+1)}\sum_{k}k^{2}P_{k}+\frac{(N+1)(6N^{2}-N+4)}{12(2N+1)}.$
(25)

##  4 Topology Statistics

###  4.1 The Topology Ensemble

We now want to perform averages also over the topology, for the moments of the
ground-state energy $E_{0}$, the excitation energy $D$, and the total energy
$E$.

$N$ | Pattern | Prob. | $\\{P_{k}\\}$ | $\\{Q_{k}\\}$ | $\\{a_{i}\\}$ | $\\{g_{i}\\}$ | $\langle E_{0}\rangle$ | $\langle D\rangle$
---|---|---|---|---|---|---|---|---
1 | AD | 1 | [1] | [1,1] | 1,2,3,… | (1) | 1/2 | 0
2 | AADD | 2/3 | [1,1] | [1,2,1] | 2,6,12,… | (2) | 2/2 | 0
2 | ADAD | 1/3 | [2] | [2,2] | 1,4,9,… | (1,1) | 1/2 | 1/2
3 | AAADDD | 3/10 | [1,1,1] | [1,2,2,1] | 6,24,60,… | (6) | 9/6 | 0
3 | AADADD | 3/10 | [1,2] | [1,3,2] | 4,18,48,… | (4,2) | 7/6 | 1/3
3 | AADDAD | 3/10 | [2,1] | [2,3,1] | 2,12,36,… | (2,4) | 5/6 | 2/3
3 | ADADAD | 1/10 | [3] | [3,3] | 1,8,27,… | (1,4,1) | 3/6 | 3/3
4 | AAAADDDD | 4/35 | [1,1,1,1] | [1,2,2,2,1] | 24,120,360,… | (24) | 8/4 | 0
4 | AAADADDD | 4/35 | [1,1,2] | [1,2,3,2] | 18,96,300,… | (18,6) | 7/4 | 1/4
4 4 |  AAADDADD AADAADDD | $\left.\parbox{19.91692pt}{4/35 4/35}\right\\}$ | [1,2,1] | [1,3,3,1] | 12,72,240,… | (12,12) | 6/4 | 2/4
4 | AAADDDAD | 4/35 | [2,1,1] | [2,3,2,1] | 6,48,180,… | (6,18) | 5/4 | 3/4
4 | AADADADD | 4/35 | [1,3] | [1,4,3] | 8,54,192,… | (8,14,2) | 5/4 | 3/4
4 4 |  AADDAADD AADADDAD | $\left.\parbox{19.91692pt}{2/35 4/35}\right\\}$ | [2,2] | [2,4,2] | 4,36,144,… | (4,16,4) | 4/4 | 4/4
4 | AADDADAD | 4/35 | [3,1] | [3,4.1] | 2,24,108,… | (2,14,8) | 3/4 | 5/4
4 | ADADADAD | 1/35 | [4] | [4,4] | 1,16,81,256,… | (1,11,11,1) | 2/4 | 6/4
Table 1: Inequivalent topologies for $N\leq 4$. For each topology, the third
column gives the probability of occurrence, while the next two give its
$P_{k}$ and $Q_{k}$ sequence, respectively. In the following two columns the
$a_{i}$ and $g_{i}$ sequences are given, whereas the last two columns give the
average ground-state energy and excitation, respectively. Note that different
topologies might have identical characteristics in terms of $P_{k}$, etc.

In table 1 a complete list of topologies for $N$ up to four is given, along
with various characteristics. Note the similarity in symmetry-properties
between the $P$, $Q$ and $g$ sequences.

To compute $<P_{k}>$, we must analyze the number of ways $m[P]$ to obtain a
given $P$ sequence.

  * •

The lowest level lines define $P_{1}$ non-overlapping subgroups; these are
cyclically indistinguishable, so the naive multiplicity must be divided by
$P_{1}$.

  * •

For the $P_{k}$ indistinguishable lines on a higher level $k$, there are
$P_{k-1}$ possible lines on the previous level to put them above. This can be
done in $\left(\begin{array}[]{c}P_{k}+P_{k-1}-1\\\ P_{k}\end{array}\right)$
ways.

  * •

In addition, there are $2N$ cyclic rotations of the $AD$ sequence that yield
the same topology.

  * •

The $P$ sequence must sum up to $N$, and $P_{1}$ must not vanish.

This gives for the multiplicity $m[P]$

$m[P]=\frac{2N}{P_{1}}\delta_{N-\sum
P_{k}}\prod_{k=1}^{\infty}\left(\begin{array}[]{c}P_{k}+P_{k+1}-1\\\
P_{k+1}\end{array}\right).$ (26)

The total number of possible topologies for a fixed $N$ is simply the number
of possible $AD$ sequences, i.e. the number of distinct orderings of a
sequence of $N$ $A$’s and $N$ $D$’s. This is given by
$\left(\begin{array}[]{c}2N\\\ N\end{array}\right)$. This should match the
total multiplicity $M_{N}\equiv\sum_{[P]}m([P])$, which can be recursively
obtained from $m[P]$ by expanding the Kronecker delta in terms of a complex
integral along a small contour $C$ around the origin:

$\delta_{N-\sum P_{k}}=\oint_{C}\frac{dz}{2\pi iz^{N+1}}\prod z^{P_{k}}.$ (27)

By iteratively using the identity

$z^{a}\sum_{b=0}^{\infty}\left(\begin{array}[]{c}a+b-1\\\
a-1\end{array}\right)w^{b}=z^{a}\sum_{b=0}^{\infty}\left(\begin{array}[]{c}-a\\\
b\end{array}\right)(-w)^{b}=\left(\frac{z}{1-w}\right)^{a}.$ (28)

when summing $m[P]\prod z^{P_{k}}$ over the $P_{k}$ in order of decreasing
$k$, we obtain convergence of the power factor to $\omega^{P_{k}}$, with
$\omega$ given by

$\omega=\frac{z}{1-\omega}\;\Rightarrow\;\omega=\frac{1}{2}\left(1-\sqrt{1-4z}\right).$
(29)

Thus, the total multiplicity is given by

$M_{N}=2N\oint_{C}\frac{dz\;z^{-N}}{2\pi
iz}\sum_{P_{1}}\frac{\omega(z)^{P_{1}}}{P_{1}}=2N\oint_{C}\frac{dz\;z^{-N}}{2\pi
iz}\left(-log(1-\omega(z))\right)=\left(\begin{array}[]{c}2N\\\
N\end{array}\right),$ (30)

obtained by extracting the $z^{N}$ coefficient of $-\log(1-\omega(z))$ as
$\frac{1}{2N}\left(\begin{array}[]{c}2N\\\ N\end{array}\right)$, by doing the
integral in terms of $w$ using $z=w(1-w)$.

###  4.2 Distribution of $P_{1}$

The number of non-overlapping groups of lines is given by the number of lines
at level 1, i.e. $P_{1}$. In fact, $w(z)$ can be seen as a generating
function,

$w(z)=\sum_{n=1}^{\infty}m^{(1)}_{n}z^{n}=-\frac{1}{2}\sum_{n=1}^{\infty}(4z)^{n}\left(\begin{array}[]{c}n-\frac{3}{2}\\\
n\end{array}\right)=\frac{1}{2}\sum_{n=1}^{\infty}(2z)^{n}\frac{(2n-3)!!}{n!},$
(31)

for the number of ways ($m^{(1)}_{n}$) to arrange $n$ lines in such a group.
The distribution of $P_{1}$ is easy to compute by slightly modifying eq. (30):
skipping the summation over $P_{1}$ and dividing by $M_{N}$ yields

$\mbox{Prob}(P_{1})=\frac{2N}{\left(\begin{array}[]{c}2N\\\
N\end{array}\right)}\oint_{C}\frac{dz\;z^{-N}}{2\pi
iz}\frac{\omega(z)^{P_{1}}}{P_{1}}=\frac{N!(2N-1-P_{1})!}{(2N-1)!(N-P_{1})!}$
(32)

for $P_{1}>0$. For large $P_{1}\ll N$, this is close to $2^{-P_{1}}$.

The grouping of lines corresponds to a grouping of the flights, with a
matching number of arr’s and dep’s in each group. In a ground-state pairing,
flights in different groups are never paired, and within a group, the arr has
to precede its paired dep. This fact can be used to simplify also restricted
problems.

###  4.3 Moments of $P_{k}$

Similarly, $<P_{k}>$ can be obtained from realizing that inserting a factor
$P_{k}$ in the sum over $P_{k}$ gives a factor $P_{k-1}\omega/(1-\omega)$.
This yields in the end an extra factor $P_{1}(\omega/(1-\omega))^{k-1}$ in the
$P_{1}$ sum:

$\displaystyle M_{N}<P_{k}>$ $\displaystyle=$ $\displaystyle
2N\oint_{C}\frac{dz}{2\pi
iz^{N+1}}\left(\frac{\omega}{1-\omega}\right)^{k-1}\sum_{P_{1}>1}\omega^{P_{1}}$
(33) $\displaystyle=$ $\displaystyle 2N\oint_{C}\frac{dz}{2\pi
iz^{N+1}}\left(\frac{\omega}{1-\omega}\right)^{k}$ (34) $\displaystyle=$
$\displaystyle 2N\oint_{C}\frac{d\omega(1-2\omega)}{2\pi
i\omega^{N-k+1}(1-\omega)^{N+k+1}},$ (35)

where the last expression is a reformulation in terms of a loop integral
around $\omega=0$. In a similar way $<P_{k}P_{l}>$ etc. can be computed. We
obtain

$\displaystyle<P_{k}>$ $\displaystyle=$
$\displaystyle\frac{2N}{M_{N}}\left\\{\left(\begin{array}[]{c}2N-1\\\
N-k\end{array}\right)-\left(\begin{array}[]{c}2N-1\\\
N-k-1\end{array}\right)\right\\}=\frac{2k}{M_{N}}\left(\begin{array}[]{c}2N\\\
N-k\end{array}\right),$ (42) $\displaystyle<P_{k}P_{m}>$ $\displaystyle=$
$\displaystyle\frac{2N}{M_{N}}\left\\{\left(\begin{array}[]{c}2N\\\
N-m\end{array}\right)-\left(\begin{array}[]{c}2N\\\
N-k-m\end{array}\right)\right\\},\;k\leq m.$ (47)

From these expressions, we can derive the following particular averages,
needed to compute the various energy moments over the topology:

$\displaystyle\left\langle\sum_{k}kP_{k}\right\rangle$ $\displaystyle=$
$\displaystyle\frac{N}{2}4^{N}\left(\begin{array}[]{c}2N\\\
N\end{array}\right)^{-1}\approx\frac{N}{2}\sqrt{\pi
N}\left(1+\frac{1}{8N}+\ldots\right),$ (50)
$\displaystyle\left\langle\sum_{k}k^{2}P_{k}\right\rangle$ $\displaystyle=$
$\displaystyle N^{2},$ (51)
$\displaystyle\left\langle\left(\sum_{k}kP_{k}\right)^{2}\right\rangle$
$\displaystyle=$ $\displaystyle\frac{N^{2}(5N+1)}{6},$

where the approximate form in the top equation is valid for large $N$.

###  4.4 Full Energy Averages

We are now ready to compute the final averages also over the topology.
Inserting the results of eqs. (51) into eqs. (3.4), we obtain for the moments
of the ground-state energy $E_{0}$ of a random problem:

$\displaystyle<E_{0}>$ $\displaystyle=$
$\displaystyle\frac{1}{2}4^{N}\left(\begin{array}[]{c}2N\\\
N\end{array}\right)^{-1}-\frac{1}{2}\;\approx\;\frac{1}{2}\sqrt{\pi
N}-\frac{1}{2},$ (54) $\displaystyle<E_{0}^{2}>$ $\displaystyle=$
$\displaystyle\frac{5N}{6}-\frac{1}{2}4^{N}\left(\begin{array}[]{c}2N\\\
N\end{array}\right)^{-1}+\frac{1}{2}\;\approx\;\frac{5N}{6},$ (57)
$\displaystyle<E_{0}^{2}>_{c}$ $\displaystyle\approx$
$\displaystyle\frac{(10-3\pi)N}{12},$

where $\langle ab\rangle_{c}=\langle ab\rangle-\langle a\rangle\langle
b\rangle$ is the connected moment; the approximate forms are valid for large
$N$.

Similarly, from eqs. (3.3) we get for the moments of the excitation energy
$D$, for a random pairing in a random topology,

$\displaystyle<D>$ $\displaystyle=$
$\displaystyle\frac{N+1}{2}-\frac{1}{2}4^{N}\left(\begin{array}[]{c}2N\\\
N\end{array}\right)^{-1}\;\approx\;\frac{N}{2}-\frac{1}{2}\sqrt{\pi
N}+\frac{1}{2},$ (60) $\displaystyle<D^{2}>$ $\displaystyle=$
$\displaystyle\frac{N^{2}}{4}+\frac{17N}{12}+\frac{1}{3}-\frac{N+1}{2}4^{N}\left(\begin{array}[]{c}2N\\\
N\end{array}\right)^{-1}$ $\displaystyle\approx$
$\displaystyle\frac{N^{2}}{4}-\frac{N}{2}\sqrt{\pi N}+\frac{17N}{12},$
$\displaystyle<D^{2}>_{c}$ $\displaystyle\approx$
$\displaystyle\frac{(11-3\pi)N}{12}.$ (64)

Averaging the combined moment, eq. (24), over topology yields

$\displaystyle<E_{0}D>$ $\displaystyle=$
$\displaystyle\frac{N+2}{4}4^{N}\left(\begin{array}[]{c}2N\\\
N\end{array}\right)^{-1}-\frac{13N+5}{12}\;\approx\;\frac{N}{4}\sqrt{\pi
N}-\frac{13N}{12},$ (67) $\displaystyle<E_{0}D>_{c}$ $\displaystyle=$
$\displaystyle-\frac{5N+1}{6}+\frac{1}{4}16^{N}\left(\begin{array}[]{c}2N\\\
N\end{array}\right)^{-2}\;\approx\;-\frac{(10-3\pi)N}{12}.$ (70)

Combining this with the moments of $E_{0}$ and $D$, we get for the moments of
the total energy, $E=E0+D$, the following simple results

$\displaystyle<E>$ $\displaystyle=$ $\displaystyle\frac{N}{2},$ (71)
$\displaystyle<E^{2}>$ $\displaystyle=$
$\displaystyle\frac{N^{2}}{4}+\frac{N}{12},$ (72) $\displaystyle<E^{2}>_{c}$
$\displaystyle=$ $\displaystyle\frac{N}{12},$ (73)

which can be understood by noting that a random pairing in a random topology
corresponds to a set of $N$ lines with random endpoints. Then each line has a
uniform length distribution between 0 and 1, and $E$ is their total length.

Computing the corresponding standard deviations, we have for a typical random
problem (and a random pairing, for $D$ and $E$) at large $N$,

$\displaystyle E_{0}$ $\displaystyle\sim$ $\displaystyle\frac{1}{2}\sqrt{\pi
N}\;\pm\;\frac{1}{2}\sqrt{\left(\frac{10}{3}-\pi\right)N},$ (74)
$\displaystyle D$ $\displaystyle\sim$
$\displaystyle\frac{N}{2}-\frac{1}{2}\sqrt{\pi
N}\;\pm\;\frac{1}{2}\sqrt{\left(\frac{11}{3}-\pi\right)N},$ (75)
$\displaystyle E$ $\displaystyle\sim$
$\displaystyle\frac{N}{2}\;\pm\;\frac{1}{2}\sqrt{\frac{N}{3}},$ (76)

and we see that $E$ scales as $\sqrt{N}$, while $E$ and $D$ scale as $N$,
while the standard deviation scales as $\sqrt{N}$ in each case.

Of interest are also the correlations at large $N$, given to order $N$ by

$\displaystyle<E_{0}D>_{c}$ $\displaystyle\sim$
$\displaystyle-\frac{N}{12}(10-3\pi),$ $\displaystyle<E_{0}E>_{c}$
$\displaystyle\sim$ $\displaystyle 0,$ (77) $\displaystyle<DE>_{c}$
$\displaystyle\sim$ $\displaystyle\frac{N}{12},$

indicating that $E$ and $E_{0}$ become uncorrelated for a random pairing of a
large random problem.

###  4.5 Statistics for Individual Degeneracies

An interesting but more difficult thing to compute is the average fraction
$\gamma_{n}=<g_{n}>/N!$ of pairings having a given excitation energy $D=n$, in
particular the ground-state fraction $\gamma_{0}$, which gives the average
probability of hitting a ground-state by chance. Since $g_{n}$ is simply
related to $a_{n}$, we will start by considering

$\alpha_{n}=\frac{\left(\begin{array}[]{c}2N\\\
N\end{array}\right)}{2N}<a_{n}>,$ (78)

in terms of which $\gamma_{n}$ can be expressed as

$\gamma_{n}=\frac{N!}{(2N-1)!}\sum_{m}(-)^{m}\left(\begin{array}[]{c}N+1\\\
m\end{array}\right)\alpha_{n-m}.$ (79)

We then have

$\alpha_{n}=\frac{1}{2N}\sum_{[P]}m[P]\prod_{k}(n+k)^{P_{k}}=\sum_{[P]}\frac{1}{P_{1}}\delta_{N-\sum
P_{k}}\prod_{k\geq 1}(n+k)^{P_{k}}\left(\begin{array}[]{c}P_{k}+P_{k+1}-1\\\
P_{k+1}\end{array}\right).$ (80)

A generating function for the $N$-dependence of $\alpha_{n}$ is then

$A_{n}(z)\equiv\sum_{N}z^{N}\alpha_{n}(N)=\sum_{[P]}\frac{1}{P_{1}}\prod_{k\geq
1}[z(n+k)]^{P_{k}}\left(\begin{array}[]{c}P_{k}+P_{k+1}-1\\\
P_{k+1}\end{array}\right).$ (81)

Again, starting the $P_{k}$ summation at a large $k_{0}$ and proceeding in
order of decreasing $k$, gives convergence of the full sum in the limit
$k_{0}\to\infty$. Denoting by $\omega_{k+n-1}^{P_{k}}$ the result of summing
above a certain $k$, we have by eq. (28),

$\omega_{k-1}=\frac{kz}{1-w_{k}},$ (82)

and the final result is

$A_{n}(z)=\sum_{P_{1}}\frac{1}{P_{1}}\omega_{n}^{P_{1}}=-\log\left(1-\omega_{n}(z)\right).$
(83)

The recursion relation (82) can now be linearized by assuming
$\omega_{k}=p_{k}/q_{k}$, which can be solved e.g. by

$\displaystyle q_{k-1}$ $\displaystyle=$ $\displaystyle q_{k}-p_{k},$ (84)
$\displaystyle p_{k-1}$ $\displaystyle=$ $\displaystyle kzq_{k},$ (85)

which, upon eliminating $p_{k}$ gives

$kzq_{k}-q_{k-1}+q_{k-2}=0.$ (86)

By partial integration, it is simple to prove that the following sequence of
integrals solves the recursion relation (86) for $q_{k}$,

$q_{k}=\mbox{Im}\int_{0}^{\infty}\frac{(-iy)^{k}}{k!}\exp\left(-z\frac{y^{2}}{2}+iy\right)dy\;,\;k\geq
0,\;z>0,$ (87)

which can be extended to negative $k$ by recursion. The $q_{k}$ can be
expanded in $z$ as

$q_{k}=\left\\{\begin{array}[]{l}{}\sum_{m=0}^{\infty}\left(\begin{array}[]{c}2m+k\\\
2m\end{array}\right)(2m-1)!!z^{m}\;,\;k\geq 0,\\\
\sum_{m=0}^{[-(k+1)/2]}\left(\begin{array}[]{c}-k-1\\\
2m\end{array}\right)(2m-1)!!z^{m}\;,\;k<0,\end{array}\right.$ (88)

where $[\,]$ denotes integer part. Note that for non-negative $k$ the series
is an asymptotic one, while for negative $k$ it is finite. In particular, we
have $q_{-1}=1$. An independent solution to (86) is given by

$\hat{q}_{k}(z)=\frac{z^{-k}}{k!}q_{-k-1}(-z),$ (89)

but this solution is irrelevant, having the wrong large-$k$ behaviour.

In terms of $q_{k}$ we now have

$\displaystyle 1-\omega_{k}$ $\displaystyle=$ $\displaystyle q_{k-1}/q_{k},$
(90) $\displaystyle A_{n}$ $\displaystyle=$
$\displaystyle\log\left(q_{n}\right)-\log\left(q_{n-1}\right),$ (91)

which should then be expanded in powers of $z$ to yield $\alpha_{n}(N)$ as the
coefficient of $z^{N}$. In particular, we have $A_{0}=\log(q_{0})$, where
$q_{0}=1+z+3z^{2}+15z^{3}+\ldots$.

In table 2, results are displayed for the average degeneracy of the two lowest
states, based on an expansion of $A_{0}$ and $A_{1}$ in powers of $z$.

$N$ | $K_{N}=\frac{1}{2}\left(\begin{array}[]{c}2N\\\ N\end{array}\right)$ | $K_{N}<g_{0}>$ | $K_{N}<g_{1}>$ | $<g_{0}>$ | $<g_{1}>$ | $\gamma_{0}$ | $\gamma_{1}$
---|---|---|---|---|---|---|---
1 | 1 | 1 | 0 | 1.00000 | 0.00000 | 1.000000 | 0.000000
2 | 3 | 5 | 1 | 1.66667 | 0.33333 | 0.833333 | 0.166667
3 | 10 | 37 | 22 | 3.70000 | 2.20000 | 0.616667 | 0.366667
4 | 35 | 353 | 411 | 10.0857 | 11.7429 | 0.420238 | 0.489286
5 | 126 | 4081 | 7676 | 32.3889 | 60.9206 | 0.269907 | 0.507672
6 | 462 | 55205 | 149741 | 119.491 | 324.115 | 0.165960 | 0.450159
7 | 1716 | 854197 |  3.09875$\times 10^{6}$ | 497.784 | 1805.80 | 0.098767 | 0.358294
8 | 6435 |  1.4876$\times 10^{7}$ |  6.84187$\times 10^{7}$ | 2311.74 | 10632.3 | 0.057335 | 0.263697
9 | 24310 |  2.88019$\times 10^{8}$ |  1.61447$\times 10^{9}$ | 11847.7 | 66411.7 | 0.032649 | 0.183013
10 | 92378 |  6.13891$\times 10^{9}$ |  4.07031$\times 10^{10}$ | 66454.3 | 440614. | 0.018313 | 0.121421
11 | 352716 |  1.42882$\times 10^{11}$ |  1.09496$\times 10^{12}$ | 405092. |  3.10436$\times 10^{6}$ | 0.010148 | 0.077771
12 |  1.35208$\times 10^{6}$ |  3.60668$\times 10^{12}$ |  3.13708$\times 10^{13}$ |  2.66751$\times 10^{6}$ |  2.32019$\times 10^{7}$ | 0.005569 | 0.048438
13 |  5.20030$\times 10^{6}$ |  9.81584$\times 10^{13}$ |  9.55147$\times 10^{14}$ |  1.88755$\times 10^{7}$ |  1.83672$\times 10^{8}$ | 0.003031 | 0.029496
14 |  2.00583$\times 10^{7}$ |  2.86562$\times 10^{15}$ |  3.08337$\times 10^{16}$ |  1.42865$\times 10^{8}$ |  1.5372$\times 10^{9}$ | 0.001639 | 0.017633
Table 2: Results for the average degeneracy of the lowest energy-states for
various system sizes $N$. The second column gives an integer normalization
factor $K_{N}$. Dividing the integers in the next two columns by $K_{N}$
yields the average number of ground-states $<g_{0}>$ and first excited states
$<g_{1}>$, respectively. Dividing these by $N!$ yields $\gamma_{0}$ and
$\gamma_{1}$.

It is easy to check for small $N$, using table 1, that $g_{m}/\sum g_{k}$ for
$m=0,1$, averaged over topologies with the proper probabilities, indeed agrees
with $\gamma_{m}$ of table 2, obtained from the expansion of $A_{k}$.

The result for $\gamma_{0}$ strongly indicates an asymptotic behaviour of
$\gamma_{0}\sim 2N2^{-N}$. This corresponds to an exponential decrease with
$N$ of the average probability for a random pairing to hit a ground-state.
However, the average number of ground-states grows faster than exponentially:
$<g_{0}>\sim 2N\;N!\;2^{-N}$.

This abundance of ground-states indicates that a corresponding restricted
problem might well have a solution with a locally minimal waiting-time, if the
restrictions are not too severe; this is used in ref. [1] to simplify the
solution of a set of restricted crew scheduling problems.

By reducing the state-space of a restricted problem to the set of ground-
states of the corresponding unrestricted problem, the average information
gained at each airport is given by $\log\gamma_{0}$, which for a large airport
roughly yields $N\log 2$. Summing this over several airports yields a total
information gain scaling as $N_{f}\log 2$, with $N_{f}$ the total number of
flights. Partly, this is due to a grouping of flights, which contributes an
average information gain of $N_{a}\log 2$, with $N_{a}$ the number of
airports.

##  5 Conclusions

We have performed a statistical analysis of an ensemble of random unrestricted
local crew scheduling problems, formulated in terms of mapping arrivals onto
departures at a single airport so as to minimize waiting-time.

For the ground-state energy $E_{0}$ of a large random problem, we find that
both the average and the fluctuations scale like $\sqrt{N}$. For a random
pairing of such a problem, on the other hand, the excitation $D$ and the total
energy $E=E_{0}+D$ both grow linearly with $N$, with fluctuations scaling like
$\sqrt{N}$.

The individual degeneracies of the lowest energy states for random problems
are such that the average probability for hitting an optimal pairing by chance
decreases like $2N2^{-N}$ for large $N$. Since the total number of pairings
grows like $N!$, the average number of ground-states grows very fast with
system size.

The results and the methods of analysis described in this paper are useful
when designing efficient algorithms for the crew scheduling problems with
global restrictions, by providing means for estimating the difficulty of a
given problem, and for understanding and simplifying its structure.

The optimal crew waiting-time for a restricted problem is bounded from below
by the ground-state energy of the corresponding unrestricted problem, which is
useful for gauging algorithmic performance for problems of realistic size. Due
to a faster than exponential growth of the number of ground-states with
problem size, this bound is often saturated. This can be used to simplify a
restricted problem: By insisting on the local ground-state energy, airports
can be split into several parts. For a large random problem, this results on
the average in a reduction of state-space size by a factor of two for each
airport.

Some of the calculations in this paper are based on novel methods of a general
nature, that may have applications also in other contexts with a similar
topological structure.

## References

  * [1] M. Lagerholm, C. Peterson and B. Söderberg, “Airline Crew Scheduling using Potts Mean Field Techniques”, LU TP 97-10 (submitted to Operations Research).