snr short = SNR, long = signal-to-noise ratio, sar short = SAR, long =
synthetic aperture radar, insar short = InSAR, long = interferometric
synthetic aperture radar, minlp short = MINLP, long = mixed-integer non-linear
program, long-plural-form = mixed-integer non-linear programs uav short = UAV,
long = unmanned aerial vehicle, long-plural-form = unmanned aerial vehicles
fdma short = FDMA, long = frequency-division multiple-access, 3d short = 3D,
long = three-dimensional, 2d short = 2D, long = two-dimensional, dem short =
DEM, long = digital elevation model, gs short = GS, long = ground station,
long-plural-form = ground stations los short = LOS, long = line-of-sight, sca
short = SCA, long = successive convex approximation, ao short = AO, long =
alternating optimization, aoi short = AoI, long = area of interest, nesz short
= NESZ, long = noise equivalent sigma zero, wrt short = w.r.t., long = with
respect to, rhs short = r.h.s., long = right-hand side, lhs short = l.h.s.,
long = left-hand side, hoa short = HoA, long = height of ambiguity,

# Robust Trajectory and Resource Optimization for Communication-assisted UAV
SAR Sensing
††thanks: Mohamed-Amine Lahmeri, Walid R. Ghanem, and Robert Schober are with
the Institute for Digital Communications, Friedrich-Alexander University
Erlangen-Nürnberg (FAU), Germany (emails:{amine.lahmeri, walid.ghanem,
robert.schober}@fau.de). Christina Bonfert is with the Institute of Microwave
Engineering, Ulm, Germany (email: christina.knill@uni-ulm.de). This paper was
presented in part at the IEEE Global Communications Conference, Rio De
Janeiro, Brazil, Dec. 2022 [1]. This work was supported in part by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) GRK 2680 –
Project-ID 437847244.

Mohamed-Amine Lahmeri, Walid R. Ghanem, Christina Bonfert, and Robert Schober

###### Abstract

In this paper, we investigate joint 3-dimensional (3D) trajectory planning and
resource allocation for rotary-wing unmanned aerial vehicle (UAV) synthetic
aperture radar (SAR) sensing. To support emerging real-time SAR applications
and enable live mission control, we incorporate real-time communication with a
ground station (GS). The UAV’s main mission is the mapping of large areas of
interest (AoIs) using an onboard SAR system and transferring the unprocessed
raw radar data to the ground in real time. We propose a robust trajectory and
resource allocation design that takes into account random UAV trajectory
deviations. To this end, we model the UAV trajectory deviations and study
their effect on the radar coverage. Then, we formulate a robust non-convex
mixed-integer non-linear program (MINLP) such that the UAV 3D trajectory and
resources are jointly optimized for maximization of the radar ground coverage.
A low-complexity sub-optimal solution for the formulated problem is presented.
Furthermore, to assess the performance of the sub-optimal algorithm, we derive
an upper bound on the optimal solution based on monotonic optimization theory.
Simulation results show that the proposed sub-optimal algorithm achieves
close-to-optimal performance and not only outperforms several benchmark
schemes but is also robust with respect to UAV trajectory deviations.

###### Index Terms:

Optimization, communication, sensing, unmanned aerial vehicles, synthetic
aperture radar, successive convex approximations, monotonic optimization.

##  I INTRODUCTION

Synthetic aperture radar (sar) is a moving radar system that uses the flight
path of a moving platform to simulate a large aperture [1, 2]. In classical
radar systems, the spatial resolution depends on the ratio of the wavelength
of the radar signal to the physical length of the antenna. This means that
spatial resolution can be improved by using large antenna arrays, which are
however often not practical. sar overcomes this hurdle by using the motion of
the moving platform to form a synthetic aperture where observations in
different time slots, corresponding to different positions of the trajectory,
are jointly processed to create the final image of the target area [3, 4]. sar
can provide high-resolution 2d and even 3d images [5, 6]. The benefits of sar
also include its resilience to harsh weather conditions and its ability to
operate at day and night. In addition, sar is non-invasive, which is useful
for applications where penetration or direct contact with the aoi is not
possible or not desirable, such as imaging of glaciers and wetlands [5].
Thanks to these features, sar sensing has been employed in a variety of
applications, such as environmental monitoring, agriculture, and surveillance.

sar sensing has traditionally relied on airborne and spaceborne systems such
as satellites and aircraft. However, these types of platforms lack flexibility
and their revisit time for a given aoi is often large [2]. In this context,
the recent growth of the drone industry has motivated the use of sar onboard
small and lightweight uav, which is promising and challenging at the same
time. On the one hand, sar sensing benefits from the high flexibility and the
easy deployment of small drones at relatively short distances from aoi.
Consequently, uav-based sar sensing, which is now possible at low cost,
facilitates short revisit times and enables various new applications, such as
the monitoring of highly dynamic environments. On the other hand, the drones’
limited hardware capacity is confronted with highly computationally complex
sar image processing algorithms[7]. This makes real-time sar, which is an
emerging sar application [8, 9], difficult to realize onboard the flying
platform. In this regard, one potential solution is to transmit the sar raw
data to the ground for processing, which requires a reliable communication
link between the uav and the ground [10, 11, 12]. In real-time sar, the radar
raw data received for each sensing pulse is processed with minimum delay such
that the sar image frame can be generated as a function of the received sar
raw data. This facilitates immediate adjustments based on the processed sar
data. Moreover, uav instability presents another challenge facing uav-sar
systems. In fact, lightweight drones are vulnerable to external perturbations
such as wind and atmospheric turbulence [13]. In this context, robust
trajectory and resource allocation design wrt uav jittering were investigated
for uav-based communication systems [14, 15]. Another challenge is the drones’
limited battery capacity [16]. Optimizing the uav’s energy consumption
requires a careful trajectory design as energy consumption and sar sensing
performance are tightly coupled.

In general, trajectory design and resource allocation for uav have been
exhaustively studied in the literature [17, 18] and different optimization
techniques ranging from classical optimization [19] to learning-based methods
[20, 21] have been used to tackle the problem. However, for the uav-sar use
case, trajectory design and resource allocation have not been comprehensively
addressed. In fact, existing experimental works on uav-sar systems only
investigate the practical deployment of the sar sensors onboard small drones
and report results of experimental measurement campaigns with intuitive pre-
planned flight paths with fixed velocity and altitude [22, 23, 24].
Furthermore, existing experimental research works do not account for real-time
sar applications and are limited to offline image processing without
considering a potential wireless communication link to the ground. Among the
few works targeting the optimization of uav-sar systems, the authors in [25]
jointly optimized a focus measure for the sar images and the estimation of the
platform trajectory. However, while the estimation error of the actual
trajectory was minimized based on a nominal trajectory, the trajectory itself
was not optimized in [25], and communication with the ground was not
considered. Furthermore, for a bistatic uav-sar system, the authors in [26]
investigated trajectory optimization for a cellular-connected bistatic uav-sar
system, where the aoi was illuminated by a ground base station, and the energy
consumption was minimized. Yet, the proposed solution is limited to bistatic
systems with a stationary ground transmitter and a moving receiver, making it
inapplicable for active monostatic systems, where both transmitter and
receiver are mounted on the same moving platform. The authors in [27, 28]
investigated the use of a geosynchronous sar system, where the aoi was
illuminated by a satellite and echoes were collected by a passive sensing uav.
For instance, in [27], a multi-objective optimization problem for minimization
of the distance traveled by the uav-sar was formulated and solved. However,
the trajectory design of passive bistatic uav-sar systems is fundamentally
different from that of active monostatic uav-sar systems.

In this paper, we investigate the joint 3d trajectory and resource allocation
design for communication-assisted uav-based sar sensing. Different from other
works, we focus on active monostatic sar sensing, where the uav-sar enables
real-time sar applications by transmitting the collected radar data
instantaneously to the ground for processing. To this end, the uav 3d
trajectory and resources are jointly optimized to maximize the radar ground
coverage under communication, radar, and energy constraints. Our main
contributions are summarized as follows:

  * •

The joint 3d trajectory and resource allocation algorithm design is formulated
as a non-convex minlp. Our problem formulation takes into account the
constraints imposed by both the radar and communication subsystems.

  * •

The unavoidable uav trajectory deviations are modeled based on their
statistics and their effect on the uav-sar ground radar coverage is unveiled.

  * •

Exploiting successive convex approximation (SCA), we propose a low-complexity
sub-optimal solution for the formulated problem, which is robust wrt uav
trajectory deviations. Furthermore, based on monotonic optimization theory, we
provide an upper bound on the optimal solution to assess the performance of
the proposed sub-optimal algorithm.

  * •

Our simulation results show that the proposed sub-optimal algorithm achieves
close-to-optimal performance. In addition, we demonstrate the superiority of
the proposed solution in terms of radar coverage compared to several benchmark
schemes and confirm its robustness against random trajectory deviations.

We note that this article extends the corresponding conference version [1]. In
[1], only a non-robust sub-optimal algorithm was provided, and random uav
trajectory deviations were not accounted for. Moreover, the gap between the
sub-optimal scheme and the optimal solution was not assessed in [1].

The remainder of this article is organized as follows. In Section II, we
present the system model for the considered uav-sar system. In Section III, we
model random uav trajectory deviations and propose a concept for robust
trajectory design. In Section IV, we first formulate an optimization problem
for the joint 3d trajectory and resource allocation design for uav-sar
systems. Then, a robust sub-optimal solution and a corresponding upper bound
are provided. In Section V, we evaluate the performance of the proposed scheme
using the derived upper bound and other benchmark schemes. Finally, Section VI
concludes the paper and provides some directions for future research.

Notations: In this paper, lower-case letters $x$ refer to scalar numbers,
boldface lower-case letters $\mathbf{x}$ denote vectors, and boldface upper-
case letters $\mathbf{X}$ refer to matrices. $\\{a,...,b\\}$ denotes the set
of all integers between $a$ and $b$. $|\cdot|$ denotes the absolute value
operator, whereas $||\cdot||_{2}$ refers to the $L_{2}$-norm. Operators
$[x]^{+}$ and $\lceil x\rceil$ correspond to $\max(0,x)$ and the ceiling
function of $x$, respectively. For given sets $\mathcal{X}_{1}$ and
$\mathcal{X}_{2}$, $\mathcal{X}_{1}\setminus\mathcal{X}_{2}$ is the set that
comprises all elements of $\mathcal{X}_{1}$ that are not in $\mathcal{X}_{2}$.
$\mathbb{N}$ represents the set of natural numbers. For arbitrary integer $k$,
$\mathbb{N}_{k}$ is the set $\\{1,..,k\\}$. $\mathbb{R}^{N\times 1}$
represents the set of $N\times 1$ vectors with real-valued entries,
$\mathbb{R}^{N\times 1}_{+}$ denotes the set of $N\times 1$ vectors with non-
negative real entries, and $\mathbb{R}^{N\times M}$ denotes the set of real
$N\times M$ matrices. For $a,b\in\mathbb{N}$, $a$ mod $b$ is the remainder of
the division of $a$ by $b$. For vector $\mathbf{x}$, vector $\mathbf{x}^{T}$
stands for its transpose. For vectors
$\mathbf{x}=[x(1),...,x(N)]^{T}\in\mathbb{R}^{N\times 1}$ and
$\mathbf{y}=[y(1),...,y(N)]^{T}\in\mathbb{R}^{N\times 1}$,
$\mathbf{x}\leq\mathbf{y}$ and $\mathbf{x}<\mathbf{y}$ define element-wise
order relations on vectors $\mathbf{x}$ and $\mathbf{y}$, i.e., $x(n)\leq
y(n)$ and $x(n)<y(n)$, $\forall n$, respectively. For
$\mathbf{A}\in\mathbb{R}^{N\times N}$, $\mathbf{A}\succeq\mathbf{0}$ indicates
that $\mathbf{A}$ is a positive semi-definite matrix. For a given random
variable $X\in\mathbb{R}$, $F_{X}$ denotes its cumulative distribution
function (CDF) and $X\sim\mathcal{N}(\mu,\sigma)$ means that $X$ is Gaussian
distributed with mean $\mu$ and standard deviation $\sigma$.
$\mathbb{P}(\cdot)$ denotes the probability operator. For a real-valued
function $f(x)$, $f^{\prime}(x)$ denotes the derivative of $f$. $\rm
erf(\cdot)$ is the Gaussian error function. In Table I, we introduce the
definition of most important variables used in this paper.

TABLE I: Definition of most important variables.

Variable Definition  Variable Definition  $L$ Length of the aoi
$\mathrm{SNR}_{n}(\mathbf{P}_{\mathrm{sar}},\mathbf{u})$ Instantaneous radar
snr  $N$ Number of azimuth sweeps  $R_{\rm sl}$ Synchronization and
localization data size  $M$ Number of time slots per azimuth scan
$R_{\mathrm{min},n}(\mathbf{u})$ Instantaneous sar data rate  $z_{\rm max}$
and $z_{\rm min}$ uav maximum and minimum altitudes
$R_{n}(\mathbf{P}_{\mathrm{com}},\mathbf{u})$ Instantaneous uav data rate
$\mathbf{u}(n)$ uav 3d location in time slot $n$ $B_{c}$ Communication
bandwidth  $q(n)$ Remaining battery capacity in time slot $n$ $\gamma$
Reference channel power gain  $P_{\rm sar}(n)$ Instantaneous sar transmit
power  $\mathbf{g}$ gs 3d location  $P_{\rm com}(n)$ Instantaneous
communication transmit power  $d_{n}(\mathbf{u})$ uav-gs distance in time slot
$n$ $P_{\rm sar}^{\rm max}$ Maximum sar transmit power  $v$ uav velocity
$P_{\rm com}^{\rm max}$ Maximum communication transmit power  $q_{\rm start}$
Initial uav battery capacity  $\mathbf{P}_{\rm sar}$ sar transmit power vector
$P_{\rm prop}$ uav propulsion power  $\mathbf{P}_{\rm com}$ Communication
transmit power vector
$E_{n}(\mathbf{P}_{\mathrm{sar}},\mathbf{P}_{\mathrm{com}})$ Total energy
consumed during time slot $n$ $\tau_{p}$ and $\mathrm{PRF}$ Pulse duration and
repetition frequency  $S_{n}(\mathbf{u})$ uav sar swath width in time slot $n$
$\theta_{d}$ and $\theta_{\rm 3db}$ Radar depression angle and beamwidth
$C(\mathbf{u})$ Total sar ground coverage  $B_{r}$ Radar bandwidth  $r$
Coverage reliability level  $f$ Radar center frequency
$\bm{\delta}_{\mathbf{p}_{N}}^{r}$ Near-range compensation vector
$\mathrm{SNR_{min}}$ Minimum radar snr  $\bm{\delta}_{\mathbf{p}_{F}}^{r}$
Far-range compensation vector

##  II System Model

Figure 1: System model for uav-sar system with real-time data offloading to a
GS.

We consider a rotary-wing uav equipped with a side-looking radar system
designed to perform sar. The main mission of the uav is to scan a certain aoi
on the ground. To enable live mission control of the uav-sar system, real-time
sar processing is facilitated by transmitting the sensing data to a gs in
real-time.

###  II-A Trajectory Design

As a stripmap sar mode [2] is employed, we assume that the aoi is rectangular
with a given length $L$. Without loss of generality, we place the origin of
the adopted coordinate system at the corner of the aoi, where the uav mission
starts, see Figure 1. To enable real-time processing of the sar images, we
assume that frequency-domain algorithms are employed due to their
computational efficiency[4, 8, 9]. These algorithms impose the use of a linear
uav trajectory with fixed velocity $v>0$ and zero acceleration[29]. The
$x$-axis and $y$-axis define the range and azimuth direction, respectively,
see Figure 1. We discretize time such that the length of the area scanned in
$y$-direction in time slot $n$ of duration $\delta_{t}$, denoted by
$\Delta_{s}$, is given by:

$\Delta_{s}=\delta_{t}v.$ (1)

The uav trajectory consists of $N$ linear azimuth scans, where each scan
corresponds to a traversal in the $y$-direction from one edge of the aoi
(i.e., $y=0$) to the other (i.e., $y=L$) and vice versa, see Figure 1. Here,
the use of a linear trajectory is imposed by the sar stripmap mode [2]. In
each azimuth scan, $M$ uniform pulses are transmitted from the radar
system111Here, $\Delta_{s}$ is set such that $M\Delta_{s}=L$., making the
total number of time slots equal to $NM$. The drone location in time slot
$n\in\mathbb{N}_{NM}$ is defined by $\mathbf{u}(n)=[x(n),y(n),z(n)]^{T}$,
where $x(n)$ and $y(n)$ represent the position in the $xy$-plane, while $z(n)$
denotes the drone altitude. The uav trajectory
$\mathbf{u}=[\mathbf{u}(1),...,\mathbf{u}(NM)]$ is designed such that
successive azimuth scans are ideally adjacent with no gaps in the coverage. At
the end of each azimuth scan, the uav makes its turns and moves to the
appropriate $x$-position for the next scan as defined later in (4) and (5).
The uav altitude can be adjusted from scan to scan, but must be fixed during a
given azimuth scan to have a fixed width of the on-ground footprint of the
sensing beam along the range direction, i.e., fixed swath width. Consequently,
the radar power is also adjusted from one azimuth scan to the next, such that
the minimum required sensing snr is achieved, but is fixed during a given scan
. To account for these limitations of practical radar systems, we define two
sets of time slots, denoted by $\mathcal{A}$ and $\mathcal{A}^{c}$, as
follows:

$\displaystyle\mathcal{A}=\\{k\in\mathbb{N}_{NM}|(k-1)\;\mathrm{mod}\;M\neq
0\\},$ (2)
$\displaystyle\mathcal{A}^{c}=\mathbb{N}_{NM}\setminus\mathcal{A}=\\{1+(k-1)M,k\in\mathbb{N}_{N}\\}.$
(3)

Sets $\mathcal{A}$ and $\mathcal{A}^{c}$ will be used to (i) force a back-and-
forth linear motion and (ii) set all of the relevant radar parameters, such as
the radar power and altitude, constant during a given scan [30, 23]. To this
end, set $\mathcal{A}^{c}$ contains the $N$ indices of the first time slots at
the beginning of each azimuth scan, where the radar parameters can be
adjusted, whereas set $\mathcal{A}$ contains the indices of all other time
slots, where the parameters are fixed. The following constraints on the drone
trajectory are imposed:

$\displaystyle\mathrm{C1}:x(1)=-c_{1}\;z(1),$ (4)
$\displaystyle\mathrm{C2}:x(n)=x(n-1)+c_{2}\;z(n-1)-c_{1}\;z(n),\forall
n\in\mathcal{A}^{c}\setminus\\{1\\},$ (5)
$\displaystyle\mathrm{C3}:x(n+1)=x(n),\forall n\in\mathcal{A},$ (6)
$\displaystyle\mathrm{C4}:z(n+1)=z(n),\forall n\in\mathcal{A},$ (7)

Figure 2: Illustration of a Boustrophedon-shaped trajectory as imposed by
constraints $\mathrm{C1-C4}$ for $N=3$, $M=6$, and a set of altitudes
$\\{1,2,2\\}$ m.

where $c_{1}=\tan(\theta_{1})$, $c_{2}=\tan(\theta_{2})$,
$\theta_{1}=\theta_{d}-\frac{\Theta_{3\mathrm{dB}}}{2}$,
$\theta_{2}=\theta_{d}+\frac{\Theta_{3\mathrm{dB}}}{2}$,
$\Theta_{3\mathrm{dB}}$ is the radar antenna 3dB beam width, and $\theta_{d}$
is the radar depression angle, see Figure 1. Constraint $\mathrm{C2}$ ensures
that two successive radar azimuth scans of the ground are adjacent without any
gap between them. During each scan, constraint $\mathrm{C3}$ imposes a linear
motion, while constraint $\mathrm{C4}$ ensures a fixed altitude during a scan.
In a nutshell, constraints $\mathrm{C1-C4}$ give rise to a Boustrophedon-
shaped trajectory, as illustrated in Figure 2. Forcing a linear motion for
stripmap uav-sar does not mean that the trajectory is completely fixed. In
fact, in each turn, the drone’s $z$-position (altitude) and $x$-position can
be adjusted. Only the $y$-position is fixed such that back-and-forth stripmap
mode motion is guaranteed. The $y$-position is defined as follows:

$y(1)=0,\;y(n)=y(n-1)+c(n-1)\Delta_{s},\forall
n\in\mathbb{N}_{NM}\setminus\\{1\\},$ (8)

where vector $\mathbf{c}=[c(1),...,c(NM)]^{T}\in\mathbb{R}^{NM\times 1}$ is
defined as follows:

$\mathbf{c}=[\overbrace{\underbrace{1,...,1}_{\text{$M$
times}},\underbrace{-1,...,-1,}_{\text{$M$ times}}...]^{T}}^{\text{$N\times M$
times}}.$ (9)

###  II-B Radar Imaging

Since we operate at relatively low altitudes, the maximal imageable area on
the ground is limited by the antenna beam width. In time slot $n$, the width
in $x$-direction of the on-ground footprint of the sensing beam, also referred
to as sar swath width, is denoted by $S_{n}$ and given as follows:

$S_{n}(\mathbf{u})=\left(c_{2}-c_{1}\right)z(n).$ (10)

Following the trajectory $\mathbf{u}$, the total ground area covered by the
drone is the sum of all infinitesimal small swath widths, and therefore can be
approximated222The approximation is due to the elliptical shape of the beam
footprint on the ground and becomes negligible for large $M$.  as follows:

$C(\mathbf{u})=\sum_{n=1}^{NM}\Delta_{s}S_{n}(\mathbf{u}).$ (11)

The data rate produced by the sar sensing data in time slot $n$ can be
expressed as [31]:

$R_{\mathrm{min},n}(\mathbf{u})=\;B_{r}\left(\frac{2\;z(n)\;\Omega}{c}+\tau_{p}\right)\mathrm{PRF},$
(12)

where $B_{r}$ is the radar bandwidth, $\tau_{p}$ is the radar pulse duration,
$c$ is the speed of light, $\rm PRF$ is the radar pulse repetition frequency,
and $\Omega$ is given by:

$\Omega=\frac{\cos(\theta_{1})-\cos(\theta_{2})}{\cos(\theta_{1})\cos(\theta_{2})}.$
(13)

Furthermore, the achieved radar snr in time slot $n$ is given as follows [32]:

$\mathrm{SNR}_{n}(\mathbf{P}_{\mathrm{sar}},\mathbf{u})=\frac{P_{\mathrm{sar}}(n)\;G_{t}\;G_{r}\;\lambda^{3}\;\sigma_{0}\;c\;\tau_{p}\mathrm{PRF}\;\sin^{2}(\theta_{d})}{(4\pi)^{4}\;z^{3}(n)\;{k_{B}}\;T_{o}\;NF\;\;B_{r}\;L_{\mathrm{tot}}\;v},$
(14)

where
$\mathbf{P}_{\mathrm{sar}}=[P_{\mathrm{sar}}(1),...,P_{\mathrm{sar}}(NM)]^{T}$,
$P_{\mathrm{sar}}(n)$ is the radar transmit power in time slot $n$, $G_{t}$
and $G_{r}$ are the radar antenna gains for transmission and reception,
respectively, $\lambda$ is the radar wavelength, $\sigma_{o}$ is the
backscattering coefficient, ${k_{B}}$ is Boltzmann’s constant, $T_{o}$ is the
equivalent noise temperature, $NF$ is the system noise figure, and
$L_{\mathrm{tot}}$ represents the total radar losses.

###  II-C uav Data Link

We denote the location of the gs by $\mathbf{g}=[g_{x},g_{y},g_{z}]^{T}$. In
time slot $n$, the distance between the drone and the gs is given by:

$\displaystyle
d_{n}(\mathbf{u})=\sqrt{(x(n)-g_{x})^{2}+(y(n)-g_{y})^{2}+(z(n)-g_{z})^{2}}.$
(15)

We assume a los communication link to the gs. Thus, based on the free space
path loss model333 The offline optimization framework proposed in this paper
can be extended to other communication channels with small-scale and large-
scale fading components by an adequate design of the link margin that accounts
for the channel power fluctuations [33]., the instantaneous throughput of the
backhaul link between the uav and the gs is given by:

$R_{n}(\mathbf{P}_{\mathrm{com}},\mathbf{u})=B_{c}\;\log_{2}\left(1+\frac{P_{\mathrm{com}}(n)\;\gamma}{d^{2}_{n}(\mathbf{u})}\right),$
(16)

where
$\mathbf{P}_{\mathrm{com}}=[P_{\mathrm{com}}(1),...,P_{\mathrm{com}}(NM)]^{T}$,
$P_{\mathrm{com}}(n)$ is the power allocated for communication in time slot
$n$, $B_{c}$ is the communication bandwidth, and $\gamma$ is the channel power
gain at a reference distance of 1 m divided by the noise variance. To
guarantee successful real-time backhauling of the sensing data collected by
the drone to the gs, the following inequality must be satisfied at all times:

$\displaystyle R_{n}(\mathbf{P}_{\mathrm{com}},\mathbf{u})\geq
R_{\mathrm{min},n}(\mathbf{u})+R_{\mathrm{sl}},\forall n\in\mathbb{N}_{NM},$
(17)

where $R_{\mathrm{sl}}$ is the fixed data rate needed to transmit the
synchronization and localization data necessary for successful radar image
processing. In this work, we assume that the communication system operates in
a licensed frequency band that is different from the one used for sensing,
thus, avoiding interference that could lead to transmission failures. Lastly,
the total energy consumed by the uav in time slot $n$ is given by:

$E_{n}(\mathbf{P}_{\mathrm{com}},\mathbf{P}_{\mathrm{sar}})=\delta_{t}\left(\mathbf{P}_{\mathrm{com}}(n)+\mathbf{P}_{\mathrm{sar}}(n)+P_{\mathrm{prop}}\right),$
(18)

where $P_{\mathrm{prop}}$ is the propulsion power required for maintaining the
operation of the drone and is given by [34]:

$P_{\mathrm{prop}}=P_{0}\left(1+\frac{3v^{2}}{U^{2}_{\mathrm{tip}}}\right)+P_{I}\left(\sqrt{1+\frac{v^{4}}{4v_{0}^{4}}}-\frac{v^{2}}{2v_{0}^{2}}\right)^{\frac{1}{2}}+\frac{1}{2}d_{0}\rho
sAv^{3}.$ (19)

Here, $P_{0}$ and $P_{I}$ are two constants defined in [34],
$U_{\mathrm{tip}}$ is the tip speed of the rotor blade, $d_{0}$ is the
fuselage drag ratio, $\rho$ is the air density, $s$ is the rotor solidity, $A$
is the rotor disc area, $v_{0}=\sqrt{\frac{W_{u}}{2\rho A}}$ is the mean rotor
induced velocity in hover, and $W_{u}$ is the uav weight in Newton. Note that,
in this work, the propulsion power is constant since the drone’s velocity,
$v$, is fixed to allow for efficient signal processing of the sar images.

##  III Robust uav Trajectory Design Under Random Position Deviations

In practice, it is difficult for lightweight drones to strictly follow a
designed trajectory due to perturbations caused by weather conditions,
mechanical vibrations, and localization errors [13]. These errors may result
in areas not covered by sar. In this section, we first develop a model for uav
trajectory deviations and quantify their impact on the footprint of the
sensing beam. Then, we propose a statistics-based robust trajectory design
that suitably widens and adjusts the sensing beam footprint. This creates
overlapping scans and allows covering the aoi without gaps.

###  III-A Trajectory Deviation Model

|  |
---|---|---
(a) | (b) | (c)

|  |
---|---|---
(d) | (e) | (f)

Figure 3:  Robust uav trajectory design. (a) Successive azimuth scans for an
ideal trajectory. (b) Effect of in-range error on the ground sar swath. (c)
Effect of altitude error on the ground sar swath. (d) Combined effect of in-
range and altitude errors on the ground sar swath. (e) Uncovered areas due to
uav flight deviations. (f) Proposed robust uav-sar trajectory.

If the actual trajectory is different from the ideal linear trajectory
$\mathbf{u}$, described by $\rm C1-C4$, the actual footprint of the sensing
beam is also affected in terms of both its on-ground position and width. Let
position vectors
$(\mathbf{x}^{a},\mathbf{y}^{a},\mathbf{z}^{a})\in\mathbb{R}^{NM\times 1}$
describe the actual trajectory, denoted by
$\mathbf{u}^{a}=[\mathbf{u}^{a}(1),...,\mathbf{u}^{a}(NM)]$, where the drone’s
actual position in time slot $n$ is given by
$\mathbf{u}^{a}(n)=[x^{a}(n),y^{a}(n),z^{a}(n)]^{T}$. In what follows, we
quantify the change in the actual footprint of the sensing beam due to random
uav trajectory deviations.

####  III-A1 Cross-range error

The deviation in $y$-direction, denoted by
$\bm{\Delta}_{y}=\mathbf{y}^{a}-\mathbf{y}=$
$[\Delta_{y}(1),...,\Delta_{y}(NM)]^{T}\in\mathbb{R}^{NM\times 1}$, has the
least effect on the sar ground coverage. Although this error can affect the
sar resolution by influencing the illumination time of a given reference point
on the ground, it has no effect on the radar coverage thanks to the radar beam
width in $y$-direction. In fact, the main idea behind sar systems is to use
this direction to create the synthetic aperture such that any reference point
on the ground is illuminated multiple times from different angles along the
$y$-axis [2]. Therefore, in what follows, we neglect the effect of this error
on the radar coverage and focus on uav trajectory deviations in the
$xz$-plane.

####  III-A2 In-range error

This error is caused by the uav position deviation in $x$-direction and is
denoted by vector $\bm{\Delta}_{x}=\mathbf{x}^{a}-\mathbf{x}=$
$[\Delta_{x}(1),...,\Delta_{x}(NM)]^{T}$ $\in$ $\mathbb{R}^{NM\times 1}$, see
Figure 3b. As shown in the figure, this error causes near-range and far-range
shifts of the ground footprint of the sensing beam along the $x$-axis. To this
end, in the $xz-$plane and in time slot $n$, we define point
$\mathbf{p}_{N}(n)=[x(n)+c_{1}z(n),y(n),0]^{T}$ representing the intersection
between the near-range slant distance line and the $x$-axis and point
$\mathbf{p}_{F}(n)=[x(n)+c_{2}z(n),y(n),0]^{T}$ representing the intersection
between the far-range slant distance line and the $x$-axis. Both points are
defined wrt the drone’s ideal trajectory, see Figure 3b. Any deviation from
the ideal uav position in the $xz$-plane results in a near-range error (shift
of point $\mathbf{p}_{N}(n)$ along the $x$-axis), denoted by
$\bm{\Delta}_{\mathbf{p}_{N}}=[\Delta_{\mathbf{p}_{N}}(1),...,\Delta_{\mathbf{p}_{N}}(NM)]^{T}\in\mathbb{R}^{NM\times
1}$, and a far-range error (shift of point $\mathbf{p}_{F}(n)$ along the
$x$-axis), denoted by
$\bm{\Delta}_{\mathbf{p}_{F}}=[\Delta_{\mathbf{p}_{F}}(1),...,\Delta_{\mathbf{p}_{F}}(NM)]^{T}\in\mathbb{R}^{NM\times
1}$, of the sensing beam footprint. A deviation in $x$-direction, causes
identical near-range and far-range errors such that, in time slot $n$,
$\Delta_{\mathbf{p}_{F}}(n)=\Delta_{\mathbf{p}_{N}}(n)=\Delta_{x}(n)$.
Therefore, an in-range error does not affect the width of the sar swath, but
results in a translation of the footprint of the sensing beam along the
$x$-axis as shown in Figure 3b.

####  III-A3 Altitude error

A deviation of the trajectory in $z$-direction, given by
$\bm{\Delta}_{z}=\mathbf{z}^{a}-\mathbf{z}=[\Delta_{z}(1),...,\Delta_{z}(NM)]^{T}\in\mathbb{R}^{NM\times
1}$, results in non-equal far-range and near-range errors. In fact, an
altitude error in time slot $n$, denoted by $\Delta_{z}(n)$, results in
$\Delta_{\mathbf{p}_{N}}(n)=c_{1}\Delta_{z}(n)$ and
$\Delta_{\mathbf{p}_{F}}(n)=c_{2}\Delta_{z}(n)$. In contrast to in-range
errors, altitude errors induce a change in both the width and the position of
the sensing beam footprint, see Figure 3c.

As shown in Figure 3d, in time slot $n$, the combined in-range and altitude
errors result in a near-range error of the footprint of the sensing beam given
by:

$\Delta_{\mathbf{p}_{N}}(n)=\Delta_{x}(n)+c_{1}\Delta_{z}(n).$ (20)

Similarly, in time slot $n$, the total far-range error is given by:

$\Delta_{\mathbf{p}_{F}}(n)=\Delta_{x}(n)+c_{2}\Delta_{z}(n).$ (21)

####  III-A4 Trajectory deviation model

In this work, based on measurements in [35, 36], we propose to model the uav-
sar trajectory deviations as independent Gaussian random variables such that,
in each time slot $n$, $\Delta_{x}(n)\sim\mathcal{N}(o_{x},\sigma)$ and
$\Delta_{z}(n)\sim\mathcal{N}(o_{z},\sigma)$, $\forall n$. The standard
deviation, denoted by $\sigma$, depends on the stability of the uav in the
presence of external perturbations caused, e.g., by adverse weather
conditions. The offsets in $x$\- and $z$-direction, denoted by $o_{x}$ and
$o_{z}$, respectively, represent a possible constant error between the uav’s
actual trajectory and its ideal intended trajectory and depend on the uav
localization accuracy. For the special case, $o_{x}=o_{z}=0$, the uav
trajectory deviations reduce to random fluctuations around the ideal
trajectory, also referred to as uav jittering [14].

###  III-B Proposed Robust Trajectory Design

Both positive near-range errors, i.e., $\Delta_{\mathbf{p}_{N}}(n)>0$, and
negative far-range errors, i.e., $\Delta_{\mathbf{p}_{F}}(n)<0$, result in
uncovered areas, whereas the reverse result in an overlap between successive
scans, see Figure 3e. For robust trajectory design, uncovered areas are
avoided by tolerating an overlap between adjacent azimuth scans. This can be
achieved by shifting the actual footprint of the sensing beam such that the
near-range and far-range errors are compensated. Since excessive overlap
between scans reduces the sar coverage, the overlap should be limited to the
degree needed for compensation of uav trajectory deviations. To this end, we
define the coverage reliability level $r$, $0\leq r\leq 1$, as the probability
of having no gaps between successive sar azimuth scans. The compensation that
is required in the near range of the sensing beam footprint to achieve
coverage reliability level $r$ is denoted by
$\bm{\delta}_{\mathbf{p}_{N}}^{r}=[\delta_{\mathbf{p}_{\scriptscriptstyle
N}}^{r}(1),...,\delta_{\mathbf{p}_{N}}^{r}(NM)]^{T}$ $\in\mathbb{R}^{NM\times
1}$, and is obtained as:

$\delta_{\mathbf{p}_{N}}^{r}(n)\stackrel{{\scriptstyle}}{{=}}\operatorname*{arg\,min}_{s}\Big{\\{}|s|\Big{|}s\in\mathbb{R},\mathbb{P}\left(\Delta_{\mathbf{p}_{N}}(n)+s\leq
0\right)\geq r\Big{\\}}.$ (22)

Similarly, the far-range compensation, which we denote by
$\bm{\delta}_{\mathbf{p}_{F}}^{r}=$
$[\delta_{\mathbf{p}_{F}}^{r}(1),...,\delta_{\mathbf{p}_{F}}^{r}(NM)]^{T}$
$\in\mathbb{R}^{NM\times 1}$, is obtained as:

$\delta_{\mathbf{p}_{F}}^{r}(n)\stackrel{{\scriptstyle}}{{=}}\operatorname*{arg\,min}_{s}\Big{\\{}|s|\Big{|}s\in\mathbb{R},\mathbb{P}\left(\Delta_{\mathbf{p}_{F}}(n)+s\geq
0\right)\geq r\Big{\\}}.$ (23)

Note that, on the one hand, to avoid coverage holes, coverage reliability
values close to 1 are desirable. On the other hand, if $r$ is chosen too close
to 1, this might result in a reduction of the total coverage, as a larger
swath overlap is required to compensate for random trajectory deviations.
Based on our simulations, $r=0.95$ provides a good compromise between coverage
reliability and total coverage.
Next, given the trajectory deviation statistics, we derive analytical
expressions for the near-range and far-range compensations.

######  Proposition 1.

The near-range and far-range compensations given in (22) and (23),
respectively, can be obtained in time slot $n$ as follows:

$\displaystyle\delta_{\mathbf{p}_{N}}^{r}(n)$
$\displaystyle=-\Big{[}\mathrm{erf}^{-1}(2r-1)\sigma\sqrt{2(1+{c_{1}}^{2})}+o_{x}+c_{1}o_{z}\Big{]}^{+},$
(24) $\displaystyle\delta_{\mathbf{p}_{F}}^{r}(n)$
$\displaystyle=\Big{[}\mathrm{erf}^{-1}(2r-1)\sigma\sqrt{2(1+{c_{2}}^{2})}-o_{x}-c_{2}o_{z}\Big{]}^{+}.$
(25)

###### Proof.

Please refer to Appendix A. ∎

######  Remark 1.

Note that $\delta_{\mathbf{p}_{N}}^{r}(n)\leq 0$, whereas
$\delta_{\mathbf{p}_{F}}^{r}(n)\geq 0$, $\forall n\in\mathbb{N}_{NM}$, which
means that the actual footprint of the radar beam is widened, see Figure 3f.
This creates overlapping scans that prevent gaps between adjacent scans.

The proposed near-range and far-range compensations can be used to adjust the
uav position. In the following proposition, we translate the near-range and
far-range compensation of the beam footprint to an in-range ($x$-position)
adjustment, denoted by
$\bm{\delta}_{x}^{r}=[\delta_{x}^{r}(1),...,\delta_{x}^{r}(NM)]^{T}\in\mathbb{R}^{NM\times
1}$, and an altitude ($z$-position) adjustment, denoted by
$\bm{\delta}_{z}^{r}=[\delta_{z}^{r}(1),...,\delta_{z}^{r}(NM)]^{T}\in\mathbb{R}^{NM\times
1}$, as shown in Figure 3f.

######  Proposition 2.

Given the near-range and far-range compensation of the beam footprint,
$\bm{\delta}_{\mathbf{p}_{N}}^{r}$ and $\bm{\delta}_{\mathbf{p}_{F}}^{r}$, the
corresponding in-range and altitude adjustment of the uav position are
respectively given by:

$\displaystyle\bm{\delta}_{x}^{r}=\bm{\delta}_{\mathbf{p}_{N}}^{r}-c_{1}\;\bm{\delta}_{z}^{r},$
(26)
$\displaystyle\bm{\delta}_{z}^{r}=\frac{\bm{\delta}_{\mathbf{p}_{F}}^{r}-\bm{\delta}_{\mathbf{p}_{N}}^{r}}{c_{2}-c_{1}}.$
(27)

###### Proof.

Based on the geometry of the problem, In time slot $n$, given
$\delta_{\mathbf{p}_{N}}^{r}(n)\leq 0$ and $\delta_{\mathbf{p}_{F}}^{r}(n)\geq
0$, we shift the uav position along the $x$-axis by $\delta_{x}^{r}(n)\leq 0$,
and along the $z$-axis by $\delta_{z}^{r}(n)\geq 0$, see Figure 3f. The
required uav position adjustments can be derived based on the geometry of the
problem by solving the following system of equations:

$\displaystyle\begin{cases}c_{1}=\frac{c_{1}z(n)+\delta_{\mathbf{p}_{N}}^{r}(n)-\delta_{x}^{r}(n)}{z(n)+\delta_{z}^{r}(n)},\\\
c_{2}=\frac{c_{2}z(n)+\delta_{\mathbf{p}_{F}}^{r}(n)-\delta_{x}^{r}(n)}{z(n)+\delta_{z}^{r}(n)},\end{cases}$
(28)

where the first and second lines of (28) are obtained from triangles ($ABC$)
and ($ABD$) in Figure 3f, respectively. ∎

For robust uav-sar trajectory design, the uav position is adjusted in each
time slot to achieve a coverage reliability level $r$ in the presence of uav
trajectory deviations. Consequently, the robust trajectory is denoted by
$\textbf{u}^{r}=[\textbf{u}^{r}(1),...,\textbf{u}^{r}(NM)]$, where
$\mathbf{u}^{r}(n)$ is the adjusted position in time slot $n$ and is given by:

$\displaystyle\mathbf{u}^{r}(n)$
$\displaystyle=\left[x^{r}(n),y(n),z^{r}(n)\right]^{T}=\left[x(n)+\delta_{x}^{r}(n),y(n),z(n)+\delta_{z}^{r}(n)\right]^{T}.$
(29)

##  IV Robust Resource Allocation and Trajectory Optimization Framework

In this section, we formulate an optimization problem for joint 3d trajectory
and resource allocation design. The formulated problem is solved offline where
random uav trajectory deviations are accounted for via the robust trajectory
given in (29). In particular, for a given coverage reliability level $r$, we
provide a low-complexity sub-optimal solution of the formulated problem based
on SCA. Moreover, to assess the performance of the proposed sub-optimal
scheme, we derive an upper bound for the optimal solution based on monotonic
optimization theory [41].

###  IV-A Problem Formulation

Our objective is the maximization of the drone ground coverage $C$ while
satisfying the constraints imposed by the radar and communication subsystems.
We assume Gaussian trajectory deviation as described in Section III-A. To
achieve a coverage reliability level $r$, the uav follows the robust
trajectory $\mathbf{u}^{r}$ provided in (29). We optimize the drone trajectory
$\\{\mathbf{z},\mathbf{x},N\\}$ and the power allocation
$\\{\mathbf{q},\mathbf{P_{\mathrm{com}}},\mathbf{P_{\mathrm{sar}}}\\}$, where
vectors $\mathbf{z}$, $\mathbf{x}$, $\mathbf{q}$, $\mathbf{P_{\mathrm{com}}}$,
and $\mathbf{P_{\mathrm{sar}}}$, represent the collections of all
$z(n),x(n),q(n),P_{\rm com}(n),$ and $P_{\rm sar}(n),\forall
n\in\mathbb{N}_{NM}$, respectively, and $q(n)$ is the energy remaining in the
drone’s battery in time slot $n$. To this end, the following optimization
problem is formulated:

$\displaystyle\\!(\mathrm{P.1}):\max_{\mathbf{z},\mathbf{x},\mathbf{q},\mathbf{P_{\mathrm{sar}}},\mathbf{P_{\mathrm{com}}},N}C(\mathbf{u}^{r})$
$\displaystyle\text{s.t.}\hskip 5.69054pt\mathrm{C1-C4,}$
$\displaystyle\mathrm{C5:}\;P_{\mathrm{sar}}(n+1)=P_{\mathrm{sar}}(n),\forall
n\in\mathcal{A},$ $\displaystyle\mathrm{C6}:\;z_{\mathrm{min}}\leq
z^{r}(n)\leq z_{\mathrm{max}},\forall n\in\mathbb{N}_{NM},$
$\displaystyle\mathrm{C7}:\mathrm{SNR}_{n}(\mathbf{P}_{\mathrm{sar}},\mathbf{u}^{r})\geq\mathrm{SNR_{min}},\forall
n\in\mathbb{N}_{NM},$
$\displaystyle\mathrm{C8}:R_{n}(\mathbf{P}_{\mathrm{com}},\mathbf{u}^{r})\geq
R_{\mathrm{min},n}(\mathbf{u}^{r})+R_{\mathrm{sl}},\forall
n\in\mathbb{N}_{NM},$ $\displaystyle\mathrm{C9}:0\leq P_{\mathrm{sar}}(n)\leq
P_{\mathrm{sar}}^{\mathrm{max}},0\leq P_{\mathrm{com}}(n)\leq
P_{\mathrm{com}}^{\mathrm{max}},\forall n\in\mathbb{N}_{NM},$
$\displaystyle\mathrm{C10}:q(1)=q_{\mathrm{start}},\;q(n)\geq 0,\forall
n\in\mathbb{N}_{NM},$
$\displaystyle\mathrm{C11}:q(n+1)=q(n)-E_{n}(\mathbf{P}_{\mathrm{sar}},\mathbf{P}_{\mathrm{com}}),\forall
n\in{\mathbb{N}_{NM-1}},$ $\displaystyle\mathrm{C12}:N\in\mathbb{N}.$

Constraints $\mathrm{C1-C4}$ define the shape of the uav ideal trajectory as
described in (4)-(7). Constraint $\mathrm{C5}$ ensures constant radar power
along each azimuth scan. Constraint $\mathrm{C6}$ determines the allowed range
for the uav’s operational altitude when the drone follows the proposed robust
trajectory. Constraint $\mathrm{C7}$ ensures the minimum snr,
$\mathrm{SNR_{min}}$, required for sar imaging when the drone follows the
robust trajectory $\mathbf{u}^{r}$ given by (29). Constraint $\mathrm{C8}$
guarantees successful real-time transmission of the sar data collected by the
uav to the gs. Constraint $\mathrm{C9}$ ensures that the radar and
communication transmit powers are non-negative and do not exceed the maximum
permissible levels. Constraint $\mathrm{C10}$ specifies the energy available
in the uav battery, $q_{\mathrm{start}}$, at the start of the mission, and
ensures that the energy remaining in the battery during operation, $q(n)$,
does not become negative. In addition, constraint $\mathrm{C11}$ updates the
uav battery level in time slot $n+1$ by subtracting the total energy consumed
during time slot $n$, denoted by
$E_{n}(\mathbf{P}_{\mathrm{sar}},\mathbf{P}_{\mathrm{com}})$, from $q(n)$.
Constraint $\mathrm{C12}$ defines the number of azimuth scans realized by the
uav which is an integer.

Problem $\rm(P.1)$ is a non-convex minlp due to non-convex constraint
$\mathrm{C8}$ and integer variable constraint $\mathrm{C12}$. The latter does
not only affect the objective function but also the dimension of all other
optimization variables. Therefore, $N$ cannot be optimized jointly with the
remaining variables. Furthermore, non-convex minlp are usually very difficult
to solve. Nevertheless, in the following, we provide a low-complexity sub-
optimal solution to problem $\rm(P.1)$ based on SCA.

###  IV-B Low-complexity Sub-optimal Solution of Problem $\mathrm{(P.1)}$

In this subsection, we provide a low-complexity sub-optimal solution to
problem $\rm(P.1)$ based on sca. First, we fix the number of azimuth scans $N$
and denote the resulting problem by $\mathrm{(P.2)}$:

$\displaystyle\\!\mathrm{(P.2)}:\max_{\mathbf{z},\mathbf{x},\mathbf{q},\mathbf{P_{\mathrm{sar}}},\mathbf{P_{\mathrm{com}}}}C(\mathbf{u}^{r})$
$\displaystyle\text{s.t.}\hskip 8.53581pt\mathrm{C1-C11.}$

However, even for a given $N$, problem $\mathrm{(P.2)}$ is still non-convex
and difficult to solve. Thus, we approximate problem $\rm(P.2)$ as a convex
problem and use iterative sca to solve it. Then, we perform a search for the
optimal number of azimuth scans $N^{*}$. As a first step, we transform
constraint $\rm C7$, which involves a cubic function of drone altitude $z$,
using second-order cone programming. To this end, we introduce slack variables
$\psi(n)$ stacked in vector
$\bm{\psi}=[\psi(1),...,\psi(NM)]^{T}\in\mathbb{R}^{NM\times 1}$, and
reformulate $\rm C7$ as follows:

$\displaystyle\mathrm{C7a}$ $\displaystyle:\begin{bmatrix}\psi(n)&z^{r}(n)\\\
z^{r}(n)&1\end{bmatrix}\succeq\mathbf{0},\forall n\in\mathbb{N}_{NM},$ (30)
$\displaystyle\mathrm{C7b}$ $\displaystyle:\begin{bmatrix}P_{\rm
sar}(n)\beta&\psi(n)\\\ \psi(n)&z^{r}(n)\end{bmatrix}\succeq\mathbf{0},\forall
n\in\mathbb{N}_{NM},$ (31) $\displaystyle\mathrm{C7c}$
$\displaystyle:\psi(n)\geq 0,\forall n\in\mathbb{N}_{NM},$ (32)

where
$\beta=\frac{G_{t}\;G_{r}\;\lambda^{3}\;\sigma_{0}\;c\;\tau_{p}\mathrm{PRF}\;\sin^{2}(\theta_{d})}{(4\pi)^{4}\;{k_{B}}\;T_{o}\;NF\;\;B_{r}\;L_{\mathrm{tot}}\;v\;\mathrm{SNR_{min}}}$.
Furthermore, non-convex constraint $\rm C8$ can be rewritten as follows:

${\rm C8}:\left(A2^{\alpha z^{r}(n)}-1\right)d^{2}_{n}(\mathbf{u}^{r})\leq
P_{\mathrm{com}}(n)\gamma,\forall n\in\mathbb{N}_{NM},$ (33)

where $A=2^{\frac{B_{r}\;\tau_{p}\;\mathrm{PRF}}{B_{c}}}$ and
$\alpha=\frac{2\;\Omega\;B_{r}\;\mathrm{PRF}}{c\;B_{c}}$. In time slot $n$,
let function $H_{n}(\mathbf{u}^{r})$ be equal to the left-hand side of (33).
Function $H_{n}$ is non-convex wrt to optimization variables $x$ and $z$.
Therefore, in each sca iteration $j\in\mathbb{N}$, we use a Taylor series
approximation around points $\\{z^{(j)}(n),x^{(j)}(n)\\}\in\mathbb{R}$ to
convexify constraint $\rm C8$. This is achieved by replacing function $H_{n}$
with a global overestimate, denoted by $\overline{H}_{n}$. To this end, in
time slot $n$, we first expand function $H_{n}$ as follows:

$\displaystyle
H_{n}(\mathbf{u}^{r})=h_{1}(z(n))h_{2}(z(n))+h_{1}(z(n))h_{3}(x(n))+h_{1}(z(n))(y(n)-g_{y})^{2},$
(34)

where functions $h_{1},h_{2},$ and $h_{3}$ are convex differentiable functions
that are given by:

$\displaystyle h_{1}(z(n))$ $\displaystyle=A2^{\alpha z^{r}(n)}-1,$ (35)
$\displaystyle h_{2}(z(n))$ $\displaystyle=(z^{r}(n)-g_{z})^{2},$ (36)
$\displaystyle h_{3}(x(n))$ $\displaystyle=(x^{r}(n)-g_{x})^{2}.$ (37)

Yet, terms $h_{1}(z(n))h_{2}(z(n))$ and $h_{1}(z(n))h_{3}(x(n))$ contain
coupling and are non-convex w.r.t. optimization variables $x$ and $z$. To
overcome this obstacle, we rewrite them as follows:

$\displaystyle h_{1}(z(n))h_{2}(z(n))$
$\displaystyle=\frac{1}{2}\big{(}h_{1}(z(n))+h_{2}(z(n))\big{)}^{2}-\frac{1}{2}h_{1}^{2}(z(n))-\frac{1}{2}h_{2}^{2}(z(n)),$
(38) $\displaystyle h_{1}(z(n))h_{3}(x(n))$
$\displaystyle=\frac{1}{2}\big{(}h_{1}(z(n))+h_{3}(x(n))\big{)}^{2}-\frac{1}{2}h_{1}^{2}(z(n))-\frac{1}{2}h_{3}^{2}(x(n)).$
(39)

Equations (38) and (39) are differences of convex functions, where functions
$f_{i}=\frac{1}{2}h_{i}^{2},i\in\\{1,2,3\\}$, represent an obstacle to solving
problem $\rm(P.2)$. Therefore, we replace them with their global
underestimates around points $\\{z^{(j)}(n),x^{(j)}(n)\\}$ based on Taylor
series approximation:

$\displaystyle\underline{f_{i}}(z(n))$
$\displaystyle=f_{i}(z^{(j)}(n))+f_{i}^{\prime}(z^{(j)}(n))(z(n)-z^{(j)}(n)),\forall
i\in\\{1,2\\},$ (40) $\displaystyle\underline{f_{3}}(x(n))$
$\displaystyle=f_{3}(x^{(j)}(n))+f_{3}^{\prime}(x^{(j)}(n))(x(n)-x^{(j)}(n)),\forall
n.$ (41)

Consequently, the global overestimate for function $\overline{H}_{n}$ around
$\\{z^{(j)}(n),x^{(j)}(n)\\}$ is given by:

$\displaystyle\overline{H}_{n}(\mathbf{u}^{r})=$
$\displaystyle\frac{1}{2}\big{(}h_{1}(z(n))+h_{2}(z(n))\big{)}^{2}+\frac{1}{2}\big{(}h_{1}(z(n))+h_{3}(x(n))\big{)}^{2}-$
$\displaystyle
2\underline{f_{1}}(z(n))-\underline{f_{2}}(z(n))-\underline{f_{3}}(x(n))+h_{1}(z(n))(y(n)-g_{y})^{2},\forall
n.$ (42)

In the final step, we introduce two new slack vectors
$\mathbf{t}=[t(1),...,t(NM)]^{T}\in\mathbb{R}^{NM\times 1}$,
$\mathbf{o}=[o(1),...,o(NM)]^{T}\in\mathbb{R}^{NM\times 1}$, and two new
constraints $\rm\widetilde{C8a}$ and $\rm\widetilde{C8b}$ to reformulate
constraint $\rm C8$ as follows:

$\displaystyle\widetilde{\rm
C8}:\frac{1}{2}t^{2}(n)+\frac{1}{2}o^{2}(n)-2\underline{f_{1}}(z(n))-\underline{f_{2}}(z(n))-\underline{f_{3}}(x(n))+$
$\displaystyle\hskip 19.91692pth_{1}(z(n))(y(n)-g_{y})^{2}\leq P_{\rm
com}(n)\gamma,\forall n\in\mathbb{N}_{NM},$ (43) $\displaystyle\widetilde{\rm
C8a}:h_{1}(z(n))+h_{2}(z(n))\leq t(n),\forall n\in\mathbb{N}_{NM},$ (44)
$\displaystyle\widetilde{\rm C8b}:h_{1}(z(n))+h_{3}(z(n))\leq o(n),\forall
n\in\mathbb{N}_{NM}.$ (45)

Let $\mathbf{z}^{(j)}\in\mathbb{R}^{NM\times 1}$ and
$\mathbf{x}^{(j)}\in\mathbb{R}^{NM\times 1}$ be the collections of all $z(n)$
and $x(n),\forall n$, respectively. Problem $\mathrm{(P.2)}$ can be
approximated around the solution $\\{\mathbf{z}^{(j)},\mathbf{x}^{(j)}\\}$ by
the following convex optimization problem:

$\displaystyle\\!\mathrm{(P.3)}:\max_{\mathbf{z},\mathbf{x},\mathbf{q},\mathbf{P_{\mathrm{sar}}},\mathbf{P_{\mathrm{com}}},\bm{\psi},\mathbf{t},\mathbf{o}}C(\mathbf{u}^{r})$
$\displaystyle\text{s.t.}\hskip
8.53581pt\mathrm{C1-C6,C7a,C7b,C7c,\widetilde{C8},\widetilde{C8a},\widetilde{C8b},C9-C11}.$

Algorithm 1 Successive Convex Approximation (SCA)

1:Initialization: Set iteration index $j=1$, point
$\\{\mathbf{z}^{(1)},\mathbf{x}^{(1)}\\}$, coverage $C^{(1)}$ by solving
$\rm(P.3)$ around point $\\{\mathbf{z}^{(1)},\mathbf{x}^{(1)}\\}$, and error
tolerance $0<\epsilon\ll 1$.

2:repeat

3:Set $j=j+1$ $\triangleright$ Increment iteration index

4:Get coverage $C^{(j)}$ and solution
$\\{\mathbf{z},\mathbf{x},\mathbf{q},\mathbf{P_{\mathrm{sar}}},\mathbf{P_{\mathrm{com}}}\\}$
by solving $\rm(P.3)$ around point
$\\{\mathbf{z}^{(j-1)},\mathbf{x}^{(j-1)}\\}$

5:Set
$\mathbf{z}^{(j)}=\mathbf{z},\mathbf{x}^{(j)}=\mathbf{x}$$\triangleright$
Update solution

6:until $\big{|}\frac{C^{(j)}-C^{(j-1)}}{C^{(j)}}\big{|}\leq\epsilon$
$\triangleright$ Verify convergence

7:return solution
$\\{\mathbf{z}^{r},\mathbf{x}^{r},\mathbf{q},\mathbf{P_{\mathrm{sar}}},\mathbf{P_{\mathrm{com}}}\\}$

To solve $\mathrm{(P.2)}$, we use SCA such that in each iteration $j$ of
Algorithm 1, we update the solution $\\{\mathbf{z}^{(j)},\mathbf{x}^{(j)}\\}$
by solving problem $\mathrm{(P.3)}$ using conventional solvers such as CVX
[37]. To solve problem $\rm(P.1)$, we first solve problem $\rm(P.2)$ for every
fixed and feasible number of drone scans $N$, then we select the minimum
number of azimuth scans that results in the maximum coverage to save drone
energy. In the following proposition, we show that the number of feasible $N$
is finite, and therefore, the exhaustive search is guaranteed to yield the
optimal number of azimuth scans $N$.

######  Proposition 3.

The feasible number of azimuth scans $N$ is finite. The optimal number of
azimuth scans $N^{*}$ is upper bounded by
$\frac{1}{M}\left(\frac{q_{\mathrm{start}}}{\delta_{t}P_{\mathrm{prop}}}+1\right)$.

###### Proof:

Please refer to Appendix B. ∎

Algorithm 1 has a polynomial time complexity [38]. The exhaustive search is of
constant time complexity as the total number of iterations is less than
$\Big{\lceil}\frac{1}{M}\left(\frac{q_{\mathrm{start}}}{\delta_{t}P_{\mathrm{prop}}}+1\right)\Big{\rceil}$
according to Proposition 3. Therefore, the overall time complexity of the
proposed robust scheme is still polynomial. Moreover, Algorithm 1 converges to
a local optimum of problem $\rm(P.1)$ as sca methods are known to converge
fast to locally optimal solutions[39]. Next, to assess the performance of the
proposed robust sub-optimal scheme, we derive an upper bound on the optimal
solution of problem $\rm(P.1)$ using monotonic optimization theory [40].

|
---|---
(a) | (b)
Figure 4:  Lower bound on the uav-gs distance. (a) Proposed lower bound. (b)
Accuracy of the proposed lower bound.

###  IV-C Upper Bound on the Optimal Solution of Problem $\mathrm{(P.1)}$

To provide an upper bound to problem $\mathrm{(P.1)}$, we start by fixing the
total number of azimuth scans $N$, which leads to problem $\mathrm{(P.2)}$.
However, problem $\rm(P.2)$ is non-monotonic w.r.t. optimization variables
$\mathbf{x}$, $\mathbf{z}$, and $\mathbf{P_{\rm com}}$, due to constraints
$\rm C2$ and $\rm C8$. Equality constraints $\rm C3-C5$ also represent an
obstacle to applying monotonic optimization theory. In what follows, we
provide an upper bound on the optimal solution to problem $\rm(P.2)$. To this
end, we introduce an optimization problem, denoted by
$\mathrm{(\overline{P.2})}$, the optimal solution of which presents an upper
bound on the optimal solution to problem $(\mathrm{P.2})$. Then, we use
monotonic optimization theory to provide the optimal solution to problem
$\mathrm{(\overline{P.2})}$.

####  IV-C1 Upper Bound Formulation

During the $M$ time slots of a given azimuth scan, the uav-gs distance
slightly varies depending on the aoi-gs distance. In each azimuth scan, we
impose an upper bound on the uav throughput $R_{n}$ by lower-bounding the uav-
gs distance $d_{n}$. In fact, for a given azimuth scan, we consider the
minimum possible uav-gs distance in the considered scan, see Figure 4a.
Following the proposed lower-bound, the uav-gs distance is constant during a
scan and only varies from scan to scan. This approach allows us to focus on
optimizing the system parameters from scan to scan, instead of optimizing the
system along all time slots, which reduces the problem dimension from $N\times
M$ to $N$. Hereinafter, each $n\in\mathbb{N}_{N}$, refers to the start of a
different azimuth scan. The $y$-positions that correspond to the minimum uav-
gs distance for each azimuth scan, denoted by
$\overline{\mathbf{y}}=[\overline{y}(1),...,\overline{y}(N)]^{T}\in\mathbb{R}^{N\times
1}$, are given by:

$\overline{y}(n)=y(i),i=\operatorname*{arg\,min}\limits_{1+(n-1)M\leq k\leq
nM}(y(k)-g_{y})^{2}.$ (47)

In Figure 4b, we plot the exact uav-gs distance with solid lines and the
proposed lower bound on the uav-gs distance with circle markers. The figure
confirms the tightness of the proposed lower bound for different placements of
the gs. The parameter values used in Figure 4b are specified in Table II. As
each time slot now corresponds to a different azimuth scan, we define vectors
$\overline{\bm{\delta}_{x}^{r}}$ and $\overline{\bm{\delta}_{z}^{r}}$ that
contain the corrections in $x$\- and $z$-direction, respectively, for each of
these scans as follows:

$\displaystyle\overline{\bm{\delta}^{r}_{x}}=\left[\delta_{x}^{r}(1),\delta_{x}^{r}(1+M),...,\delta_{x}^{r}(1+(N-1)M)\right]^{T}\in\mathbb{R}^{N\times
1},$ (48)
$\displaystyle\overline{\bm{\delta}_{z}^{r}}=\left[\delta_{z}^{r}(1),\delta_{z}^{r}(1+M),...,\delta_{z}^{r}(1+(N-1)M)\right]^{T}\in\mathbb{R}^{N\times
1},$ (49)

Based on this approach, the following $N$-dimensional optimization problem
denoted by $\rm(\overline{P.2})$ is obtained as an upper bound for
optimization problem $\rm(P.2)$:

$\displaystyle\\!(\overline{\mathrm{P.2}}):\max_{\mathbf{z},\mathbf{q},\mathbf{P_{\mathrm{sar}}},\mathbf{P_{\mathrm{com}}}}\hskip
8.53581pt\frac{L}{\Delta_{s}}C(\overline{\mathbf{u}^{r}})$ s.t.
$\displaystyle\overline{\mathrm{C6}}:\;z_{\mathrm{min}}\leq\overline{z^{r}}(n)\leq
z_{\mathrm{max}},\forall n\in\mathbb{N}_{N},$
$\displaystyle\overline{\mathrm{C7}}:\mathrm{SNR}_{n}(\mathbf{P}_{\rm
sar},\overline{\mathbf{u}^{r}})\geq\mathrm{SNR_{min}},\forall
n\in\mathbb{N}_{N},$
$\displaystyle\overline{\mathrm{C8}}:R_{n}(\mathbf{P}_{\rm
com},\overline{\mathbf{u}^{r}})\geq
R_{\mathrm{min},n}(\overline{\mathbf{u}^{r}})+R_{\mathrm{sl}},\forall
n\in\mathbb{N}_{N},$ $\displaystyle\overline{\mathrm{C9}}:\;0\leq
P_{\mathrm{sar}}(n)\leq P_{\mathrm{sar}}^{\mathrm{max}},0\leq
P_{\mathrm{com}}(n)\leq P_{\mathrm{com}}^{\mathrm{max}},\forall
n\in\mathbb{N}_{N},$
$\displaystyle\overline{\mathrm{C10}}:q(1)=q_{\mathrm{start}},\;q(n)\geq
0,\forall n\in\mathbb{N}_{N},$
$\displaystyle\overline{\mathrm{C11}}:q(n+1)=q(n)-M\delta_{t}P_{\mathrm{tot}}(n),\forall
n\in{\mathbb{N}_{N-1}},$

where $\overline{z^{r}}(n)=z(n)+\overline{\delta_{z}^{r}}(n)$, and the uav
position in each azimuth scan $n$, denoted by $\overline{\mathbf{u}^{r}}(n)$,
is given by:

$\displaystyle\overline{\mathbf{u}^{r}}(n)=\left[l^{r}_{n}(\mathbf{z}),\overline{y}(n),\overline{z^{r}}(n)\right]^{T},\forall
n\in\mathbb{N}_{N},$ (50)

where
$l^{r}_{n}(\mathbf{z})=(c_{2}-c_{1})\sum\limits_{k=1}^{n-1}z(k)-c_{1}\;z(n)+\overline{\delta_{x}^{r}}(n)$.
Problem $\mathrm{(\overline{P.2})}$ is an $N$-dimensional non-convex non-
monotonic minlp. Next, we show that $(\overline{\mathrm{P.2}})$ is an upper
bound to problem $(\mathrm{P.2})$, then, we provide the optimal solution of
$\rm\overline{(P.2)}$ using monotonic optimization.

######  Theorem 1.

Let $\overline{\mathbf{u}^{r}}$ and $\mathbf{u}^{r}$ be the optimal solutions
to problems $(\overline{\mathrm{P.2}})$ and $\mathrm{(P.2)}$, respectively,
then $\frac{L}{\Delta_{s}}C(\overline{\mathbf{u}^{r}})\geq
C(\mathbf{u}^{r}),\forall r\in[0,1].$ In other words, the optimal solution to
problem $(\overline{\mathrm{P.2}})$ provides an upper bound for the optimal
solution to problem $\mathrm{(P.2)}$.

###### Proof.

Please refer to Appendix C. ∎

####  IV-C2 Reformulation as Monotonic Optimization Problem

To provide the optimal solution to problem $\rm(\overline{P.2})$, we write the
problem in the canonical form of a monotonic optimization problem given by
[41, 42]:

$\max\\{\mathcal{C}(\mathbf{x})|\mathbf{x}\in\mathcal{G}\cap\mathcal{H}\\},$
(51)

where $\mathcal{C}:\mathbb{R}^{N\times 1}_{+}\rightarrow\mathbb{R}$ is the
cost function which is increasing w.r.t. $\mathbf{x}\in\mathbb{R}^{N\times
1}$, $\mathcal{G}\subset[\mathbf{0},\;\mathbf{b}]\subset\mathbb{R}^{N\times
1}_{+}$ is a compact normal set with non-empty interior, and $\mathcal{H}$ is
a conormal set on the box $[\mathbf{0},\;\mathbf{b}]$444Rigorous definitions
of increasing functions, boxes, normal sets, and conormal sets are provided in
the monotonic optimization theory guides [41, 42].. As a first step, we
rewrite constraint $\overline{{\rm C8}}$ as follows:

$\overline{{\rm
C8}}:\left(A2^{\alpha\overline{z^{r}}(n)}-1\right)d^{2}_{n}(\overline{\mathbf{u}^{r}})-P_{\mathrm{com}}(n)\gamma\leq
0,\forall n\in\mathbb{N}_{N}.$ (52)

Constraint $\overline{\mathrm{C8}}$ is non-convex and non-monotonic w.r.t.
optimization variable $\mathbf{z}$ and $\mathbf{P_{\rm com}}$. Thus, we
transform $\overline{\mathrm{C8}}$ into a set of monotonic constraints by
first writing it as a difference of monotonic functions.

######  Proposition 4.

Constraint $\overline{\rm C8}$ is a difference of two monotonic functions
$f(z(n))-g(z(n),P_{\rm com}(n))$, where $f(z(n))$ is given by:

$\displaystyle
f(z(n))=\left(2^{\alpha\overline{z^{r}}(n)}-1\right)\Bigg{(}(c_{2}-c_{1})^{2}\left(\sum_{k=1}^{n-1}z(k)\right)^{2}+2z(n)\Big{(}c_{1}[g_{x}]^{+}-c_{1}\overline{\delta_{x}^{r}}(n)+\overline{\delta_{z}^{r}}(n)\Big{)}+$
$\displaystyle\left(1+{c_{1}}^{2}\right)z^{2}(n)+2(c_{2}-c_{1})\sum_{k=1}^{n-1}z(k)[-g_{x}]^{+}+\left(g_{x}-\overline{\delta_{x}^{r}}(n)\right)^{2}+\left(g_{z}-\overline{\delta_{z}^{r}}(n)\right)^{2}+\left(\overline{y}(n)-g_{y}\right)^{2}\Bigg{)},$
(53)

and $g(z(n),P_{\rm com}(n))$ is given by :

$\displaystyle g(z(n),P_{\rm com}(n))=\gamma P_{\rm
com}(n)+\left(2^{\alpha\overline{z^{r}}(n)}-1\right)\times\Bigg{(}2(c_{1}[-g_{x}]^{+}+g_{z})z(n)+$
$\displaystyle(c_{2}-c_{1})\left(2[g_{x}]^{+}+c_{1}\;z(n)\right)\sum_{k=1}^{n-1}z(k)\Bigg{)},n\in\mathbb{N}_{N}.$
(54)

###### Proof.

The result can be easily obtained by verifying that functions $f$ and $g$ are
increasing wrt $z$ and $f(z(n))-g(z(n),P_{\rm com}(n))$ is equal to the left-
hand side of (52). ∎

Constraint $\overline{\rm C8}$ is a difference of monotonic functions $f$ and
$g$, thereby, as proven in [41], we can equivalently represent it with the
following constraints:

$\displaystyle\overline{\mathrm{C8a}}:f(z(n))+t(n)\leq f(z_{\rm max}),\forall
n\in\mathbb{N}_{N},$ (55) $\displaystyle\overline{\mathrm{C8b}}:g(z(n),P_{\rm
com}(n))+t(n)\geq f(z_{\rm max}),\forall n\in\mathbb{N}_{N},$ (56)
$\displaystyle\overline{\mathrm{C8c}}:0\leq t(n)\leq f(z_{\rm
max})-f(0),\forall n\in\mathbb{N}_{N},$ (57)

where $\mathbf{t}=[t(1),...,t(N)]^{T}\in\mathbb{R}^{N\times 1}$ is a vector of
auxiliary optimization variables. Problem $(\overline{\rm P.2})$ can now be
written in the canonical form of a monotonic optimization problem as follows:

$\displaystyle\\!(\overline{\rm
P.3}):\max_{\mathbf{z},\mathbf{t},\mathbf{q},\mathbf{P_{\mathrm{sar}}},\mathbf{P_{\mathrm{com}}}}\frac{L}{\Delta_{s}}C(\overline{\mathbf{u}^{r}})$
$\displaystyle\text{s.t.}\hskip
8.53581pt\mathrm{\overline{C6},\overline{C7},\overline{C8a},\overline{C8b},\overline{C8c},\overline{C9},\mathrm{\overline{C10}},\overline{C11}},$

where the feasible set is the intersection of normal set $\mathcal{G}$ that is
spanned by constraints ($\rm\overline{C6}$, $\rm\overline{C7}$,
$\rm\overline{C8a}$, $\rm\overline{C8c}$, $\rm\overline{C9}$,
$\rm\overline{C10}$, $\rm\overline{C11}$) and conormal set $\mathcal{H}$ that
is spanned by constraints $\rm(\overline{C6},\overline{C8b})$. The optimal
solution of problem $(\overline{\rm P.3})$ lies at the upper boundary of the
normal set $\mathcal{G}$, denoted by $\partial^{+}\mathcal{G}$ [42]. This
region is not known a priori but the sequential polyblock approximation
algorithm can be used to approach this set from above.

####  IV-C3 Polyblock Approximation Algorithm

Algorithm 2 Polyblock Outer Approximation Algorithm

1:Initialization: Set polyblock $\mathcal{P}_{1}$ with vertex $\mathbf{v}_{\rm
max}=[\mathbf{z}_{\rm max},\mathbf{t}_{\rm max},\mathbf{P_{\rm com}^{\rm
max}}]^{T}$, vertex set $\mathcal{T}^{(1)}=\\{\mathbf{v}_{\rm max}\\}$, small
positive number $\epsilon\geq 0$, objective function value $\rm
CBV_{0}=-\infty$, and iteration number $k=0$.

2:repeat

3: $k=k+1$.$\triangleright$ Increment iteration index

4: Select
$\mathbf{v}_{k}=\operatorname*{arg\,max}_{\mathbf{v}\in\mathcal{T}^{(k)}}(C(\mathbf{v}))$$\triangleright$
Select the best vertex in set $\mathcal{T}^{(k)}$

5: Compute the projection $\Phi(\mathbf{v}_{k})$ on $\partial\mathcal{G}^{+}$
using Algorithm 3

6: if $\Phi(\mathbf{v}_{k})=\mathbf{v}_{k}$, i.e.,
$\mathbf{v}_{k}\in\mathcal{G}$ then

7: $\mathbf{s}_{k}=\mathbf{v}_{k}$ and
$\mathrm{CBV}_{k}=C(\Phi(\mathbf{v}_{k}))$

8: else

9: if $\Phi(\mathbf{v}_{k})\in\mathcal{G}\cap\mathcal{H}$ and
$C(\Phi(\mathbf{v}_{k}))\geq\mathrm{CBV}_{k-1}$ then $\triangleright$ Better
solution is found

10: Let the current best solution $\textbf{s}_{k}=\Phi(\mathbf{v}_{k})$ and
$\mathrm{CBV}_{k}=C(\Phi(\mathbf{v}_{k}))$ $\triangleright$ Update best
solution

11: else

12: $\textbf{s}_{k}=\textbf{s}_{k-1}$ and
$\mathrm{CBV}_{k}=\mathrm{CBV}_{k-1}$

13: endif

14: Set
$\mathcal{T}^{(k+1)}=\\{\mathcal{T}^{(k)}\setminus\mathcal{T}^{*}\\}\cup\\{\mathbf{v}^{i}=\mathbf{v}+({(\Phi(\mathbf{v}_{k}))}_{i}-v_{i})\mathbf{e}_{i}|\mathbf{v}\in\mathcal{T}^{*},i\in\\{1,..,3N\\}\\}$,

15: where
$\mathcal{T}^{*}=\\{\mathbf{v}\in\mathcal{T}^{(k)}|\mathbf{v}>\Phi(\mathbf{v}_{k})\\}$

16: Remove from $\mathcal{T}^{(k)}$ improper vertices and vertices
$\\{\mathbf{v}\in\mathcal{T}^{(k+1)}|\mathbf{v}\notin\mathcal{H}\\}$
$\triangleright$ Discard unnecessary vertices from set $\mathcal{T}^{(k+1)}$

17: endif

18: until $|C(\mathbf{v}_{k})-\mathrm{CBV}_{k}|\leq\epsilon$$\triangleright$
Verify convergence

19: return solution $\mathbf{s}_{k}$

In Algorithm 2, we provide all steps of the proposed polyblock outer
approximation algorithm that generates an $\epsilon$-approximate optimal
solution to problem $\rm(\overline{P.3})$. In addition, Figure 5 depicts the
operation of the polyblock algorithm where, for simplicity of illustration, we
consider a case with only two dimensions $x$ and $y$. The optimization
variables of problem $(\overline{\rm P.3})$ are divided into outer
optimization variables, that are collected in vertex vector
$\mathbf{v}=\big{[}{\overline{\mathbf{z}^{r}}\;}^{T}$, $\mathbf{t}^{T}$,
${(\mathbf{P}_{\rm com})}^{T}\big{]}^{T}\in\mathbb{R}^{3N\times 1}$, and inner
optimization variables, $\mathbf{q}$ and $\mathbf{P_{\mathrm{sar}}}$.
Hereinafter, as the objective function depends only on $\mathbf{z}$, for ease
of notation, $C(\mathbf{v})$ refers to
$\frac{L}{\Delta_{s}}C([v(1),...,v(N)]^{T})$.

|  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 5: Illustration of polyblock outer approximation algorithm. (a)
Feasible set for a monotonic maximization problem where the normal and
conormal sets are denoted by $\mathcal{G}$ and $\mathcal{H}$, respectively.
(b) Projection of $\mathbf{v}_{\rm max}$ on normal set $\mathcal{G}$. (c)
Projection of the best vertex
$\mathbf{v}^{1}=\operatorname*{arg\,max}\limits_{\mathbf{v}\in\mathcal{T}^{(1)}}\\{C(\mathbf{v})\\}$
on normal set $\mathcal{G}$ and removal of improper vertex $\mathbf{v}^{4}$
(dominated by $\mathbf{v}^{3}$). (d) Representation of a tighter polyblock
obtained after several iterations.

As shown in Figure 5a, first, we initialize a polyblock $\mathcal{P}^{(1)}$
with corresponding vertex set $\mathcal{T}^{(1)}$, which initially contains
vertex $\mathbf{v}_{\rm max}=[\mathbf{z}_{\rm max},\mathbf{t}_{\rm
max},\mathbf{P_{\rm com}^{\rm max}}]^{T}\in\mathbb{R}^{3N\times 1}$, where
$\mathbf{z}_{\rm max}=[z_{\rm max}-\overline{\delta_{z}^{r}}(1),...,z_{\rm
max}-\overline{\delta_{z}^{r}}(N)]\in\mathbb{R}^{1\times N}$, $\mathbf{t}_{\rm
max}=[f(z_{\rm max})-f(0),...,f(z_{\rm max})-f(0)]\in\mathbb{R}^{1\times N}$,
and $\mathbf{P_{\rm com}^{\rm max}}=[P_{\rm com}^{\rm max},...,P_{\rm
com}^{\rm max}]\in\mathbb{R}^{1\times N}$. The projection of the vertex
$\mathbf{v}_{\rm max}$ onto the normal set $\mathcal{G}$, denoted by
$\Phi(\mathbf{v}_{\rm max})\in\mathbb{R}^{3N\times 1}$, is calculated and used
to construct a tighter polyblock $\mathcal{P}^{(2)}$, see Figure 5b. The new
vertex set is calculated as
$\mathcal{T}^{(2)}=\\{\mathcal{T}^{(1)}\setminus\mathbf{v}_{\rm
max}\\}\cup\\{\mathbf{v}^{i}=\mathbf{v}_{\rm
max}+\left({\left(\Phi(\mathbf{v}_{\rm max})\right)}_{i}-{(\mathbf{v}_{\rm
max})}_{i}\right)\mathbf{e}_{i},i\in\\{1,...,3N\\}\\}$, where
${\left(\Phi(\mathbf{v}_{\rm max})\right)}_{i}$ and ${(\mathbf{v}_{\rm
max})}_{i}$ are the $i^{th}$-elements of vectors $\Phi(\mathbf{v}^{\rm max})$
and $\mathbf{v}_{\rm max}$, respectively. Here, $\mathbf{e}_{i}$ is the unit
vector whose $i^{th}$ entry is equal to $1$ and all other entries are equal to
$0$. Then, in iteration $k$, the vertex that maximizes the objective function
of problem $\mathrm{(\overline{P.3})}$ is selected to create the next
polyblock, i.e.,
$\mathbf{v}_{k}=\operatorname*{arg\,max}\limits_{\mathbf{v}\in\mathcal{T}^{(k)}}\\{C(\mathbf{v})\\}$.
Once a feasible projection (i.e.,
$\Phi(\mathbf{v}_{k})\in\mathcal{G}\cap\mathcal{H}$) that improves the
objective function of problem $\rm(\overline{P.3})$ is located, the
corresponding vector $\mathbf{s}_{k}=\Phi(\mathbf{v}_{k})$ and its objective
function $\mathrm{CBV}_{k}=C(\Phi(\mathbf{v}_{k}))$ are recorded. To speed up
convergence and save memory, improper vertices555Given a vertex set
$\mathcal{T},$ improper vertices are vertices that are dominated by other
vertices in $\mathcal{T}.$ See Figure 5c and [41] for further details. and
vertices that are not in set $\mathcal{H}$ are removed in each iteration,
e.g., vertex $\textbf{v}^{4}$ in Figure 5c. The same procedure is repeated
multiple times, such that a sequence of polyblocks converges from above to the
feasible set, i.e.,
$\mathcal{P}^{(1)}\supset\mathcal{P}^{(2)}\supset...\supset\mathcal{P}^{(k)}\supset\mathcal{G}\cap\mathcal{H}$,
see Figure 5d. The termination criterion of the proposed algorithm is
$|C(\mathbf{v}_{k})-\mathrm{CBV}_{k}|\leq\epsilon$, where $\epsilon$ is the
error tolerance.

Algorithm 3 Bisection Search to Compute the Upper Boundary Point

1:Initialization: Set vertex $\mathbf{v}$, $\lambda_{\rm min}=0$,
$\lambda_{\rm max}=1$, and error tolerance $0<\epsilon\ll 1$.

2:repeat

3: $\lambda=\frac{(\lambda_{\rm max}+\lambda_{\rm min})}{2}$ $\triangleright$
Bisection middle point

4: if $\lambda\mathbf{v}$ is feasible for problem $\mathrm{(\overline{P.3})}$
then

5: Set $\lambda_{\rm min}=\lambda$ $\triangleright$ Search in
[$\lambda,\lambda_{\rm max}$]

6: else

7: Set $\lambda_{\rm max}=\lambda$$\triangleright$ Search in [$\lambda_{\rm
min},\lambda$]

8: end if

9:until $\lambda_{\rm max}-\lambda_{\rm min}\leq\epsilon$$\triangleright$
Verify convergence

10:return $\Phi(\mathbf{v})=\lambda\mathbf{v}$

In Line $4$ of Algorithm 2, we compute the projection of vertex
$\mathbf{v}_{k}=\operatorname*{arg\,max}\limits_{\mathbf{v}\in\mathcal{T}^{(k)}}C(\mathbf{v})$
on $\partial\mathcal{G}^{+}$, denoted by
$\Phi(\mathbf{v}_{k})\in\mathbb{R}^{3N\times 1}$. This projection is not
trivial and is equivalent to solving a one-dimensional optimization problem:
$\max\limits_{0\leq\lambda\leq 1}\\{\lambda\mathbf{v}_{k}\in\mathcal{G}\\}$.
To compute $\Phi(\mathbf{v}_{k})$, we perform a bisection search using
Algorithm 3. In each iteration of Algorithm 3, checking if a vector belongs to
the normal set $\mathcal{G}$ is equivalent to solving problem $(\overline{\rm
P.3})$ when fixing the outer optimization variables, i.e., vectors
$\mathbf{z}$, $\mathbf{t}$, and $\mathbf{P}_{\rm com}$. This results in the
following feasibility problem:

$\displaystyle\\!\mathrm{(\overline{P.4})}:\max_{\mathbf{q},\mathbf{P_{\mathrm{sar}},}}\hskip
8.53581pt1$ $\displaystyle\text{s.t.}\hskip
8.53581pt\mathrm{\overline{C7},\overline{C9},\overline{C10},\overline{C11}}.$

Problem $\mathrm{(\overline{P.4})}$ is a convex optimization problem and a
feasible solution can be obtained by, e.g., CVX [37]. Once the projection
point is computed, the new polyblock is generated as described previously and
these steps are repeated until the optimal solution is located.

The computational complexity of the upper bound on the optimal solution,
provided by Algorithm 2 and Algorithm 3, is exponential in the number of
vertex elements, i.e., the time complexity order is $O(2^{3N})$ [42]. If
problem $\rm\overline{(P.2)}$ is feasible, then Algorithm 2 converges to its
global optimal solution as the polyblock approximation algorithm has been
proven to converge to an $\epsilon$-optimal solution [42]. Despite the high
complexity of the derived upper bound, it can be used to assess the
performance of the proposed low-complexity sub-optimal scheme.

TABLE II:  System parameters [16, 43, 34]. Parameter | Value  | Parameter | Value
---|---|---|---
$L$ | 60 m  | $\mathbf{g}=[g_{x},g_{y},g_{z}]^{T}$ |  $[0,0,5]^{T}$ m
$M$ | 100  | $v$ | 5 m/s
$z_{\rm max}$ and $z_{\rm min}$ | 100 m and 2 m  | $q_{\rm start}$ | 19.44 Wh
$P_{\rm sar}^{\rm max}$ | 46 dBm  | $P_{\rm prop}$ | 450 W
$P_{\rm com}^{\rm max}$ | 40 dBm  | $P_{0}$ | 79.86 W
$\theta_{d}$ and $\theta_{\rm 3db}$ | 45o and 30o | $P_{I}$ | 420.6 W
$\tau_{p}$ and $\mathrm{PRF}$ | 1 $\mu$s and 100 Hz  | $W_{u}$ | 56.5 N
$B_{r}$ | 100 MHz  | $U_{\mathrm{tip}}$ | 120 m/s
$B_{c}$ | 100 MHz  | $\rho$ | 1.225 kg/m3
$f$ | 2 GHz  | $A$ | 0.503 m2
$\mathrm{SNR}_{\rm min}$ | 20 dB  | $d_{0}$ | 0.6
$R_{\rm sl}$ | 1 kbit  | $s$ | 0.05
$\gamma$ | 20 dB  | $\beta$ |  $10^{4}$ W-1 m3
|  |  |

##  V Simulation Results and Discussion

In this section, we present simulation results for the proposed robust uav 3d
trajectory and resource allocation algorithm, where we adopted the parameter
values provided in Table II, unless specified otherwise. The minimum uav
altitude $z_{\rm min}$ was set to 2 m to avoid crashes due to collisions with
obstacles at very low altitudes. For all robust schemes, we adopt the Gaussian
model for the uav trajectory deviations where in-range position errors follow
$\mathcal{N}(o_{x}=1\;\rm{m},\sigma=0.3\;\rm{m})$ and altitude errors follow
$\mathcal{N}(o_{z}=-1\;\rm{m},\sigma=0.3\;\rm{m})$.

###  V-A Benchmark Schemes

To evaluate the performance of the proposed robust scheme, we compare it with
the following benchmark schemes:

  * •

Benchmark scheme 1: This is a non-robust scheme, where we solve problem
$\rm(P.1)$ assuming no uav trajectory deviations, i.e., $\sigma=0$, $o_{x}=0$,
and $o_{z}=0$. This benchmark scheme was presented in [1] where
$\mathbf{u}^{r}=\mathbf{u}$.

  * •

Benchmark scheme 2: This is a robust scheme, where we solve problem $\rm(P.1)$
assuming the communication power $P_{\rm com}$ is fixed across all time slots
and optimized along with the other optimization variables. The coverage
reliability level is $r=0.95$.

  * •

Benchmark scheme 3: This is also a robust scheme, with $r=0.95$, where we
solve problem $\rm(P.1)$ assuming the sar power $P_{\rm sar}$ is fixed across
all time slots but optimized along with the other optimization variables.

We note that all benchmark schemes have polynomial time complexity, similar to
the proposed sub-optimal scheme.

###  V-B Convergence and Accuracy of the Proposed Scheme

|
---|---
(a) | (b)
Figure 6: Convergence of algorithm 1 (sub-optimal) and algorithm 2 (upper
bound) for $N=1$ (a) and $N=3$ (b).  Figure 7:  Comparison of the radar
coverage of the considered schemes. Figure 8: Total missed area of the
proposed robust scheme compared to the non-robust scheme for different
trajectory deviation levels. Shaded area represents one standard deviation
while solid line represents the mean value averaged over $10^{4}$
realizations.

In Figure 6, we analyze the convergence rate of the proposed sub-optimal and
optimal (upper bound) schemes for different numbers of azimuth sweeps $N$. We
also assess the gap between the optimal and sub-optimal solutions. The results
are plotted as a function of the iteration number of Algorithm 1 and Algorithm
2, respectively. Figure 6 shows that the proposed scheme converges fast to the
upper bound on the optimal solution given by $C(z^{*})\approx 9.154\times
10^{3}$ m2 for $N=1$, and $C(z^{*})\approx 2.496\times 10^{4}$ m2 for $N=3$.
After convergence, the gap between the upper bound on the optimal solution and
the sub-optimal solution is negligible which confirms that the proposed sub-
optimal scheme achieves close-to-optimal performance with low complexity.

In Figure 7, we depict the sensing coverage as a function of the number of
azimuth sweeps $N$. The figure illustrates how the optimal $N^{*}$, leading to
the best possible coverage, can be obtained based on a finite exhaustive
search. Furthermore, we compare the proposed scheme with the benchmark
schemes. In fact, the figure illustrates the superior performance of the
proposed robust scheme over robust benchmark schemes 2 and 3. The gain is
attributed to the optimization of the communication power and radar transmit
power, respectively. In comparison to non-robust benchmark scheme 1, our
proposed robust scheme achieves comparable or improved coverage depending on
the system parameters. Yet, the advantage of the proposed robust scheme over
the non-robust scheme is illustrated in Figure 8, which shows the total area
missed by the uav (i.e., the total area of all holes in the total coverage)
averaged over $10^{4}$ realizations and for different levels of random
trajectory deviations. For the non-robust design, the missed area increases
with the number of total azimuth sweeps, as more area is missed by the uav-sar
due to random trajectory deviations. However, the proposed robust design keeps
the missed area small and constant for different numbers of azimuth sweeps to
ensure that the required $r=95\%$ coverage reliability is achieved. This
confirms that, in terms of total missed area, the robust scheme outperforms
the non-robust scheme, as the latter may cause significant coverage holes and
thus, is not suitable in practice.

###  V-C Robustness of the Proposed Design Against uav Trajectory Deviations

c  |  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 9:  Robust versus non-robust uav trajectory design. (a) Non-robust 3d
uav trajectory for optimized coverage. (b) Actual 2d ground sar coverage
following the trajectory shown in (a). (c) Robust uav trajectory for optimized
coverage. (d) Actual 2d ground sar coverage following the trajectory shown in
(c).

In Figure 9, we investigate the effect of uav trajectory deviations on sar
mapping for low-altitude radar applications (i.e., $z_{\rm max}=10$ m). In
Figure 9a, the dashed red line denotes the non-robust uav trajectory provided
by benchmark scheme 1, whereas the solid blue line shows the actual uav flight
path. In Figure 9b, we show the resulting uav-sar ground coverage with green-
and yellow-colored rectangles denoting successive azimuth scans. Following an
ideal path planning, the expected covered area is measured by the dashed red
rectangle. However, the actual covered area is different from the expected
result. Figure 9b shows how small flight deviations can affect the uav
coverage, where the total uncovered area for benchmark scheme 1 was $225$ m2
out of $C(\mathbf{z})=1.859\times 10^{3}$ m2 (12%). In Figure 9c, we present
the proposed robust uav trajectory with a solid green line. Figure 9d
highlights the tradeoff between coverage performance and robustness for low-
altitude sensing applications. Although the total coverage is comparable, the
proposed robust scheme is highly effective in avoiding uncovered areas caused
by random uav trajectory deviations, as any location in the aoi is guaranteed
to be covered with a 95% probability.

Figure 10: Total radar coverage versus maximum communication power.

###  V-D Real-time Communication

In Figure 10, we plot the total radar ground coverage as a function of the
uav’s maximum communication power for different placements of the gs and
different required radar $\mathrm{SNR}$ values. This figure highlights the
impact of communication on the sensing performance. In particular, the
communication transmit power affects the uav trajectory (see (16)), and
therefore also the sensing performance, as the trajectory of the drone is
linked to the ground coverage and the sensing snr through (11) and (14),
respectively. Figure 10 shows that, when the gs is far from the aoi, i.e., for
$\mathbf{g}=[-120,0,25]$ m, the total coverage is reduced as the minimum
required data rate cannot be achieved at distances far from the gs due to the
increased path loss. This can be compensated by increasing the maximum uav
transmit power. In fact, for sufficiently high maximum uav transmit powers,
the radar snr constraint becomes the performance bottleneck causing the total
coverage to saturate. For instance, as shown in Figure 10, a maximum uav
communication power of $P_{\rm com}^{\rm max}=7.78$ dBW results in the best
possible radar coverage when $\rm SNR_{\rm min}=30$ dB and
$\mathbf{g}=[-120,0,25]$ m. In a nutshell, a reliable data backhaul link is
crucial for enabling real-time on-ground processing of sar data.

###  V-E Impact of UAV Battery Capacity on Sensing Performance

Figure 11:  Total radar coverage versus UAV battery capacity for different
sensing requirements.

In Figure 11, we depict the sar coverage as a function of the uav battery
capacity for different required sensing snr, $\mathrm{SNR}_{\mathrm{min}}$.
The figure shows that the optimum sar coverage increases with the battery
capacity and saturates for a certain value denoted by
$q_{\mathrm{start}}^{*}$, which is the minimum uav battery capacity needed to
achieve maximum coverage for a given required sensing snr. The coverage can
not be further improved by increasing $q_{\mathrm{start}}$ beyond
$q_{\mathrm{start}}^{*}$ due to communication constraint $\mathrm{C8}$, which
eventually becomes the performance bottleneck for large uav-gs distances. In
addition, Figure 11 reveals that for more stringent sensing requirements, not
only does the total coverage decrease, but the energy needed to achieve
maximum coverage also increases. This is because the flying altitude has to be
reduced to satisfy more stringent snr requirements leading to less coverage.
However, when the uav is forced to fly at lower altitudes, a larger number of
sweeps have to be performed to maximize the coverage, which results in a
higher energy consumption. For example, for $\mathrm{SNR}_{\mathrm{min}}$ = 10
dB, $C(z^{*})=5.9\times 10^{4}$ m2 is obtained with a 22 Wh battery capacity,
whereas for $\mathrm{SNR}_{\mathrm{min}}$ = 30 dB, only $C(z^{*})=3.8\times
10^{4}$ m2 is achieved requiring almost twice the battery capacity, i.e., 41
Wh. This further shows that the limited energy available to small uav, which
is accounted for by by constraints $\mathrm{C9-C11}$, has a significant impact
on sensing performance.

##  VI Conclusion

In this paper, we investigated robust joint 3d trajectory and resource
allocation design for uav-sar systems, where real-time transmission of radar
data to the ground was taken into account. Random uav trajectory deviations
were modeled and their effect on the ground radar coverage was derived based
on their statistics. Moreover, the robust 3d trajectory and resource
allocation design was formulated as a mixed-integer non-linear optimization
problem for maximization of the uav radar coverage. A low-complexity sub-
optimal solution based on sca was provided. To assess the performance of the
proposed algorithm, we derived an upper bound for the optimal solution based
on monotonic optimization theory. Simulation results confirmed that the
proposed robust sub-optimal solution achieves close-to-optimal performance.
The superior performance of the proposed scheme compared to three benchmark
schemes was established and the robustness of the proposed solution wrt uav
trajectory deviations was confirmed. Finally, we highlighted the impact of gs
placement and the power available for communication on uav radar performance.
In this paper, we assumed a deterministic los link between the uav and the gs.
An interesting topic for future work is the extension to other types of air-
to-ground channels, including Ricean fading channels, based on online
optimization techniques, see, e.g., [39]. Furthermore, while we used classical
optimization methods and focused on single-uav sar sensing, the extension to
multi-uav sar sensing [44] and the exploration of learning-based methods for
the trajectory design and resource allocation are promising research
directions.

##  Appendix A Proof of Proposition 1

Assuming Gaussian uav position deviations, i.e.,
$\Delta_{x}(n)\sim\mathcal{N}(o_{x},\sigma)$ and
$\Delta_{z}(n)\sim\mathcal{N}(o_{z},\sigma)$, and using (20) and (21), we can
show that
$\Delta_{\mathbf{p}_{N}}(n)\sim\mathcal{N}\big{(}o_{x}+c_{1}o_{z},\sigma\sqrt{1+c_{1}^{2}}\big{)}$
and
$\Delta_{\mathbf{p}_{F}}(n)\sim\mathcal{N}\big{(}o_{x}+c_{2}o_{z},\sigma\sqrt{1+c_{2}^{2}}\big{)},\forall
n\in\mathbb{N}_{NM}$. In time slot $n$, let set
$\mathcal{I}_{N}=\\{s\in\mathbb{R}\hskip 5.69054pt\big{|}\hskip
5.69054pt\mathbb{P}\left(\Delta_{\mathbf{p}_{N}}(n)+s\leq 0\right)\geq r\\}$,
where $r$ is the coverage reliability. Based on the definition of the near-
range and far-range compensation provided in (22) and (23), respectively, we
can write:

$\displaystyle F_{\Delta_{\mathbf{p}_{N}}(n)}(-s)$ $\displaystyle\geq r,$ (58)
$\displaystyle 1-F_{-\Delta_{\mathbf{p}_{N}}(n)}(s)$ $\displaystyle\geq r,$
(59)
$\displaystyle\frac{1}{2}\left(1-\mathrm{erf}\left(\frac{s+o_{x}+c_{1}o_{z}}{\sigma\sqrt{2(1+c_{1}^{2}})}\right)\right)$
$\displaystyle\geq r,$ (60)
$\displaystyle\mathrm{erf}\left(\frac{s+o_{x}+c_{1}o_{z}}{\sigma\sqrt{2(1+c_{1}^{2}})}\right)$
$\displaystyle\leq 1-2r.$ (61)

Since the inverse of the error function, denoted by $\rm erf^{-1}(\cdot)$, is
an increasing function, we have:

$\displaystyle s\leq h_{\mathbf{p}_{N}}(r),$ (62)

where
$h_{\mathbf{p}_{N}}(r)=\mathrm{erf}^{-1}(1-2r)\sigma\sqrt{2(1+c_{1}^{2}})-o_{x}-c_{1}o_{z}$.
Then, the required near-range compensation
$\delta_{\mathbf{p}_{N}}^{r}(n)=\operatorname*{arg\,min}(|s|)$ is given by:

$\displaystyle\delta_{\mathbf{p}_{N}}^{r}(n)=\begin{cases}0&\mathrm{if}\hskip
5.69054pth_{\mathbf{p}_{N}}(r)\geq 0,\\\
h_{\mathbf{p}_{N}}(r)&\mathrm{if}\hskip 5.69054pth_{\mathbf{p}_{N}}(r)\leq
0,\end{cases}$ (63)

which can be equivalently written as follows:

$\displaystyle\delta_{\mathbf{p}_{N}}^{r}(n)=\min\big{(}0,h_{\mathbf{p}_{N}}(r)\big{)}=-[-h_{\mathbf{p}_{N}}(r)]^{+}.$
(64)

This gives the final result in (24). The far-range compensation, denoted by
$\delta_{\mathbf{p}_{N}}^{r}(n)$, is derived using the same methodology.

##  Appendix B Proof of Proposition 3

To show that the set of feasible azimuth scans is finite, we prove that $N$ is
upper bounded. From constraints $\rm C9$ and $\rm C11$, we obtain the
following inequality:

$\delta_{t}P_{\mathrm{prop}}\leq q(n)-q(n+1){,\forall n\in\mathbb{N}_{NM-1}.}$
(65)

The summation of both sides of (65) yields:

$\sum_{n=1}^{NM-1}\delta_{t}P_{\mathrm{prop}}\leq q(1)-q(NM)\leq q_{\rm
start}.$ (66)

The result in (66) can be rigorously proven using induction applied to the
telescoping series $q(n)-q(n+1)$ and constraint $\rm C10$. Since the drone
flies at a constant velocity $v$ and $P_{\rm prop}$ is constant, we obtain an
upper bound on $N$ as follows:

$N\leq\frac{1}{M}\left(\frac{q_{\mathrm{start}}}{\delta_{t}P_{\mathrm{prop}}}+1\right).$
(67)

Problem $\rm(P.1)$ is infeasible for $N$ exceeding the upper bound provided in
(67) as this would violate constraint $\rm C11$.

##  Appendix C Proof of Theorem 1

In this proof, we show that problem $\rm(\overline{P.2})$ is an upper bound to
problem $\rm(P.2)$. To this end, we introduce an intermediate problem, denoted
by $\rm(P.2^{\prime})$, that provides an upper bound for problem $\rm(P.2)$,
then we prove that problems $\rm(P.2^{\prime})$ and $\rm(\overline{P.2})$ are
equivalent. This means that $\rm(\overline{P.2})$ is an upper bound to the
original problem $\rm(P.2)$. We start by rewriting constraint $\rm C8$ as
follows:

$\displaystyle{\rm C8}:g_{1}(n)\leq-g_{2}(n)(y(n)-g_{y})^{2},\forall
n\in\mathbb{N}_{NM},$ (68)

where functions $g_{1}$ and $g_{2}$ are respectively given by:

$\displaystyle g_{1}(n)$
$\displaystyle=g_{2}(n)\Big{(}(x^{r}(n)-g_{x})^{2}+(z^{r}(n)-g_{z})^{2}\Big{)}-P_{\mathrm{com}}(n)\gamma,$
(69) $\displaystyle g_{2}(n)$ $\displaystyle=A2^{\alpha z^{r}(n)}-1,\forall
n\in\mathbb{N}_{NM}.$ (70)

As proposed in (47), we lower bound the term $(y(n)-g_{y})^{2}$ in (68) by
taking its minimum value for each azimuth scan. Since $g_{2}(n)\geq 0,\forall
n\in\mathbb{N}_{NM}$, by bounding $(y(n)-g_{y})^{2}$, the right-hand side
value in (68) increases. This is equivalent to relaxing constraint $\rm C8$,
which we denote now by $\rm\overline{C8}$:

$\displaystyle{\rm\overline{C8}}:g_{1}(n)\leq-
g_{2}(n)\min\limits_{1+(n-1)M\leq k\leq nM}(y(k)-g_{y})^{2},\forall
n\in\mathbb{N}_{N}.$ (71)

Let us denote problem $\rm(P.2)$ after relaxing constraint $\rm C8$ by
$\rm(P.2^{\prime})$. It is easy to show that $\rm(P.2^{\prime})$ provides an
upper bound to problem $\rm(P.2)$. In fact, any feasible solution to
$\rm(P.2)$ is also feasible for $\rm(P.2^{\prime})$ as the solution
necessarily satisfies ${\rm\overline{C8}}$. Now, let us denote the optimal
altitudes to problem $\rm(P.2)$ and $\rm(P.2^{\prime})$ by
$\mathbf{z}\in\mathbb{R}^{NM\times 1}$ and
$\mathbf{z}_{2}\in\mathbb{R}^{NM\times 1}$, respectively, and their robust
trajectories by $\mathbf{u}^{r}$ and $\mathbf{u}_{2}^{r}$, respectively. Since
the objective function is increasing wrt altitude vector $\mathbf{z}$, then,
relaxing $\rm C8$ results in the following inequality:

$C(\mathbf{u}^{r})\leq C(\mathbf{u}^{r}_{2}),\forall r\in[0,1].$ (72)

In the remainder of this proof, we show that problems $\rm(P.2^{\prime})$ and
$(\overline{\mathrm{P.2}})$ are equivalent.

We apply a series of transformations to transform $NM$-dimensional problem
$\rm(P.2^{\prime})$ to the $N$-dimensional problem
$(\overline{\mathrm{P.2}})$. In fact, based on $\mathrm{\overline{C8}}$, where
the uav-gs distance varies from scan to scan, we remark that all optimization
variables of problem $\rm(P.2^{\prime})$ are also constant for a given azimuth
scan and only change from scan to scan. This allows us to reduce the dimension
of the problem as optimization along the azimuth direction has no impact.
Therefore, constraints $\rm C3-C5$ are dropped as they define equalities for
time slots of the same azimuth scans. The rest of the constraints are now
specified for each scan and denoted by $\overline{\rm C6}-\overline{\rm C11}$.
We use $N$-dimensional optimization variables, where based on
${\rm\overline{C8}}$ and constraints $\rm C1$ and $\rm C2$, optimization
vector $\overline{\mathbf{x}}\in\mathbb{R}^{N\times 1}$ can be dropped from
the optimization problem $(\mathrm{P.2^{\prime}})$ as it is deterministic and
can be calculated from vector $\overline{\mathbf{z}}\in\mathbb{R}^{N\times 1}$
as follows:

$\displaystyle\overline{x}(n)=(c_{2}-c_{1})\sum_{k=1}^{n-1}\overline{z}(k)-c_{1}\overline{z}(n),n\in\mathbb{N}_{N}.$
(73)

Finally, we obtain the following relation between the objective functions of
problems $(\mathrm{P.2^{\prime}})$ and $(\overline{\mathrm{P.2}})$:

$\displaystyle C(\mathbf{u}_{2}^{r})$
$\displaystyle=\sum_{n=1}^{NM}\Delta_{s}\left(z_{2}(n)+\delta^{r}_{z}(n)\right)\left(c_{2}-c_{1}\right)$
(74)
$\displaystyle\stackrel{{\scriptstyle\mathrm{C4}}}{{=}}L\sum_{k=1}^{N}\Big{(}z_{2}(1+(1-k)M)+\delta^{r}_{z}(1+(1-k)M)\Big{)}\left(c_{2}-c_{1}\right)$
(75)
$\displaystyle=\frac{L}{\Delta_{s}}\sum_{n=1}^{N}\Delta_{s}\overline{z^{r}}(n)\left(c_{2}-c_{1}\right)=\frac{L}{\Delta_{s}}C(\overline{\mathbf{u}^{r}}).$
(76)

Based on (72) and (76), problem $\rm(P.2^{\prime})$ is equivalently
transformed to problem $\rm(\overline{P.2})$, with $C(\mathbf{u}^{r})\leq
C(\mathbf{u}_{2}^{r})=\frac{L}{\Delta_{s}}C(\overline{\mathbf{u}^{r}})$, which
concludes the proof.

## References

  * [1] M.-A. Lahmeri, W. Ghanem, C. Knill, and R. Schober, “Trajectory and resource optimization for UAV synthetic aperture radar,” in _Proc. IEEE Global Commun. Conf._ , Dec. 2022, pp. 897–903.
  * [2] A. Moreira, P. Prats-Iraola, M. Younis, G. Krieger, I. Hajnsek, and K. P. Papathanassiou, “A tutorial on synthetic aperture radar,” _IEEE Geosci. Remote Sens. Mag._ , vol. 1, no. 1, pp. 6–43, Mar. 2013.
  * [3] D. C. Munson, J. D. O’Brien, and W. K. Jenkins, “A tomographic formulation of spotlight-mode synthetic aperture radar,” _Proc. IEEE_ , vol. 71, no. 8, pp. 917–925, Aug. 1983.
  * [4] L. Ulander, H. Hellsten, and G. Stenstrom, “Synthetic-aperture radar processing using fast factorized back-projection,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 39, no. 3, pp. 760–776, Jul. 2003.
  * [5] P. A. Rosen, S. Hensley, I. R. Joughin, F. K. Li, S. N. Madsen, E. Rodriguez, and R. M. Goldstein, “Synthetic aperture radar interferometry,” _Proc. IEEE_ , vol. 88, no. 3, pp. 333–382, Mar. 2000.
  * [6] X. X. Zhu and R. Bamler, “Very high resolution spaceborne SAR tomography in urban environment,” _IEEE Trans. Geosci. Remote Sens._ , vol. 48, no. 12, pp. 4296–4308, Jun. 2010.
  * [7] F. Natterer, _The Mathematics of Computerized Tomography_. SIAM, 2001.
  * [8] M. Wielage, F. Cholewa, C. Fahnemann, P. Pirsch, and H. Blume, “High performance and low power architectures: GPU vs. FPGA for fast factorized backprojection,” in _Proc. Int. Symp. Comput. Netw. (CANDAR)_ , Nov. 2017, pp. 351–357.
  * [9] R. Que, O. Ponce, R. Scheiber, and A. Reigber, “Real-time processing of SAR images for linear and non-linear tracks,” in _Proc. Int. Radar Symp. (IRS)_ , May 2016, pp. 1–4.
  * [10] Y. Zeng, R. Zhang, and T. J. Lim, “Wireless communications with unmanned aerial vehicles: opportunities and challenges,” _IEEE Commun. Mag._ , vol. 54, no. 5, pp. 36–42, May 2016.
  * [11] M.-A. Lahmeri, M. A. Kishk, and M.-S. Alouini, “Artificial intelligence for UAV-enabled wireless networks: A survey,” _IEEE Open J. Commun. Soc._ , vol. 2, pp. 1015–1040, Apr. 2021.
  * [12] M. Z. Chowdhury, M. Shahjalal, S. Ahmed, and Y. M. Jang, “6G wireless communication systems: Applications, requirements, technologies, challenges, and research directions,” _IEEE Open J. Commun. Soc._ , vol. 1, pp. 957–975, 2020.
  * [13] D. Xu, Y. Sun, D. W. K. Ng, and R. Schober, “Robust resource allocation for UAV systems with UAV jittering and user location uncertainty,” in _Proc. IEEE Global Commun. Conf._ IEEE, Dec. 2018, pp. 1–6.
  * [14] W. Wang and W. Zhang, “Jittering effects analysis and beam training design for UAV millimeter wave communications,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 5, pp. 3131–3146, May 2022.
  * [15] D. Xu, Y. Sun, D. W. K. Ng, and R. Schober, “Multiuser MISO UAV communications in uncertain environments with no-fly zones: Robust trajectory and resource allocation design,” _IEEE Trans. Commun._ , vol. 68, no. 5, Jan. 2020.
  * [16] M.-A. Lahmeri, M. A. Kishk, and M.-S. Alouini, “Laser-powered UAVs for wireless communication coverage: A large-scale deployment strategy,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 1, pp. 518–533, Jan. 2022.
  * [17] L. Yang, J. Qi, J. Xiao, and X. Yong, “A literature review of UAV 3D path planning,” in _Proc. of the 11th World Congress on Intelligent Control and Automation_ , Jun. 2014, pp. 2376–2381.
  * [18] Y. Huang _et al._ , “Cognitive UAV communication via joint maneuver and power control,” _IEEE Trans. Commun._ , vol. 67, no. 11, pp. 7872–7888, 2019\.
  * [19] 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.
  * [20] H. Bayerlein, M. Theile, M. Caccamo, and D. Gesbert, “Multi-UAV path planning for wireless data harvesting with deep reinforcement learning,” _IEEE Open J. Commun. Soc._ , vol. 2, pp. 1171–1187, 2021.
  * [21] A. Hajijamali Arani _et al._ , “Reinforcement learning for energy-efficient trajectory design of UAVs,” _IEEE Internet Things Journ._ , vol. 9, no. 11, pp. 9060–9070, 2022.
  * [22] M. Lort, A. Aguasca, C. Lopez-Martinez, and T. M. Marín, “Initial evaluation of SAR capabilities in UAV multicopter platforms,” _IEEE J. Sel. Top. Appl. Earth Obs. and Remote Sens._ , vol. 11, no. 1, pp. 127–140, Jan. 2018.
  * [23] M. G. Fernández, Y. Á. López, A. A. Arboleya, B. G. Valdés, Y. R. Vaqueiro, F. L.-H. Andrés, and A. P. García, “Synthetic aperture radar imaging system for landmine detection using a ground penetrating radar on board a unmanned aerial vehicle,” _IEEE Access_ , vol. 6, pp. 45 100–45 112, 2018.
  * [24] R. Burr, M. Schartel, A. Grathwohl, W. Mayer, T. Walter, and C. Waldschmidt, “UAV-borne FMCW InSAR for focusing buried objects,” _IEEE Geosci. Remote Sens. Lett._ , vol. 19, pp. 1–5, Jul. 2021.
  * [25] Z. Sjanic and F. Gustafsson, “Simultaneous navigation and synthetic aperture radar focusing,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 51, no. 2, pp. 1253–1266, Apr. 2015.
  * [26] S. Hu, X. Yuan, W. Ni, and X. Wang, “Trajectory planning of cellular-connected UAV for communication-assisted radar sensing,” _IEEE Trans. Commun._ , vol. 70, no. 9, pp. 6385–6396, Sep. 2022.
  * [27] Z. Sun, G. G. Yen, J. Wu, H. Ren, H. An, and J. Yang, “Mission planning for energy-efficient passive UAV radar imaging system based on substage division collaborative search,” _IEEE Trans. Cybern._ , vol. 53, no. 1, pp. 275–288, Jan. 2023.
  * [28] Z. Sun, J. Wu, J. Yang, Y. Huang, C. Li, and D. Li, “Path planning for GEO-UAV bistatic SAR using constrained adaptive multiobjective differential evolution,” _IEEE Trans. Geosci. Remote Sens._ , vol. 54, no. 11, pp. 6444–6457, Nov. 2016.
  * [29] H. Cruz, M. Véstias, J. Monteiro, H. Neto, and R. P. Duarte, “A review of synthetic-aperture radar image formation algorithms and implementations: a computational perspective,” _Remote Sens._ , vol. 14, no. 5, p. 1258, Mar. 2022.
  * [30] R. Bähnemann, N. Lawrance, L. Streichenberg, J. J. Chung, M. Pantic, A. Grathwohl, C. Waldschmidt, and R. Siegwart, “Under the sand: Navigation and localization of a micro aerial vehicle for landmine detection with ground-penetrating synthetic aperture radar,” _J. Field Robot._ , vol. 2, no. 1, pp. 1028–1067, Jun. 2022.
  * [31] M. Younis, “Digital beam-forming for high resolution wide swath real and synthetic aperture radar,” Ph.D. dissertation, Zugl.: Karlsruhe, Univ., Diss., 2004.
  * [32] S. Y. Kim, N. H. Myung, and M. J. Kang, “Antenna mask design for SAR performance optimization,” _IEEE Geosci. Remote Sens. Lett._ , vol. 6, no. 3, pp. 443–447, Jul. 2009.
  * [33] A. Goldsmith, _Wireless Communications_. Cambridge University Press, 2005.
  * [34] Y. Zeng, J. Xu, and R. Zhang, “Energy minimization for wireless communication with rotary-wing UAV,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 4, pp. 2329–2345, 2019.
  * [35] M. Łabowski, P. Kaniewski, and S. Konatowski, “Estimation of flight path deviations for SAR radar installed on UAV,” _Metrology and Measurement Systems_ , vol. 23, no. 3, 2016.
  * [36] L. Zhang, Z. Qiao, M.-d. Xing, L. Yang, and Z. Bao, “A robust motion compensation approach for UAV SAR imagery,” _IEEE Trans. Geosci. Remote Sens._ , vol. 50, no. 8, pp. 3202–3218, Feb. 2012.
  * [37] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” Mar. 2014. [Online]. Available: http://cvxr.com/cvx
  * [38] Q. T. Dinh and M. Diehl, “Local convergence of sequential convex programming for nonconvex optimization,” in _Recent Advances in Optimization and its Applications in Engineering: The 14th Belgian-French-German Conf. on Optimization_. Springer, 2010, pp. 93–102.
  * [39] Y. Sun, D. Xu, D. W. K. Ng, L. Dai, and R. Schober, “Optimal 3D-trajectory design and resource allocation for solar-powered UAV communication systems,” _IEEE Trans. Commun._ , vol. 67, no. 6, pp. 4281–4298, Feb. 2019\.
  * [40] M. Hua _et al._ , “3D UAV trajectory and communication design for simultaneous uplink and downlink transmission,” _IEEE Trans. Commun._ , vol. 68, no. 9, pp. 5908–5923, 2020.
  * [41] H. Tuy, “Monotonic optimization: Problems and solution approaches,” _SIAM J. Optim._ , vol. 11, no. 2, pp. 464–494, Nov. 2000.
  * [42] Y. J. A. Zhang, L. Qian, J. Huang _et al._ , “Monotonic optimization in communication and networking systems,” _Found. Trends Netw._ , vol. 7, no. 1, pp. 1–75, Oct. 2013.
  * [43] M. Schartel, R. Burr, W. Mayer, N. Docci, and C. Waldschmidt, “UAV-based ground penetrating synthetic aperture radar,” in _Proc. IEEE MTT-S Int. Conf. Microw. Intell. Mobil._ , Apr. 2018, pp. 1–4.
  * [44] M.-A. Lahmeri _et al._ , “UAV Formation Optimization for Communication-assisted InSAR Sensing,” _accepted for presentation at IEEE Intern. Conf. Commun. (ICC) 2024, arXiv preprint arXiv:2311.06959_.