snr short = SNR, long = signal-to-noise ratio, pdf short = PDF, long =
probability density function, sar short = SAR, long = synthetic aperture
radar, insar short = InSAR, long = interferometric synthetic aperture radar,
ao short = AO, long = alternating optimization, long-plural-form = alternating
optimizations mimo short = MIMO, long = multiple-input multiple-output uav
short = UAV, long = unmanned aerial vehicle, long-plural-form = unmanned
aerial vehicles fdma short = FDMA, long = frequency-division multiple-access,
1d short = 1D, long = one-dimensional, islr short = ISLR, long = integrated
sidelobe ratio, pslr short = PSLR, long = peak sidelobe ratio, 3d short = 3D,
long = three-dimensional, pso short = PSO, long = particle swarm optimization,
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, 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 , bcd short
= BCD, long = block coordinate descent, hoa short = HoA, long = height of
ambiguity,

# UAV Formation and Resource Allocation Optimization for Communication-
Assisted 3D InSAR Sensing
††thanks: This paper is presented in part at the IEEE International Conference
on Communication, Denver, USA, Jun. 2024 [1]. This work was supported in part
by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) GRK
2680 – Project-ID 437847244.

Mohamed-Amine Lahmeri1, Víctor Mustieles-Pérez12, Martin Vossiek1, Gerhard
Krieger12, and Robert Schober1
1Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Germany
2German Aerospace Center (DLR), Microwaves and Radar Institute, Weßling,
Germany

###### Abstract

In this paper, we investigate joint unmanned aerial vehicle (UAV) formation
and resource allocation optimization for communication-assisted three-
dimensional (3D) synthetic aperture radar (SAR) sensing. We consider a system
consisting of two UAVs that perform bistatic interferometric SAR (InSAR)
sensing for generation of a digital elevation model (DEM) and transmit the
radar raw data to a ground station (GS) in real time. To account for practical
3D sensing requirements, we use non-conventional sensing performance metrics,
such as the SAR interferometric coherence, i.e., the local cross-correlation
between the two co-registered UAV SAR images, the point-to-point InSAR
relative height error, and the height of ambiguity, which together
characterize the accuracy with which the InSAR system can determine the height
of ground targets. Our objective is to jointly optimize the UAV formation,
speed, and communication power allocation for maximization of the InSAR
coverage while satisfying energy, communication, and InSAR-specific sensing
constraints. To solve the formulated non-smooth and non-convex optimization
problem, we divide it into three sub-problems and propose a novel alternating
optimization (AO) framework that is based on classical, monotonic, and
stochastic optimization techniques. The effectiveness of the proposed
algorithm is validated through extensive simulations and compared to several
benchmark schemes. Furthermore, our simulation results highlight the impact of
the UAV-GS communication link on the flying formation and sensing performance
and show that the DEM of a large area of interest can be mapped and offloaded
to ground successfully, while the ground topography can be estimated with
centimeter-scale precision. Lastly, we demonstrate that a low UAV velocity is
preferable for InSAR applications as it leads to better sensing accuracy.

##  I Introduction

The widespread use of uav has revolutionized modern technology, impacting
fields like remote sensing, communication, and disaster monitoring [2].
Specifically, uav excel in remote sensing applications using a variety of
technologies, such as cameras, lidars, and radars. Among these technologies,
radar sensing has attracted particular attention due to its capacity to
operate day and night under all weather conditions. Radar sensing can be
classified into three categories based on the dimension of the extracted
information. Firstly, 1d sensing includes applications such as ranging and
detection and necessitates at least one radar antenna. Secondly, 2d sensing,
also known as radar imaging, can be realized using, e.g., mimo radar systems
employing multiple antennas [3] or by moving a single radar antenna or a mimo
array along a trajectory to create a sar or mimo-sar system, and can provide a
high-resolution image of the illuminated area [4]. Lastly, 3d sensing provides
an additional dimension besides the spatial location of a target, such as the
displacement in time, velocity, or altitude. To this end, insar can be used,
where a given area is illuminated from different angles and at least two sar
images are formed. The image pairs are first co-registered [5], i.e., aligned
such that each ground point corresponds to the same pixel in both images, and
based on the phase difference, an interferogram is created to extract
information about the topography of the target area [5, 6].

Although trajectory and resource allocation optimization for UAV-based sensing
has been extensively studied in the literature, most works focus on 1d or 2d
sensing applications. In the context of 1d sensing, joint trajectory and
beamforming optimization for uav-based integrated sensing and communication
[7] was investigated to minimize the average power consumption [8] and to
maximize the achievable communication rate for one user [9] and for multiple
users [10]. As for 2d sensing, trajectory optimization for uav-based sar was
investigated in [11, 12, 13, 14, 15, 16, 17], where the target area was
illuminated by a geosynchronous satellite in [11, 12, 13], by a uav in [14,
16, 17], and by a communication base station in [15]. In this context,
seamless uav-to-ground connectivity enables the timely transfer of essential
information to the ground for processing. In fact, the authors in [16, 17, 13]
assumed a los link between the uav-sar system and a gs, and, in [16, 17],
real-time offloading to the gs was considered, while in [13], the authors
investigated transmitting the total radar raw data to the ground. On the other
hand, 3d sensing is more challenging due to the use of multiple sensors and
the need for considering sophisticated performance metrics. Although a
corresponding performance analysis was provided in [18, 19, 20], joint uav
formation and resource allocation optimization for 3d sensing has not been
addressed in the literature yet.

insar-based 3d sensing is fundamentally different from 1d and 2d sensing. In
fact, to assess the performance of 1d sensing, conventional metrics, such as
snr [8], beam pattern gain [10, 9], and probability of detection and false
alarm [21, 22] can be adopted. For 2d sensing, image performance metrics can
be used, including sar snr [16], spatial resolution [11], and islr and pslr
[23]. However, these performance metrics are not sufficient for 3d sensing
applications, such as generating an accurate dem based on across-track
interferometry. In fact, the key performance metric for estimating
interferometric performance is coherence, which is a function of the
correlation between the pair of co-registered sar images [6]. Additional
relevant performance metrics are the hoa, which is a proportionality constant
between the interferometric phase and the terrain height and is thus related
to the sensitivity of the radar to the ground topography [6], and the relative
height accuracy, which represents the accuracy with which the altitude of
ground objects can be estimated. In this context, an interesting trade-off
arises; while a large inter-uav separation distance, i.e., large
interferometric baseline, leads to a small hoa, which in general improves the
height accuracy of dem, it also leads to a degradation of the image coherence
[20]. Furthermore, 3d insar sensing is different from cooperative sensing [21,
24] for two main reasons; (i) the coverage is defined by the region where the
ground footprints of the two radar antennas overlap and (ii) the sensing
performance is highly dependent on the uav formation. Therefore, based on the
above discussion, existing results for uav-based 1d and 2d sensing are not
applicable to interferometric 3d sensing.

In this paper, we consider a system of two uav that perform sar sensing of a
given area from different angles, such that the difference in the sensors’
phases allows for estimating the height of ground objects based on insar
techniques. To this end, we focus on maximizing the sensing coverage by
jointly optimizing the formation, velocity, and communication transmit power
of both uav while satisfying practical energy, communication, and sensing
constraints. Our main contributions can be summarized as follows:

  * •

Different from other works, we consider insar-specific sensing performance
metrics, such as insar coverage, coherence, hoa, and relative height accuracy.

  * •

A joint formation and resource allocation optimization problem is formulated
to maximize the coverage while satisfying energy, communication, and sensing
constraints.

  * •

A novel ao-based algorithm, which incorporates classical, monotonic, and
stochastic optimization techniques, is proposed to solve the formulated
challenging non-smooth and non-convex optimization problem.

  * •

Our simulations show the superior performance of the proposed solution
compared to a conventional ao algorithm and other benchmark schemes and unveil
the impact of the sensing requirements on communication and vice versa.

We note that this article extends the corresponding conference version [1]. In
contrast to this paper, in [1], the velocity of both the master and slave uav
was fixed and the radar look angle of the slave drone was also not optimized.
Furthermore, the minimum sensing data rate was assumed to be fixed, i.e.,
independent of the imageable area. Lastly, the relative point-to-point height
accuracy was not accounted for in [1].

The remainder of this paper is organized as follows. In Section II, we present
the considered system. In Section III, the formation and resource allocation
problem is formulated, and the proposed solution is presented in Section IV.
In Section VI, the performance of the proposed scheme is assessed via
simulations, and finally, conclusions are drawn in Section VII.

Notations: In this paper, lower-case letters $x$ refer to scalar variables,
while boldface lower-case letters $\mathbf{x}$ denote vectors. $\\{a,...,b\\}$
denotes the set of all integers between $a$ and $b$ and $\emptyset$ denotes an
empty set. $|\cdot|$ denotes the absolute value operator. $\mathbb{R}^{N}$
represents the set of all $N$-dimensional vectors with real-valued entries.
For a vector $\mathbf{x}=(x_{1},...,x_{N})^{T}\in\mathbb{R}^{N}$,
$||\mathbf{x}||_{2}$ denotes the Euclidean norm, whereas $\mathbf{x}^{T}$
stands for the transpose of $\mathbf{x}$. For real-valued multivariate
functions $f(\mathbf{x})$ and $g(\mathbf{x})$, $\frac{\partial
f}{\partial\mathbf{x}}(\mathbf{a})=\Big{(}\frac{\partial f}{\partial
x_{1}}(\mathbf{a}),...,\frac{\partial f}{\partial
x_{N}}(\mathbf{a})\Big{)}^{T}$ denotes the partial derivative of $f$ wrt
$\mathbf{x}$ evaluated for an arbitrary vector $\mathbf{a}$ and
$f(x)\circledast g(x)$ refers to the convolution operation. For real numbers
$a$ and $b$, $\max(a,b)$ and $\min(a,b)$ stand for the maximum and minimum of
$a$ and $b$, respectively. For a scalar $x\in\mathbb{R}$, $[x]^{+}$ refers to
$\max(0,x)$. The notation $X_{n}(\mathbf{x}_{1},...,\mathbf{x}_{i})$
highlights that $X$ depends on optimization variables
$(\mathbf{x}_{1},...,\mathbf{x}_{i})$ and time slot $n$.

##  II System Model

Figure 1: insar sensing system with two uav sar sensors and a GS for real-time
data offloading.

We consider two rotary-wing111The results presented in this paper can be
extended to fixed-wing uav by considering a corresponding propulsion power
model. uav, denoted by $U_{1}$ and $U_{2}$, that perform insar sensing over a
given area of interest. We adopt a 3d coordinate system, where the $x$-axis
represents the ground range direction, the $y$-axis represents the azimuth
direction, and the $z$-axis defines the flying altitude of the drones above
ground, see Figure 1. We discretize the total mission time $T$ into $N$
uniform time slots of duration $\delta_{t}$ such that $T=N\cdot\delta_{t}$.
Both drones follow a linear trajectory as the stripmap sar imaging mode,
typically used in insar systems, is employed [25]. Therefore, the radar
coverage along the $x$-axis, referred to as swath, is centered wrt a line
$l_{t}$ that is parallel to the $y$-axis and passes, in time slot
$n\in\\{1,...,N\\}$, through reference point
$\mathbf{p}_{t}[n]=(x_{t},y[n],0)^{T}\in\mathbb{R}^{3}$, see Figure 1. Without
loss of generality, $U_{1}$ is the farthest uav from
$\mathbf{p}_{t}[n],\forall n,$ and is referred to as the master drone, whereas
$U_{2}$ is the slave drone. Moreover, we consider across-track interferometry
[5], where both drones transmit and receive echoes and are located in the same
$xz$-plane, also referred to as the across-track plane [5], by flying with the
same velocity $\mathbf{v}_{y}=(v_{y}[1],...,v_{y}[N])^{T}\in\mathbb{R}^{N}$,
such that in time slot $n$, the instantaneous velocity vector is given by
$\mathbf{v}[n]=(0,v_{y}[n],0)^{T}\in\mathbb{R}^{3},\forall n$, [19]. The
location of $U_{i}$ in time slot $n$, $i\in\\{1,2\\}$, is denoted by
$\mathbf{q}_{i}[n]=(x_{i},y[n],z_{i})^{T}$, where the $y$-axis position vector
$\mathbf{y}=(y[1]=0,y[2],...,y[N])^{T}\in\mathbb{R}^{N}$ is given by:

$\displaystyle y[n+1]=y[n]+v_{y}[n]\delta_{t},\forall n\in\\{1,N-1\\}.$ (1)

Hereinafter, for ease of notation, we simply use
$\mathbf{q}_{i}=(x_{i},z_{i})^{T}\in\mathbb{R}^{2},\forall i\in\\{1,2\\}$, to
denote the position of $U_{i}$ in the across-track plane (i.e., $xz-$plane),
whereas $\mathbf{q}_{i}[n]$ refers to the full 3d position of uav $U_{i}$ in
time slot $n$. Furthermore, the interferometric baseline, which refers to the
distance between the two insar sensors, is given by:

$\displaystyle
b(\mathbf{q}_{1},\mathbf{q}_{2})=||\mathbf{q}_{2}-\mathbf{q}_{1}||_{2}.$ (2)

Figure 2: Illustration of the bistatic uav formation, denoted by
$\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}$, in the across-track plane.

The perpendicular baseline, denoted by $b_{\bot}$, is the magnitude of the
projection of the baseline vector perpendicular to $U_{1}$’s los to
$\mathbf{p}_{t}[n]$ and is given by:

$b_{\bot}(\mathbf{q}_{1},\mathbf{q}_{2})=b(\mathbf{q}_{1},\mathbf{q}_{2})\cos\Big{(}\theta_{1}-\alpha(\mathbf{q}_{1},\mathbf{q}_{2})\Big{)},$
(3)

where $\theta_{1}$ is the fixed angle that $U_{1}$’s los to
$\mathbf{p}_{t}[n]$ has with the vertical, i.e., the look angle of the master
uav, and $\alpha(\mathbf{q}_{1},\mathbf{q}_{2})$ is the angle between the
interferometric baseline and the horizontal plane, see Figure 2.

###  II-A Bistatic insar Coverage

Let $r_{1}$ and $r_{2}$ denote the slant range of the radars of $U_{1}$ and
$U_{2}$ wrt $\mathbf{p}_{t}[n]$, respectively. These slant ranges are
independent of time and given by:

$\displaystyle
r_{i}(\mathbf{q}_{i})=\sqrt{(x_{i}-x_{t})^{2}+z_{i}^{2}},\forall
i\in\\{1,2\\}.$ (4)

To maximize the coverage, the master uav is positioned in the across-track
plane such that its sar antenna beam footprint on the ground is centered
around $\mathbf{p}_{t}[n]$. Furthermore, without loss of generality, we assume
$r_{2}(\mathbf{q}_{2})\leq r_{1}(\mathbf{q}_{1})$ and impose the following
condition:

$\displaystyle x_{1}=x_{t}-z_{1}\tan(\theta_{1}).$ (5)

As the look angle of the master platform is fixed to $\theta_{1}$, we propose
to adjust the look angle of the slave uav, denoted by
$\theta_{2}(\mathbf{q}_{2})$, such that its beam footprint is centered around
$\mathbf{p}_{t}[n]$, i.e.,
$\theta_{2}(\mathbf{q}_{2})=\arctan\left(\frac{x_{2}-x_{t}}{z_{2}}\right)$.
The insar coverage is limited to the area in which the beam footprints of
$U_{1}$ and $U_{2}$ overlap, see Figure 2. In this work, we focus on a single
scan (i.e., single swath) for the two drones, where the usable swath can be
obtained as follows:

$\displaystyle S(\mathbf{q}_{1},\mathbf{q}_{2})=\Big{[}\min\left(S_{\rm
far}(\mathbf{q}_{1}),S_{\rm far}(\mathbf{q}_{2})\right)-\max\left(S_{\rm
near}(\mathbf{q}_{1}),S_{\rm near}(\mathbf{q}_{2})\right)\Big{]}^{+},$ (6)

where $S_{\rm
far}(\mathbf{q}_{i})=x_{i}+z_{i}\tan(\theta_{i}(\mathbf{q}_{i})+\frac{\Theta_{\mathrm{3dB}}}{2})$,
$S_{\rm
near}(\mathbf{q}_{i})=x_{i}+z_{i}\tan(\theta_{i}(\mathbf{q}_{i})-\frac{\Theta_{\mathrm{3dB}}}{2})$,
$\forall i\in\\{1,2\\}$, $\theta_{1}(\mathbf{q}_{1})=\theta_{1}$, and
$\Theta_{\mathrm{3dB}}$ is the -3 dB beamwidth in elevation. The total area
covered by the insar radar is approximated222The approximation is due to the
elliptical shape of the beam footprint on the ground and becomes negligible
for large $N$.  as follows:

$\displaystyle
C_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})=\sum^{N-1}_{n=1}S(\mathbf{q}_{1},\mathbf{q}_{2})v_{y}[n]\delta_{t}.$
(7)

###  II-B insar Performance Metrics

Next, we introduce the relevant insar sensing performance metrics, namely the
coherence, the hoa, and the relative height accuracy of the dem.

####  II-B1 insar coherence

Two important types of decorrelation that depend on the uav formation and
affect the coherence, thereby, insar performance, are snr decorrelation and
baseline decorrelation. In particular, low snr in sar data acquisition result
in a loss of coherence between the master and slave sar images during
processing. In time slot $n$, the resulting snr decorrelation is given by
[19]:

$\displaystyle\gamma_{\mathrm{SNR},n}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})=\prod_{i\in\\{1,2\\}}\frac{1}{\sqrt{1+\mathrm{SNR}^{-1}_{i,n}(\mathbf{q}_{i},\mathbf{v}_{y})}},\forall
n,$ (8)

where $\mathrm{SNR}_{i,n}$ denotes the snr achieved by $U_{i}$ in time slot
$n$ and is given by [19]:

$\mathrm{SNR}_{i,n}(\mathbf{q}_{i},\mathbf{v}_{y})=\frac{\gamma_{r,i}}{v_{y}[n]r_{i}^{3}(\mathbf{q}_{i})\sin(\theta_{i}(\mathbf{q}_{i}))},\forall
n,$ (9)

where
$\gamma_{r,i}=\frac{\sigma_{0}P_{t,i}\;G_{t}\;G_{r}\;\lambda^{3}\;c\;\tau_{p}\;\mathrm{PRF}}{4^{4}\pi^{3}k_{b}T_{\mathrm{sys}}\;B_{\mathrm{Rg}}\;F\;L_{\mathrm{atm}}\;L_{\mathrm{sys}}\;L_{\mathrm{az}}},\forall
i\in\\{1,2\\}$. Here, $\sigma_{0}$ is the normalized backscatter coefficient,
$P_{t,i}$ is $U_{i}$’s radar transmit power, $G_{t}$ and $G_{r}$ are the
transmit and receive antenna gains, respectively, $\lambda$ is the wavelength,
$c$ is the speed of light, $\tau_{p}$ is the pulse duration, $\mathrm{PRF}$ is
the pulse repetition frequency, $k_{b}$ is the Boltzmann constant,
$T_{\mathrm{sys}}$ is the receiver temperature, $B_{\mathrm{Rg}}$ is the
bandwidth of the radar pulse, $F$ is the noise figure, and $L_{\mathrm{atm}}$,
$L_{\mathrm{sys}}$, and $L_{\mathrm{az}}$ represent the atmospheric, system,
and azimuth losses, respectively. Another relevant type of decorrelation is
the baseline decorrelation. It reflects the loss of coherence caused by the
different angles used for acquisition of the two sar images. The baseline
decorrelation is given by[26]:

$\gamma_{\mathrm{Rg}}(\mathbf{q}_{2})=\frac{(2+B_{p})\sin\left(\theta_{2}(\mathbf{q}_{2})\right)-(2-B_{p})\sin(\theta_{1})}{B_{p}\big{(}\sin(\theta_{1})+\sin\left(\theta_{2}(\mathbf{q}_{2})\right)\big{)}},$
(10)

where $B_{p}=\frac{B_{\mathrm{Rg}}}{f_{0}}$ is the fractional bandwidth and
$f_{0}$ is the radar center frequency. It can be shown that large uav
formations, i.e., large interferometric baselines, result in high baseline
decorrelation and, therefore, degrade the coherence. The total coherence is
given by $\gamma=\gamma_{\rm Rg}\gamma_{\mathrm{SNR},n}\gamma_{\rm other}$
[6], where $\gamma_{\rm other}$ represents the contribution from all other
decorrelation sources, e.g., volumetric decorrelation. Here, we assume that
surface scattering is dominant such that $\gamma_{\rm other}$ is a constant
[6].

####  II-B2 Height of Ambiguity (HoA)

The hoa is defined as the height difference that results in a complete $2\pi$
cycle of the interferometric phase [19]. It is therefore related to the
sensitivity of the radar system to topographic height variations. Similar to
the baseline decorrelation, the hoa depends on the uav formation and is given
by [19]:

$\displaystyle h_{\mathrm{amb}}(\mathbf{q}_{1},\mathbf{q}_{2})=\frac{\lambda
r_{1}(\mathbf{q}_{1})\sin(\theta_{1})}{b_{\perp}(\mathbf{q}_{1},\mathbf{q}_{2})}.$
(11)

Note that a large uav baseline leads to a small hoa value, which results in
increased errors for interferometric phase unwrapping [27]. However, a very
small uav baseline, i.e., a large hoa, is not a good choice either, as we
explain in the following.

####  II-B3 Sensing Height Accuracy

A common statistics-based performance metric for assessing the accuracy of the
height estimation in dem is the point-to-point 90% relative height error given
by [19]:

$\Delta
h_{90\%}(\mathbf{q}_{1},\mathbf{q}_{2},\gamma)=h_{\mathrm{amb}}(\mathbf{q}_{1},\mathbf{q}_{2})\frac{\Delta\phi_{90\%}(\gamma)}{2\pi},$
(12)

where $\gamma$ is the total interferometric coherence and $\Delta\phi_{90\%}$
is defined as the 90% percentile of the random phase error $\Delta\phi$ and
obtained by evaluating the difference between two random variables, each of
them describing the fluctuation of the interferometric phase within one
interferogram as follows [19]:

$\int_{-\Delta\phi_{90\%}(\gamma)}^{{\Delta\phi_{90\%}(\gamma)}}p_{\Delta\phi}(\Delta\phi,\gamma)\circledast
p_{\Delta\phi}(\Delta\phi,\gamma)\cdot d\Delta\phi=0.9,$ (13)

where $p_{\Delta\phi}$ is the pdf of the phase error $\Delta\phi$ for a given
interferometric coherence $\gamma$ and is given by333In [28], the pdf is
derived by assuming Gaussian scattering-matrix statistics. Yet, the proposed
solution can be extended to other models, e.g., the K-distribution model [29].
[28]:

$p_{\Delta\phi}(\Delta\phi,\gamma)=\frac{\Gamma(n_{L}+0.5)(1-|\gamma|^{2})^{n_{L}}|\gamma|\cos(\Delta\phi)}{2\sqrt{\pi}\Gamma(n_{L})\left(1-|\gamma|^{2}\cos^{2}(\Delta\phi)\right)^{n_{L}+0.5}}+\frac{(1-|\gamma|^{2})^{n_{L}}}{2\pi}F_{G}\left(n_{L},1;0.5;|\gamma|^{2}\cos^{2}(\Delta\phi)\right).$
(14)

Here, $\Gamma$ is the Gamma function, $F_{G}$ is the Gauss hypergeometric
function [19], and $n_{L}$ is the number of independent looks. In fact, multi-
looking is a technique that averages adjacent pixels of the interferogram,
e.g., by filtering with a boxcar window, to improve phase estimation leading
to a trade-off between spatial resolution and height accuracy, see [28]. Note
that a large hoa results in a large point-to-point 90% height error $\Delta
h_{90\%}$, and thus, in a low relative height estimation accuracy. This
explains why optimizing the uav baseline, and therefore, the hoa, is an
intricate problem.

###  II-C Communication Performance

We consider real-time offloading of the radar data to a gs, where the master
and slave uav employ fdma. The instantaneous communication transmit power
consumed by uav $U_{i}$ is given by
$\mathbf{P}_{\mathrm{com},i}=(P_{\mathrm{com},i}[1],...,P_{\mathrm{com},i}[N])^{T}\in\mathbb{R}^{N},i\in\\{1,2\\}$.
We denote the location of the gs by
$\mathbf{g}=(g_{x},g_{y},g_{z})^{T}\in\mathbb{R}^{3}$ and the distance from
$U_{i}$ to the gs by
$d_{i,n}(\mathbf{q}_{i},\mathbf{v}_{y})=||\mathbf{q}_{i}[n]-\mathbf{g}||_{2},\forall
i\in\\{1,2\\},\forall n.$ We suppose that both uav fly at sufficiently high
altitudes to allow obstacle-free communication with the gs over a los link.
Thus, based on the free-space path loss model and fdma, the instantaneous
throughput from $U_{i},\forall i\in\\{1,2\\},$ to the gs is given by:

$\displaystyle
R_{i,n}(\mathbf{q}_{i},\mathbf{P}_{\mathrm{com},i},\mathbf{v}_{y})=B_{c,i}\;\log_{2}\left(1+\frac{P_{\mathrm{com},i}[n]\;\beta_{c,i}}{d_{i,n}^{2}(\mathbf{q}_{i},\mathbf{v}_{y})}\right),\forall
n,$ (15)

where $B_{c,i}$ is the fixed communication bandwidth allocated for $U_{i}$ and
$\beta_{c,i}$ is its reference channel gain444The reference channel gain is
the channel power gain at a reference distance of 1 m. divided by the noise
variance. The minimum required data rate, usually fixed in communication
systems, is significantly impacted by the formation of the two UAVs and
depends on the imageable area size. Based on [30], the sar sensing data rate
for $U_{i},\forall i\in\\{1,2\\},$ is given by:

$R_{\mathrm{min},i}(\mathbf{q}_{i})=n_{B}B_{\mathrm{Rg}}\mathrm{PRF}\left(\frac{z_{i}}{c}\left(\frac{1}{\cos\left(\theta_{i}(\mathbf{q}_{i})+\frac{\Theta_{\mathrm{3dB}}}{2}\right)}-\frac{1}{\cos\left(\theta_{i}(\mathbf{q}_{i})-\frac{\Theta_{\mathrm{3dB}}}{2}\right)}\right)+\tau_{p}\right),$
(16)

where $n_{B}$ is the number of bits per complex sample.

###  II-D Energy Consumption

The propulsion power $P_{\mathrm{prop,n}}$ consumed by the drone in time slot
$n$ is given by [31]:

$P_{\mathrm{prop},n}(\mathbf{v}_{y})=P_{0}\left(1+\frac{3v_{y}^{2}[n]}{U^{2}_{\mathrm{tip}}}\right)+P_{I}\left(\sqrt{1+\frac{v^{4}_{y}[n]}{4v_{0}^{4}}}-\frac{v^{2}_{y}[n]}{2v_{0}^{2}}\right)^{\frac{1}{2}}+\frac{1}{2}d_{0}\rho
sA_{e}v^{3}_{y}[n],\forall n.$ (17)

Here, $P_{0}=\frac{\delta_{u}}{8}\rho sA_{e}\Omega^{3}R^{3}$ and
$P_{I}=(1+k_{u})\frac{W_{u}^{\frac{3}{2}}}{\sqrt{2\rho A_{e}}}$ are two
constants, $v_{0}=\sqrt{\frac{W_{u}}{2\rho A_{e}}}$ is the mean rotor induced
velocity while hovering, $U_{\mathrm{tip}}$ is the tip speed of the rotor
blade, $d_{0}$ is the fuselage drag ratio, $\delta_{u}$ is the profile drag
coefficient, $\rho$ is the air density, $s$ is the rotor solidity, $A_{e}$ is
the rotor disc area, $\Omega$ is the blade angular velocity, $R$ is the rotor
radius, $k_{u}$ is a correction factor, and $W_{u}$ is the aircraft weight in
Newton. The total energy consumed by uav $U_{i},\forall i\in\\{1,2\\},$ is
given by:

$E_{i,N}(\mathbf{P}_{\mathrm{com},i},\mathbf{v}_{y})=\sum\limits_{n=1}^{N}\delta_{t}\Big{(}P_{\mathrm{prop},n}(\mathbf{v}_{y})+P_{t,i}+P_{\mathrm{com},i}[n]\Big{)}.$
(18)

##  III Problem Formulation

In this paper, we aim to maximize the insar coverage by jointly optimizing the
uav formation $\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}$, the instantaneous
communication transmit powers
$\\{\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,2}}\\}$, and the
instantaneous uav velocity $\mathbf{v}_{y}$, while satisfying energy,
communication, and interferometric quality-of-service constraints. To this
end, we formulate the following optimization problem:

$\displaystyle(\mathrm{P.1}):\max_{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{P}_{\mathrm{com},1},\mathbf{P}_{\mathrm{com},2},\mathbf{v}_{y}}\hskip
8.53581ptC_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})$ s.t.
$\displaystyle\mathrm{C1:}\;z_{\mathrm{min}}\leq z_{i}\leq
z_{\mathrm{max}},\forall i\in\\{1,2\\},$
$\displaystyle\mathrm{C2}:\;x_{1}=x_{t}-z_{1}\tan(\theta_{1}),$
$\displaystyle\mathrm{C3}:\;r_{2}(\mathbf{q}_{2})\leq r_{1}(\mathbf{q}_{1}),$
$\displaystyle\mathrm{C4}:\;x_{2}\leq x_{t},$
$\displaystyle\mathrm{C5}:b(\mathbf{q}_{1},\mathbf{q}_{2})\geq
b_{\mathrm{min}},$
$\displaystyle\mathrm{C6}:\gamma_{\mathrm{SNR},n}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})\geq\gamma_{\mathrm{SNR}}^{\mathrm{min}},\forall
n,$
$\displaystyle\mathrm{C7}:\gamma_{\mathrm{Rg}}(\mathbf{q}_{2})\geq\gamma_{\mathrm{Rg}}^{\mathrm{min}},$
$\displaystyle\mathrm{C8}:\;h_{\mathrm{amb}}(\mathbf{q}_{1},\mathbf{q}_{2})\geq
h_{\mathrm{amb}}^{\mathrm{min}},$ $\displaystyle\mathrm{C9}:\;\Delta
h_{90\%}^{\rm worst}(\mathbf{q}_{1},\mathbf{q}_{2})\leq\Delta
h^{\mathrm{max}},$ $\displaystyle\mathrm{C10}:0\leq P_{\mathrm{com},i}[n]\leq
P_{\mathrm{com}}^{\mathrm{max}},\forall\;i\in\\{1,2\\},\forall n,$
$\displaystyle\mathrm{C11}:R_{i,n}(\mathbf{q}_{i},\mathbf{P}_{\mathrm{com},i},\mathbf{v}_{y})\geq
R_{\mathrm{min},i}(\mathbf{q}_{i}),\forall\;i\in\\{1,2\\},\forall n,$
$\displaystyle\mathrm{C12}:E_{i,N}(\mathbf{P}_{\mathrm{com},i},\mathbf{v}_{y})\leq
E_{\mathrm{max},i},\forall\;i\in\\{1,2\\},$
$\displaystyle\mathrm{C13}:v_{\mathrm{min}}\leq v_{y}[n]\leq
v_{\mathrm{max}},\forall n,$
$\displaystyle\mathrm{C14}:\theta_{\mathrm{min}}\leq\theta_{2}(\mathbf{q}_{2})\leq\theta_{\mathrm{max}},$
$\displaystyle\mathrm{C15}:\theta_{2}(\mathbf{q}_{2})=\arctan\left(\frac{x_{2}-x_{t}}{z_{2}}\right).$

Constraint $\mathrm{C1}$ specifies the maximum and minimum allowed flying
altitude, denoted by $z_{\mathrm{max}}$ and $z_{\mathrm{min}}$, respectively.
Constraints $\mathrm{C2}$ and $\mathrm{C3}$ ensure maximum overlap between the
beam footprint of the master drone and the area of interest. Constraint
$\mathrm{C4}$ is imposed because a side-looking sar is assumed. Constraint
$\mathrm{C5}$ ensures a minimum separation distance between the uav, denoted
by $b_{\mathrm{min}}$. Constraints $\mathrm{C6}$ and $\mathrm{C7}$ ensure
minimum required thresholds for the sensing $\mathrm{SNR}$ and baseline
decorrelation, denoted by $\gamma_{\mathrm{SNR}}^{\mathrm{min}}$ and
$\gamma_{\mathrm{Rg}}^{\mathrm{min}}$, respectively. Constraint $\mathrm{C8}$
imposes a minimum hoa, denoted by $h_{\mathrm{amb}}^{\mathrm{min}}$, that
satisfies phase unwrapping requirements [6]. Constraint $\mathrm{C9}$ ensures
that the worst-case point-to-point 90% relative height error of uav formation
$\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}$, defined as $\Delta h_{90\%}^{\rm
worst}(\mathbf{q}_{1},\mathbf{q}_{2})=\Delta
h_{90\%}(\mathbf{q}_{1},\mathbf{q}_{2},\gamma_{\rm SNR}^{\rm min}\gamma_{\rm
Rg}^{\rm min}\gamma_{\rm other})$, does not exceed the maximum allowed height
error, denoted by $\Delta h^{\rm max}$. The definition of $\Delta
h_{90\%}^{\rm worst}$ is based on the fact that (i) an analytical expression
for (13) cannot be derived and (ii) the phase error variance decreases with
the coherence [28]. Therefore, (13) is solved numerically for the worst case
scenario $\gamma=\gamma_{\rm SNR}^{\rm min}\gamma_{\rm Rg}^{\rm
min}\gamma_{\rm other}$. Constraint $\mathrm{C10}$ ensures that the
communication transmit power is non-negative and does not exceed the maximum
allowed level denoted by $P_{\mathrm{com}}^{\mathrm{max}}$. Constraint
$\mathrm{C11}$ ensures that the achievable throughput of drone $U_{i}$ does
not fall below the minimum required data rate
$R_{\mathrm{min},i}(\mathbf{q}_{i}),\forall i\in\\{1,2\\}$, which depends on
the imageable area, thus the uav formation. Constraint $\mathrm{C12}$ limits
the total consumed energy of uav $U_{i}$ to its maximum battery capacity,
denoted by $E_{\mathrm{max},i},\forall i\in\\{1,2\\}$. Constraint
$\mathrm{C13}$ specifies the minimum and maximum allowed velocities of the
uav-sar system, denoted by $v_{\rm min}$ and $v_{\rm max}$, respectively.
Constraint $\mathrm{C14}$ limits the maximum and minimum slave look angles to
$\theta_{\rm max}$ and $\theta_{\rm min}$, respectively. Constraint
$\mathrm{C15}$ directs the slave look angle toward $\mathbf{p}_{t}[n],\forall
n$. We note that some of the constraints do not depend on time $n$ as some
variables, such as $x_{1}$ and $z_{2}$, are optimized but are fixed across all
time slots due to the prescribed linear insar trajectory imposed by the
stripmap mode sar operation [25].

Problem $\mathrm{(P.1)}$ is a non-smooth non-convex optimization problem for
the following reasons. First, the use of a variable slave look angle as well
as $\max$ and $\min$ functions causes the objective function $C_{N}$ to be
non-smooth and non-convex wrt $\mathbf{q}_{2}$. Second, the coupling between
different optimization variables and the presence of trigonometric functions
in the expressions for the $\mathrm{SNR}$ decorrelation, baseline
decorrelation, and hoa (see (8), (10), and (11)), make constraints
$\mathrm{C6}$, $\mathrm{C7}$, $\mathrm{C8}$, and $\mathrm{C9}$ non-convex.
Third, the use of variable minimum sensing data rates causes $\mathrm{C11}$ to
be neither monotonic nor convex. Fourth, constraints $\mathrm{C3}$ and
$\mathrm{C5}$ involve lower bounds on Euclidean distances and, thus, are also
non-convex. Lastly, the uav propulsion power, given in (17), causes constraint
$\mathrm{C12}$ to be non-convex. Overall, it is very difficult and challenging
to find a globally optimal solution to problem $\mathrm{(P.1)}$. To overcome
this challenge, in the next section, we divide problem $\mathrm{(P.1)}$ into
several more tractable and easier-to-solve sub-problems and tackle them in an
alternating manner.

##  IV Solution of the Optimization Problem

Figure 3:  Block diagram of the proposed solution to problem $\rm(P.1)$ based
on ao, sca, polyblock outer approximation, and pso.

To strike a balance between performance and complexity, we provide a low-
complexity sub-optimal solution for the formulated problem based on ao. To
this end, in each iteration $m$ of the proposed ao algorithm, problem
$\mathrm{(P.1)}$ is divided into three sub-problems, namely
$\mathrm{(P.1.a)}$, $\mathrm{(P.1.b)}$, and $\mathrm{(P.1.c)}$, that are
solved alternately, as illustrated in Figure 3. In problems $\mathrm{(P.1.a)}$
and $\mathrm{(P.1.b)}$, we optimize the positions $\mathbf{q}_{2}$ and
$\mathbf{q}_{1}$ of the slave and master uav in the across-track plane,
respectively. In sub-problem $\mathrm{(P.1.c)}$, we optimize the uav
resources, namely the instantaneous velocity $\mathbf{v}_{y}$ (which affects
the propulsion power $\mathbf{P}_{\rm prop}$) and the instantaneous
communication transmit powers, $\mathbf{P}_{\rm com,1}$ and $\mathbf{P}_{\rm
com,2}$. Each sub-problem has different characteristics, therefore, in the
following subsections, we provide the solution for each sub-problem using
different optimization techniques.

###  IV-A Slave uav Optimization

In this subsection, we aim to optimize the across-track position of the slave
uav, $\mathbf{q}_{2}$, for given resources
$\\{\mathbf{v}_{y},\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,2}}\\}$
and given master uav position $\mathbf{q}_{1}$. Therefore, problem
$\mathrm{(P.1)}$ reduces to sub-problem $\mathrm{(P.1.a)}$ given by:

$\displaystyle(\mathrm{P.1.a}):\max_{\mathbf{q}_{2}}\hskip
8.53581ptC_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})$ s.t.
$\displaystyle\mathrm{C3-C9,C14,C15},$
$\displaystyle\mathrm{C1:}\;z_{\mathrm{min}}\leq z_{2}\leq z_{\mathrm{max}},$
$\displaystyle\mathrm{C11}:R_{2,n}(\mathbf{q}_{2},\mathbf{P}_{\mathrm{com},2},\mathbf{v}_{y})\geq
R_{\mathrm{min},2}(\mathbf{q}_{2}),\forall n.$

Sub-problem $(\mathrm{P.1.a})$ is a non-smooth non-convex optimization problem
due to its objective function and constraints $\mathrm{C5}$, $\mathrm{C7-C9}$,
and $\mathrm{C11}$. The presence of trigonometric functions in constraints
$\mathrm{C7}$, $\mathrm{C8}$, and $\mathrm{C9}$, as well as the non-smoothness
of the objective function wrt $\mathbf{q}_{2}$ make it very difficult to solve
problem $(\mathrm{P.1.a})$ using classical optimization methods that are
generally based on gradient descent. Nevertheless, we provide a low-complexity
close-to-optimal solution based on stochastic optimization [32]. Motivated by
its success in solving complex non-smooth optimization problems and its fast
convergence [33], we use particle swarm optimization (PSO), which is a
population-based stochastic optimization framework that mimics the intelligent
collective behavior of animals, such as birds [34]. To solve sub-problem
$\mathrm{(P.1.a)}$, a population of $D$ particles is employed to explore the
2d across-track plane. In iteration $k$ of the pso algorithm, particle $d$,
denoted by $\mathbf{p}_{d}^{(k)}\in\mathbb{R}^{2},\forall d,$ has two
dimensions corresponding to the $x$\- and $z$-coordinates of the slave uav as
follows:

$\mathbf{p}_{d}^{(k)}=(x_{2,d}^{(k)},z_{2,d}^{(k)})^{T},\forall d,$ (19)

where $x_{2,d}^{(k)}$ and $z_{2,d}^{(k)}$ correspond to the $x$\- and
$z$-positions of particle $d$ in iteration $k$, respectively. Furthermore, the
velocity vector of particle $d$ in iteration $k$ is denoted by
$\mathbf{v}_{\mathrm{PSO},d}^{(k)}\in\mathbb{R}^{2}$ and given by:

$\mathbf{v}_{\mathrm{PSO},d}^{(k)}=\big{(}v_{\mathrm{PSO},d}^{(k)}[1],v_{\mathrm{PSO},d}^{(k)}[2]\big{)}^{T},\forall
d,$ (20)

where $v_{\mathrm{PSO},d}^{(k)}[1]$ and $v_{\mathrm{PSO},d}^{(k)}[2]$ are the
$x$\- and $z$-components of the particle velocity vector, respectively. In
iteration $k\geq 2$, the velocity of each particle is affected by its previous
velocity, given by $\mathbf{v}_{\mathrm{PSO},d}^{(k-1)}$, its previous local
experience, given by its previous local best-known position, denoted by
$\mathbf{p}_{\mathrm{best},d}^{(k-1)}\in\mathbb{R}^{2}$, and its global
experience, given by the best-known position, denoted by
$\mathbf{p}_{\mathrm{best}}^{(k-1)}\in\mathbb{R}^{2}$, as shown below [35]:

$\mathbf{v}_{\mathrm{PSO},d}^{(k)}=\stackrel{{\scriptstyle\text{Inertial
term}}}{{\overbrace{w(k-1)\mathbf{v}_{\mathrm{PSO},d}^{(k-1)}}}}+\stackrel{{\scriptstyle\text{Cognitive
term}}}{{\overbrace{c_{1}\mathcal{R}_{1}(\mathbf{p}_{\mathrm{best},d}^{(k-1)}-\mathbf{p}_{d}^{(k-1)})}}}+\stackrel{{\scriptstyle\text{Social
term}}}{{\overbrace{c_{2}\mathcal{R}_{2}(\mathbf{p}_{\mathrm{best}}^{(k-1)}-\mathbf{p}_{d}^{(k-1)})}}},\forall
k\geq 2,$ (21)

where $c_{1}$ and $c_{2}$ denote the cognitive and social learning factors,
respectively, $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ are random variables
that are uniformly distributed in $[0,1]$, and $w(k-1)\in[0,1]$ is the
inertial weight in iteration $k-1$. Different methods to adjust the inertial
weight were proposed in the literature [36]. In this work, we employ a
linearly decaying inertial weight. Moreover, the initial velocity vectors,
denoted by $\mathbf{v}_{\mathrm{PSO},d}^{(1)},\forall d$, are uniformly
distributed in $[0,v_{\rm PSO}^{\rm max}]\times[0,v_{\rm PSO}^{\rm max}]$,
where $v_{\rm PSO}^{\rm max}$ is the maximum particle velocity. Note that the
updates in (21) might cause particles to move outside of the feasible search
space. To prevent this from happening and to speed up the convergence of the
pso algorithm, we define boundaries for the search space based on constraints
$\mathrm{C1}$ and $\mathrm{C4}$ by defining a reflecting wall [37] as follows:

$\begin{dcases}v_{\mathrm{PSO},d}^{(k)}[1]&=-v_{\mathrm{PSO},d}^{(k)}[1]\hskip
14.22636pt\text{ if }x_{2,d}^{(k)}+v_{\mathrm{PSO},d}^{(k)}[1]>x_{t},\\\
v_{\mathrm{PSO},d}^{(k)}[2]&=-v_{\mathrm{PSO},d}^{(k)}[2]\hskip
14.22636pt\text{ if
}z_{2,d}^{(k)}+v_{\mathrm{PSO},d}^{(k)}[2]>z_{\mathrm{max}}\text{ or
}z_{2,d}^{(k)}+v_{\mathrm{PSO},d}^{(k)}[2]<z_{\mathrm{min}},\end{dcases},\forall
d,\forall k\geq 2.$ (22)

Next, in iteration $k$, the particle position is updated as follows:

$\mathbf{p}_{d}^{(k+1)}=\mathbf{p}_{d}^{(k)}+\mathbf{v}_{\mathrm{PSO},d}^{(k)},\forall
d,\forall k\geq 1.$ (23)

where initial particles $\mathbf{p}_{d}^{(1)},\forall d,$ are uniformly
distributed in $[x_{t}-O,x_{t}]\times[z_{\rm min},z_{\rm max}]$ according to
$\mathrm{C1}$ and $\mathrm{C4}$, where $O$ is set large enough, such that the
two drones are placed sufficiently behind $x_{t}$, see constraint
$\mathrm{C4}$. The remaining constraints of sub-problem $(\mathrm{P.1.a})$ as
well as the objective function are used to construct a fitness function
necessary to evaluate a given particle of the population. While the pso
algorithm was originally proposed to solve unconstrained optimization problems
[34], several variations have been proposed to accommodate constraints, e.g.,
[38, 35, 39]. In this work, we adopt a non-parameterized pso algorithm, which
was proposed in [39]. In particular, our fitness function $\mathcal{F}$ is
given by:

$\mathcal{F}(\mathbf{p}_{d}^{(k)})=\begin{dcases}C_{N}(\mathbf{q}_{1},\mathbf{p}_{d}^{(k)},\mathbf{v}_{y}),\hskip
14.22636pt&\text{ if }\mathbf{p}_{d}^{(k)}\in S,\\\
\min\limits_{k,d}\left(C_{N}(\mathbf{q}_{1},\mathbf{p}_{d}^{(k)},\mathbf{v}_{y})\right)+\sum\limits_{l\in\mathcal{L}}g_{l}(\mathbf{p}_{d}^{(k)}),\hskip
14.22636pt&\text{ if }\mathbf{p}_{d}^{(k)}\notin S,\\\ \end{dcases}\hskip
14.22636pt\forall d,$ (24)

where $S$ is the feasible set for sub-problem $\mathrm{(P.1.a)}$, which is
defined by constraints $\mathrm{C3}$, $\mathrm{C5-C9}$, $\mathrm{C11}$,
$\mathrm{C14}$, and $\mathrm{C15}$, and functions
$g_{l},l\in\mathcal{L}=\\{3,5,6,7,8,9,11,14\\}$, which are used as metrics to
quantify constraint violations as follows:

$\displaystyle
g_{3}(\mathbf{p}_{d}^{(k)})=\left[r_{2}(\mathbf{p}_{d}^{(k)})-r_{1}(\mathbf{q}_{1})\right]^{+},\forall
d,$ (25) $\displaystyle
g_{5}(\mathbf{p}_{d}^{(k)})=\left[b_{\mathrm{min}}-b(\mathbf{q}_{1},\mathbf{p}_{d}^{(k)})\right]^{+},\forall
d,$ (26) $\displaystyle
g_{6}(\mathbf{p}_{d}^{(k)})=\left[\gamma_{\mathrm{SNR}}^{\mathrm{min}}-\sum\limits_{n=1}^{N}\gamma_{\mathrm{SNR},n}(\mathbf{q}_{1},\mathbf{p}_{d}^{(k)},\mathbf{v}_{y})\right]^{+},\forall
d,$ (27) $\displaystyle
g_{7}(\mathbf{p}_{d}^{(k)})=\left[\gamma_{\mathrm{Rg}}^{\mathrm{min}}-\gamma_{\mathrm{Rg}}(\mathbf{p}_{d}^{(k)})\right]^{+},\forall
d,$ (28) $\displaystyle
g_{8}(\mathbf{p}_{d}^{(k)})=\left[h_{\mathrm{amb}}^{\mathrm{min}}-h_{\mathrm{amb}}(\mathbf{q}_{1},\mathbf{p}_{d}^{(k)})\right]^{+},\forall
d,$ (29) $\displaystyle g_{9}(\mathbf{p}_{d}^{(k)})=\left[\Delta h_{90\%}^{\rm
worst}(\mathbf{q}_{1},\mathbf{p}_{d}^{(k)})-\Delta
h^{\mathrm{max}}\right]^{+},\forall d,$ (30) $\displaystyle
g_{11}(\mathbf{p}_{d}^{(k)})=\left[R_{\mathrm{min},2}(\mathbf{p}_{d}^{(k)})-\sum\limits_{n=1}^{N}R_{2,n}(\mathbf{p}_{d}^{(k)},\mathbf{P}_{\mathrm{com},2},\mathbf{v}_{y})\right]^{+},\forall
d,$ (31) $\displaystyle g_{14}(\mathbf{p}_{d}^{(k)})=\left[\left(\theta_{\rm
min}-\theta_{2}(\mathbf{p}_{d}^{(k)})\right)+\left(\theta_{2}(\mathbf{p}_{d}^{(k)})-\theta_{\rm
max}\right)\right]^{+},\forall d.$ (32)

Algorithm 1 Particle Swarm Optimization Algorithm

1:For given
$\\{\mathbf{q}_{1},\mathbf{v}_{y},\mathbf{P}_{\mathrm{com},1},\mathbf{P}_{\mathrm{com},2}\\}$,
set initial iteration number $k=1$, population of $D$ particles
$(\mathbf{p}_{1}^{(1)},...,\mathbf{p}_{D}^{(1)})$ with initial random
velocities
$(\mathbf{v}_{\mathrm{PSO},1}^{(1)},...,\mathbf{v}_{\mathrm{PSO},D}^{(1)})$,
initial best local position for each particle
$(\mathbf{p}_{\mathrm{best},1}^{(1)}=\mathbf{p}_{1}^{(1)},...,\mathbf{p}_{\mathrm{best},D}^{(1)}=\mathbf{p}_{D}^{(1)})$,
initial global best swarm position
$\mathbf{p}_{\mathrm{best}}^{(1)}=\arg\max\limits_{\mathbf{p}_{1}^{(1)},...,\mathbf{p}_{D}^{(1)}}(\mathcal{F}(\mathbf{p}_{1}^{(1)}),...,\mathcal{F}(\mathbf{p}_{D}^{(1)}))$,
and initial inertial weight $w(1)=1$.

2:repeat

3: Adjust velocities $\mathbf{v}_{\mathrm{PSO},d}^{(k)},\forall d,$ based on
(22).

4: Update positions
$\mathbf{p}_{d}^{(k+1)}=\mathbf{p}_{d}^{(k)}+\mathbf{v}_{\mathrm{PSO},d}^{(k)},\forall
d$

5: Update population velocities
$\mathbf{v}_{\mathrm{PSO},d}^{(k+1)}=w(k)\mathbf{v}_{\mathrm{PSO},d}^{(k)}+$

6:
$c_{1}\mathcal{R}_{1}(\mathbf{p}_{\mathrm{best},d}^{(k)}-\mathbf{p}_{d}^{(k)})+c_{2}\mathcal{R}_{2}(\mathbf{p}_{\mathrm{best}}^{(k)}-\mathbf{p}_{d}^{(k)}),\forall
d$

7: Update inertial weight $w(k+1)=1-\frac{k}{M_{1}}$

8: Update local best known position $\mathbf{p}_{\mathrm{best},d}^{(k+1)}$ for
each

9: particle and the swarm’s best known position

10: $\mathbf{p}_{\mathrm{best}}^{(k+1)}$ based on fitness function
$\mathcal{F}$ in (24)

11: Set iteration number $k=k+1$

12:until convergence

13:return solution $\mathbf{q}_{2}=\mathbf{p}_{\rm best}^{(k-1)}$

The proposed procedure to solve sub-problem $\mathrm{(P.1.a)}$ is summarized
in Algorithm 1. According to our simulations, Algorithm 1 converges to a
close-to-optimal solution, which can be assessed using an upper bound for
$\mathrm{(P.1.a)}$;
$C_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})\leq\big{(}S_{\rm
far}(\mathbf{q}_{1})-S_{\rm
near}(\mathbf{q}_{1})\big{)}\sum\limits_{n=1}^{N-1}v_{y}[n]\delta_{t}$. Next,
we analyze the worst-case computational time complexity of Algorithm 1\. To
this end, let $M_{1}$ be the maximum number of iterations for Algorithm 1.
Then, to evaluate the fitness function for each particle, $2N+6$ computations
are required, see (25)-(32). Therefore, the initialization of Algorithm 1
requires $1+3D+D(2N+6)$ computations. Consequently, the overall worst-case
time complexity of Algorithm 1 is given by
$\mathcal{O}\big{(}1+3D+D(2N+6)+M_{1}D(4+2N+6)\big{)}$, which is equivalent to
$\mathcal{O}(M_{1}DN)$.

###  IV-B Master uav Position Optimization

Next, we focus on optimizing the master uav position in the across-track
plane, $\mathbf{q}_{1}$, for given
$\\{\mathbf{v}_{y},\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,2}},\mathbf{q}_{2}\\}$.
In this case, problem $\mathrm{(P.1)}$ reduces to sub-problem
$\mathrm{(P.1.b)}$ given by:

$\displaystyle(\mathrm{P.1.b}):\max_{\mathbf{q}_{1}}\hskip
8.53581ptC_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})$ s.t.
$\displaystyle\mathrm{C2},\mathrm{C3},\mathrm{C5},\mathrm{C6},\mathrm{C8},\mathrm{C9},$
$\displaystyle\mathrm{C1}:\;z_{\mathrm{min}}\leq z_{1}\leq z_{\mathrm{max}},$
$\displaystyle\mathrm{C11}:R_{1,n}(\mathbf{q}_{1},\mathbf{P}_{\mathrm{com},1},\mathbf{v}_{y})\geq
R_{\mathrm{min},1}(\mathbf{q}_{1}),\forall n.$

Sub-problem $(\mathrm{P.1.b})$ is non-convex due to its objective function and
constraints $\mathrm{C5}$ and $\mathrm{C11}$. Thus, finding the global optimal
solution to sub-problem $(\mathrm{P.1.b})$ is challenging. Yet, we provide the
optimal solution to sub-problem $(\mathrm{P.1.b})$ based on monotonic
optimization theory. To this end, we start by transforming $(\mathrm{P.1.b})$
into the canonical form of a monotonic optimization problem [40].

######  Proposition 1.

Sub-problem $(\mathrm{P.1.b})$ is equivalent to the following monotonic
optimization problem:

$\displaystyle(\overline{\mathrm{P.1.b}}):\max_{z_{1},t}\hskip
8.53581ptC_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})$ s.t.
$\displaystyle\mathrm{C6},$
$\displaystyle\mathrm{\overline{C1a}}:\;z_{1}-Z_{\mathrm{max}}\leq 0,$
$\displaystyle\mathrm{\overline{C1b}}:\;z_{1}-Z_{\mathrm{min}}\geq 0,$
$\displaystyle\overline{\mathrm{C11a}}:{a_{1}(z_{1})}+t\leq{a_{1}(z_{\rm
max})},$
$\displaystyle\overline{\mathrm{C11b}}:{a_{2}(z_{1})}+t\geq{a_{1}(z_{\rm
max})},$ $\displaystyle\overline{\mathrm{C11c}}:0\leq t\leq{a_{1}(z_{\rm
max})}-{a_{1}(0)},$

where $t\in\mathbb{R}$ is an auxiliary optimization variable,
$Z_{\mathrm{min}}$ is given by:

$Z_{\mathrm{min}}=\begin{dcases}\cos(\theta_{1})\max\Big{\\{}\frac{z_{\mathrm{min}}}{\cos(\theta_{1})},r_{2}(\mathbf{q}_{2}),\frac{b_{\bot}(\mathbf{q}_{2})h_{\mathrm{min}}}{\lambda\sin(\theta_{1})}\Big{\\}},&\text{if
}b_{\mathrm{min}}\leq b_{\bot}(\mathbf{q}_{2}),\vspace{5 mm}\\\
\\!\begin{aligned}
&\cos(\theta_{1})\max\Big{\\{}\frac{z_{\mathrm{min}}}{\cos(\theta_{1})},r_{2}(\mathbf{q}_{2}),\frac{b_{\bot}(\mathbf{q}_{2})h_{\mathrm{min}}}{\lambda\sin(\theta_{1})},r_{2}(\mathbf{q}_{2})\cos\big{(}\theta_{1}-\theta_{2}(\mathbf{q}_{2})\big{)}+\\\
&\sqrt{b^{2}_{\mathrm{min}}-b^{2}_{\bot}(\mathbf{q}_{2})}\Big{\\}},\end{aligned}&\text{otherwise,
}\end{dcases}$ (33)

$Z_{\mathrm{max}}$ is given by:

$Z_{\mathrm{max}}=\begin{dcases}\cos(\theta_{1})\max\Big{\\{}\frac{z_{\mathrm{max}}}{\cos(\theta_{1})},\frac{2\pi
b_{\bot}(\mathbf{q}_{2})\Delta
h^{\mathrm{max}}}{\lambda\sin(\theta_{1})\Delta\phi_{90\%}(\gamma_{\mathrm{SNR}}^{\mathrm{min}}\gamma_{\mathrm{Rg}}^{\mathrm{min}}\gamma_{\mathrm{other}})}\Big{\\}},&\text{if
}b_{\mathrm{min}}\leq b_{\bot}(\mathbf{q}_{2}),\\\ \\!\begin{aligned}
&\cos(\theta_{1})\max\Big{\\{}\frac{z_{\mathrm{max}}}{\cos(\theta_{1})},r_{2}(\mathbf{q}_{2})\cos\big{(}\theta_{1}-\theta_{2}(\mathbf{q}_{2})\big{)}-\sqrt{b^{2}_{\mathrm{min}}-b^{2}_{\bot}(\mathbf{q}_{2})},\\\
&\frac{2\pi b_{\bot}(\mathbf{q}_{2})\Delta
h^{\mathrm{max}}}{\lambda\sin(\theta_{1})\Delta\phi_{90\%}(\gamma_{\mathrm{SNR}}^{\mathrm{min}}\gamma_{\mathrm{Rg}}^{\mathrm{min}}\gamma_{\mathrm{other}})}\Big{\\}},\end{aligned}&\text{otherwise.
}\end{dcases}$ (34)

Furthermore, function $a_{1}$ is given by:

$\displaystyle a_{1}(z_{1})$
$\displaystyle=\Big{(}2^{A_{1}z_{1}+A_{2}}-1\Big{)}\Big{(}z_{1}^{2}(\tan^{2}(\theta_{1})+1)+(g_{x}-x_{t})^{2}+$
$\displaystyle g_{z}^{2}+$ $\displaystyle
2z_{1}\tan(\theta_{1})[g_{x}-x_{t}]^{+}\Big{)}+\max\limits_{n}\Big{(}P_{\mathrm{com},1}[n]{\beta_{c,1}}\Big{)}+$
$\displaystyle\max\limits_{1\leq n\leq
N}\Big{\\{}\big{(}2^{A_{1}z_{1}+A_{2}}-1\big{)}(y[n]-g_{y})^{2}\Big{\\}},$
(35)

where
$A_{1}=\frac{2n_{B}B_{\mathrm{Rg}}\mathrm{PRF}}{cB_{c,1}}\left(\frac{1}{\cos\left(\theta_{1}+\frac{\Theta_{\mathrm{3dB}}}{2}\right)}-\frac{1}{\cos\left(\theta_{1}-\frac{\Theta_{\mathrm{3dB}}}{2}\right)}\right)$,
$A_{2}=\frac{n_{B}B_{\mathrm{Rg}}\mathrm{PRF}\tau_{p}}{B_{c,1}}$, and function
$a_{2}$ is given by:

$\displaystyle{a_{2}(z_{1})}=\Big{(}2^{A_{1}z_{1}+A_{2}}-1\Big{)}\Big{(}2z_{1}g_{z}+2z_{1}\tan(\theta_{1})[x_{t}-g_{x}]^{+}\Big{)}$
(36)

.

###### Proof.

Please refer to Appendix A. ∎

According to Proposition 1, sub-problem $(\mathrm{P.1.b})$ can be written in
the canonical form of a monotonic optimization problem given by
$(\overline{\mathrm{P.1.b}})$,

|  |
---|---|---
(a) | (b) | (c)
Figure 4: 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.
Vertex $\mathbf{l}_{\rm max}$ is initially selected as the best vertex and its
projection onto the normal set, denoted by $\Phi(\mathbf{l}_{\rm max})$, is
calculated. (b) A new polyblock is constructed and the next best vertex
$\mathbf{l}^{(2)}=\operatorname*{arg\,max}_{\mathbf{l}\in\mathcal{T}^{(2)}}\\{C_{N}(\mathbf{l},\mathbf{q}_{2},\mathbf{v}_{y})\\}$
is selected for projection. (c) Illustration of a tighter polyblock obtained
after several iterations and the $\epsilon_{1}$-optimal solution
$\mathbf{l}^{*}$ is located.

where the feasible set is the intersection of normal set $\mathcal{G}$,
spanned by constraints
($\overline{\mathrm{C1a}},\mathrm{C6},\overline{\mathrm{C11a}},\overline{\mathrm{C11c}}$),
and conormal set $\mathcal{H}$, spanned by constraints
$\rm(\overline{C1b},\overline{C11b})$. The optimal solution of sub-problem
${(\overline{\mathrm{P.1.b}}})$ lies at the upper boundary of normal set
$\mathcal{G}$, denoted by $\partial^{+}\mathcal{G}$ [40], see Figure 4a. This
region is not known a priori but the sequential polyblock approximation
algorithm can be used to approach this set from above.

Algorithm 2 Polyblock Outer Approximation Algorithm

1:Initialization: Set polyblock $\mathcal{P}^{(1)}$ with vertex
$\mathbf{l}^{(1)}=\mathbf{l}_{\rm max}=(z_{\rm max},t_{\rm max})^{T}$, vertex
set $\mathcal{T}^{(1)}=\\{\mathbf{l}_{\rm max}\\}$, small positive number
$\epsilon_{1}\geq 0$, objective function value $\rm CBV^{(1)}=-\infty$,
current solution $\mathbf{s}^{(1)}=[\cdot]$, and iteration number $j=1$.

2:Remove from $\mathcal{T}^{(1)}$ vertices
$\\{\mathbf{l}\in\mathcal{T}^{(1)}|\mathbf{l}\notin\mathcal{H}\\}$

3:while $\Big{(}\mathcal{T}^{(j)}\neq\emptyset$ and
$\Big{|}C_{N}(\mathbf{l}^{(j)},\mathbf{q}_{2},\mathbf{v}_{y})-\mathrm{CBV}^{(j)}\Big{|}>\epsilon_{1}\Big{)}$

4: Select
$\mathbf{l}^{(j)}=\operatorname*{arg\,max}\limits_{\mathbf{l}\in\mathcal{T}^{(j)}}\Big{\\{}C_{N}(\mathbf{l},\mathbf{q}_{2},\mathbf{v}_{y})\Big{\\}}$

5: Compute the projection $\Phi(\mathbf{l}^{(j)})$ onto
$\partial\mathcal{G}^{+}$

6: using Algorithm 3

7: if $\Phi(\mathbf{l}^{(j)})\in\mathcal{H}$ and
$C_{N}(\Phi(\mathbf{l}^{(j)}),\mathbf{q}_{2},\mathbf{v}_{y})>\mathrm{CBV}^{(j)}$

8: then

9: Let the best solution $\textbf{s}^{(j+1)}=\Phi(\mathbf{l}^{(j)})$ and

10:
$\mathrm{CBV}^{(j+1)}=C_{N}(\Phi(\mathbf{l}^{(j)}),\mathbf{q}_{2},\mathbf{v}_{y})$

11: else

12: Set $\textbf{s}^{(j+1)}=\textbf{s}^{(j)}$ and
$\mathrm{CBV}^{(j+1)}=\mathrm{CBV}^{(j)}$

13: endif

14: Set
$\mathcal{T}^{(j+1)}=\\{\mathcal{T}^{(j)}\setminus\mathbf{l}^{(j)}\\}\cup\Big{\\{}\Big{(}\Phi(\mathbf{l}^{(j)})[1],l^{(j)}[2]\Big{)}^{T}\Big{\\}}\cup$

15: $\Big{\\{}\Big{(}l^{(j)}[1],\Phi(\mathbf{l}^{(j)})[2]\Big{)}^{T}\Big{\\}}$

16: Remove from $\mathcal{T}^{(j+1)}$ vertices
$\\{\mathbf{l}\in\mathcal{T}^{(j+1)}|\mathbf{l}\notin\mathcal{H}\\}$

17: Set $j=j+1$

18:end while

19:return solution
$\mathbf{q}_{1}=\Big{(}\underbrace{x_{t}-s^{(j)}[1]\tan(\theta_{1})}_{x_{1}},\underbrace{s^{(j)}[1]}_{z_{1}}\Big{)}^{T}$

Algorithm 3 Bisection Search Algorithm to Compute the Upper Boundary Point

1:Initialization: Given $\mathbf{l}^{(j)}$, set $\lambda_{\rm min}=0$,
$\lambda_{\rm max}=1$, and error tolerance $0<{\epsilon_{2}}\ll 1$.

2:repeat

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

4: if $\lambda\mathbf{l}^{(j)}$ is feasible for
($\overline{\mathrm{C1a}},\mathrm{C6},\overline{\mathrm{C11a}},\overline{\mathrm{C11c}}$)
then

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

6: else$\triangleright$ If $\lambda\mathbf{l}^{(j)}\notin\mathcal{G}$

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_{2}}$$\triangleright$ Convergence if $\lambda_{\rm
max}\approx\lambda_{\rm min}$

10:return $\Phi(\mathbf{l}^{(j)})=\lambda\mathbf{l}^{(j)}$

In Algorithm 2, we provide all steps of the proposed polyblock outer
approximation algorithm that generates an $\epsilon_{1}$-approximate optimal
solution to sub-problem $(\overline{\mathrm{P.1.b}})$. In addition, Figure 4
depicts the operation of the polyblock algorithm. In fact, optimization
variables $z_{1}$ and $t$ of sub-problem $(\overline{\mathrm{P.1.b}})$ are
collected in vertex vector
$\mathbf{l}=(l[1],l[2])^{T}=\big{(}z_{1},t\big{)}^{T}\in\mathbb{R}^{2}$.
Hereinafter, for ease of notation,
$C_{N}(\mathbf{l},\mathbf{q}_{2},\mathbf{v}_{y})$ refers to
$C_{N}((x_{t}-l[1]\tan(\theta_{1}),l[1])^{T},\mathbf{q}_{2},\mathbf{v}_{y})$.
As shown in Figure 4a, first, we initialize a polyblock $\mathcal{P}^{(1)}$
with corresponding vertex set $\mathcal{T}^{(1)}$ containing vertex
$\mathbf{l}_{\rm max}=(z_{\rm max},t_{\rm max})^{T}\in\mathbb{R}^{2}$, where
$t_{\rm max}={a_{1}(z_{\rm max})-a_{1}(0)}$. Then, in iteration $j$, the
vertex that maximizes the objective function of sub-problem
$(\overline{\mathrm{P.1.b}})$ is selected to create the next polyblock, i.e.,
$\mathbf{l}^{(j)}=\operatorname*{arg\,max}\limits_{\mathbf{l}\in\mathcal{T}^{(j)}}\\{C_{N}(\mathbf{l},\mathbf{q}_{2},\mathbf{v}_{y})\\}$.
The projection of vertex $\mathbf{l}^{(j)}$ onto the normal set $\mathcal{G}$,
denoted by $\Phi(\mathbf{l}^{(j)})$, is calculated and used to construct a
tighter polyblock $\mathcal{P}^{(j)}$, see Figure 4b. The new vertex set is
calculated as
$\mathcal{T}^{(j)}=\\{\mathcal{T}^{(j)}\setminus\mathbf{l}^{(j)}\\}\cup\left\\{\Big{(}\Phi(\mathbf{l}^{(j)})[1],l^{(j)}[2]\Big{)}^{T}\right\\}\cup\left\\{\Big{(}l^{(j)}[1],\Phi(\mathbf{l}^{(j)})[2]\Big{)}^{T}\right\\}$.
Once a feasible projection (i.e.,
$\Phi(\mathbf{l}^{(j)})\in\mathcal{G}\cap\mathcal{H}$) that improves the
objective function of sub-problem $(\overline{\mathrm{P.1.b}})$ is located,
the corresponding vector $\mathbf{s}^{(j+1)}=\Phi(\mathbf{l}^{(j)})$ and its
objective function
$\mathrm{CBV}^{(j+1)}=C_{N}(\Phi(\mathbf{l}^{(j)}),\mathbf{q}_{2},\mathbf{v}_{y})$
are recorded. To speed up convergence and save memory, vertices that are not
in set $\mathcal{H}$ are removed in each iteration [41]. 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}^{(j)}\supset\mathcal{G}\cap\mathcal{H}$,
see Figure 4c. The termination criteria of the proposed algorithm are (i)
$|C_{N}(\mathbf{l}^{(j)},\mathbf{q}_{2},\mathbf{v}_{y})-\mathrm{CBV}^{(j)}|\leq\epsilon_{1}$,
for convergence, where $\epsilon_{1}$ is the error tolerance, and (ii)
$\mathcal{T}^{(j)}=\emptyset$ for a non-feasible problem.

In general, Algorithm 2 converges to an $\epsilon_{1}$-optimal solution with
an exponential time complexity in the dimension of the vertex [41].
Fortunately, for the problem at hand, the vertex dimension is only 2,
therefore, the optimal solution is obtained with low computational complexity.
In particular, the time complexity of Algorithm 2 stems mainly from the
projection performed by Algorithm 3. In Algorithm 3, a bisection search in the
interval $[0,1]$ is conducted with precision ${\epsilon_{2}}$, thus, the time
complexity is $\mathcal{O}({\log_{2}}(1/{\epsilon_{2}}))$. Thus, the
complexity of Algorithm 2 is $\mathcal{O}(M_{2}{\log_{2}}(1/{\epsilon_{2}}))$,
where $M_{2}$ is the maximum number of iterations required for convergence.

###  IV-C Resource Allocation Optimization

In this subsection, we aim to optimize
$\\{\mathbf{v}_{y},\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,2}}\\}$
given uav formation $\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}$. Then, problem
$\mathrm{(P.1)}$ reduces to sub-problem $\mathrm{(P.1.c)}$ given by:

$\displaystyle(\mathrm{P.1.c}):\max_{\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,2}},\mathbf{v}_{y}}\hskip
8.53581ptC_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})$ s.t.
$\displaystyle\mathrm{C6},\mathrm{C10},\mathrm{C11},\mathrm{C12},\mathrm{C13}.$

Sub-problem $\mathrm{(P.1.c)}$ is difficult to solve due to the non-convex
propulsion power model used in constraint $\mathrm{C12}$. Yet, we provide a
low-complexity sub-optimal solution to $\mathrm{(P.1.c)}$ based on successive
convex approximation (SCA). To start with, we introduce a slack vector
$\mathbf{u}=(u[1],...,u[N])^{T}\in\mathbb{R}^{N}$ and replace the propulsion
power equivalently by:

$\displaystyle\widetilde{P_{\mathrm{prop},n}}(\mathbf{v}_{y},\mathbf{u})=P_{0}\left(1+\frac{3v_{y}^{2}[n]}{U^{2}_{\mathrm{tip}}}\right)+u[n]+\frac{1}{2}d_{0}\rho
sAv_{y}^{3}[n],\forall n,$ (37)

where the elements of slack vector $\mathbf{u}$ must satisfy the following
additional constraint:

$\mathrm{C12a}:1+\frac{v^{4}_{y}[n]}{4v_{0}^{4}}-\Big{(}\frac{u^{2}[n]}{P^{2}_{I}}+\frac{v^{2}_{y}[n]}{2v^{2}_{0}}\Big{)}^{2}\leq
0,\forall n.$ (38)

After replacing $P_{\mathrm{prop},n}$ with $\widetilde{P_{\mathrm{prop},n}}$,
constraint $\mathrm{C12}$ is denoted by $\widetilde{\mathrm{C12}}$. However,
constraint $\mathrm{C12a}$ is still non-convex due to the term
$f(\mathbf{v}_{y},\mathbf{u})=\Big{(}\frac{u^{2}[n]}{P^{2}_{I}}+\frac{v^{2}_{y}[n]}{2v^{2}_{0}}\Big{)}^{2}$.
Therefore, in iteration $i$ of the sca algorithm, we approximate function $f$
based on its Taylor series expansion around point
$\\{\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)}\\}$ as follows:

$\widetilde{f}(\mathbf{v}_{y},\mathbf{u})=f(\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)})+\nabla_{\mathbf{v}_{y}}f(\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)})^{T}(\mathbf{v}_{y}-\mathbf{v}_{y}^{(i)})+\nabla_{\mathbf{u}}f(\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)})^{T}(\mathbf{u}-\mathbf{u}^{(i)}),$
(39)

where
$\nabla_{\mathbf{v}_{y}}f(\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)})=\frac{\partial
f}{\partial\mathbf{v}_{y}}(\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)})$ and
$\nabla_{\mathbf{u}}f(\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)})=\frac{\partial
f}{\partial\mathbf{u}}(\mathbf{v}_{y}^{(i)},\mathbf{u}^{(i)})$. Function
$\widetilde{f}(\mathbf{v}_{y},\mathbf{u})$ is inserted in constraint
$\mathrm{C12a}$ leading to a new convex constraint denoted by
$\widetilde{\mathrm{C12a}}$. The resulting convex optimization problem is
given by:

$\displaystyle(\widetilde{\mathrm{P.1.c}})$
$\displaystyle:\max_{\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,2}},\mathbf{v}_{y},\mathbf{u}}\hskip
8.53581ptC_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y})$
$\displaystyle\text{s.t.
}\mathrm{C6},\mathrm{C10},\mathrm{C11},\widetilde{\mathrm{C12}},\widetilde{\mathrm{C12a}},\mathrm{C13}.$

Algorithm 4 Successive Convex Approximation Algorithm

1:For fixed $\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}$, set initial point
$\\{\mathbf{v}_{y}^{(1)},\mathbf{u}^{(1)}\\}$, iteration index $i=1$, and
error tolerance $0<{\epsilon_{3}}\ll 1$.

2:repeat

3: Determine coverage
$C_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y}),\mathbf{u},\mathbf{v}_{y},\mathbf{P}_{\mathrm{com,1}},$
and $\mathbf{P}_{\mathrm{com,2}}$ by solving $\mathrm{\widetilde{(P.1.c)}}$
around point

4: $\\{\mathbf{u}^{(i)},\mathbf{v}_{y}^{(i)}\\}$ using CVXPY

5: Set $i=i+1,\hskip 2.84526pt\mathbf{u}^{(i)}=\mathbf{u},\hskip
2.84526pt\mathbf{v}_{y}^{(i)}=\mathbf{v}_{y}$ $\triangleright$ Update solution
and iteration index

6:until
$\big{|}\frac{C_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y}^{(i)})-C_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y}^{(i-1)})}{C_{N}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{v}_{y}^{(i)})}\big{|}<{\epsilon_{3}}$

7:return solution
{$\mathbf{v}_{y},\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,2}}$}

The proposed procedure to solve sub-problem $\mathrm{(P.1.c)}$ is summarized
in Algorithm 4, where the convex approximation $\mathrm{(\widetilde{P.1.c})}$
is solved using the Python convex optimization library CVXPY [42]. It can be
shown that Algorithm 4 converges to a local optimum of sub-problem
$\mathrm{(P.1.c)}$ in polynomial computational time complexity [43]. In
particular, Algorithm 4 involves $3N-1$ optimization variables and, therefore,
has a computational time complexity of $\mathcal{O}(M_{3}(3N-1)^{3.5})$, where
$M_{3}$ is the number of iterations needed for convergence [44].

###  IV-D Solution of Problem $\rm(P.1)$

To solve problem $\rm(P.1)$, we use ao by solving sub-problems $\rm(P.1.a)$,
$\rm(P.1.b)$, and $\rm(P.1.c)$, alternately, see Figure 3. To guarantee the
convergence of the proposed ao Algorithm 5, Algorithm 1 has to be adjusted in
each iteration $m$ to ensure a non-decreasing objective function. This can be
achieved by adding one more particle to the initial pso population of size
$D$. The added particle represents the position of the slave uav in the
across-track plane position in the previous ao iteration, i.e.,
$\mathbf{q}_{2}^{(m-1)}$, and the resulting pso population is now of size
$D+1$. Consequently, Algorithm 5 is guaranteed to converge due to the fact
that the objective function is (i) non-decreasing and (ii) can be shown to be
upper-bounded by $\big{(}S_{\rm far}(\mathbf{q}_{\rm max})-S_{\rm
near}(\mathbf{q}_{\rm max})\big{)}(N-1)v_{\rm max}\delta_{t}$ , where
$\mathbf{q}_{\rm max}=(x_{t}-z_{\rm max}\tan(\theta_{1}),z_{\rm max})^{T}$
[45]. However, for complex optimization problems, characterized by multiple
local optima, ao methods are vulnerable to getting trapped in local optima
[45], which is the case for problem $\mathrm{(P.1)}$ due to the coupled
optimization variables. In particular, the objective function is increasing
wrt the uav velocity $\mathbf{v}_{y}$. However, increasing $\mathbf{v}_{y}$
results in constraining the uav formation
$\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}$ due to the limited snr decorrelation
(see constraint $\mathrm{C6}$), yielding multiple local optima depending on
the magnitude of the ao variable updates in each iteration. To avoid premature
convergence to low-quality local optima, we propose to limit the uav velocity
update in iteration $m$ of the ao algorithm by using a step size, denoted by
${\psi}$, as follows:

$\mathbf{v}^{(m)}_{y}=\begin{dcases}\mathbf{v}_{y}^{(m-1)}+{\psi}(\mathbf{v}_{y}-\mathbf{v}_{y}^{(m-1)}),\hskip
5.69054pt&\text{ if feasible for $\mathrm{(P.1.c)}$, }\\\
\mathbf{v}_{y},\hskip 5.69054pt&\text{ otherwise, $\forall m.$ }\end{dcases}$
(40)

where $\mathbf{v}^{(m)}_{y}$ is the uav velocity update in the $m^{\rm th}$
iteration of the ao algorithm, and $\\{\mathbf{v}_{y},\mathbf{P}_{\rm
com,1},\mathbf{P}_{\rm com,2}\\}$ is the solution of $\mathrm{(P.1.c)}$
generated by Algorithm 4. Note that small ${\psi}$ values result in smaller
updates for the velocity vector, whereas ${\psi}=1$ results in the classical
ao updates [45].

Algorithm 5 $\psi$-Step Alternating Optimization Algorithm

1:Set initial formation $\\{\mathbf{q}_{1}^{(1)},\mathbf{q}_{2}^{(1)}\\}$ and
initial resources
$\\{\mathbf{v}^{(1)}_{y},\mathbf{P}^{(1)}_{\mathrm{com,1}},\mathbf{P}^{(1)}_{\mathrm{com,2}}\\}$,
iteration index $m=1$, and error tolerance $0<{\epsilon_{4}}\ll 1$.

2:repeat

3: Set $m=m+1$

4: Determine coverage
$C_{N}(\mathbf{q}_{1}^{(m-1)},\mathbf{q}_{2},\mathbf{v}_{y}^{(m-1)})$ and set
$\mathbf{q}_{2}^{(m)}=\mathbf{q}_{2}$ by solving $\rm(P.1.a)$ for fixed

5:
$\\{\mathbf{v}_{y}^{(m-1)},\mathbf{P}^{(m-1)}_{\mathrm{com,1}},\mathbf{P}^{(m-1)}_{\mathrm{com,2}},\mathbf{q}_{1}^{(m-1)}\\}$
using Algorithm 1

6: Determine coverage
$C_{N}(\mathbf{q}_{1},\mathbf{q}_{2}^{(m)},\mathbf{v}_{y}^{(m-1)})$ and set
$\mathbf{q}_{1}^{(m)}=\mathbf{q}_{1}$ by solving $\rm(P.1.b)$ for fixed

7:
$\\{\mathbf{v}_{y}^{(m-1)},\mathbf{P}^{(m-1)}_{\mathrm{com,1}},\mathbf{P}^{(m-1)}_{\mathrm{com,2}},\mathbf{q}_{2}^{(m)}\\}$
using Algorithm 2

8: Determine coverage
$C_{N}(\mathbf{q}_{1}^{(m)},\mathbf{q}_{2}^{(m)},\mathbf{v}_{y})$ and
$\\{\mathbf{v}_{y},\mathbf{P}_{\mathrm{com,1}},\mathbf{P}_{\mathrm{com,1}}\\}$
by solving $\rm(P.1.c)$ for fixed

9: formation $\\{\mathbf{q}_{1}^{(m)},\mathbf{q}_{2}^{(m)}\\}$ using Algorithm
4

10: if
$\\{\mathbf{v}_{y}^{(m-1)}+{\psi}(\mathbf{v}_{y}-\mathbf{v}_{y}^{(m-1)}),\mathbf{P}_{\rm
com,1},\mathbf{P}_{\rm com,2}\\}$ is feasible for $\mathrm{(P.1.c)}$ then

11: Set
$\mathbf{v}_{y}^{(m)}=\mathbf{v}_{y}^{(m-1)}+{\psi}(\mathbf{v}_{y}-\mathbf{v}_{y}^{(m-1)}),$
$\mathbf{P}^{(m)}_{\mathrm{com,1}}=\mathbf{P}_{\mathrm{com,1}},$ and
$\mathbf{P}^{(m)}_{\mathrm{com,2}}=\mathbf{P}_{\mathrm{com,2}}$

12: else

13: Set $\mathbf{v}_{y}^{(m)}=\mathbf{v}_{y},$
$\mathbf{P}^{(m)}_{\mathrm{com,1}}=\mathbf{P}_{\mathrm{com,1}},$ and
$\mathbf{P}^{(m)}_{\mathrm{com,2}}=\mathbf{P}_{\mathrm{com,2}}$

14:until
$\big{|}\frac{C_{N}(\mathbf{q}_{1}^{(m)},\mathbf{q}_{2}^{(m)},\mathbf{v}_{y}^{(m)})-C_{N}(\mathbf{q}_{1}^{(m-1)},\mathbf{q}_{2}^{(m-1)},\mathbf{v}_{y}^{(m-1)})}{C_{N}(\mathbf{q}_{1}^{(m)},\mathbf{q}_{2}^{(m)},\mathbf{v}_{y}^{(m)})}\big{|}\leq{\epsilon_{4}}$

15:return solution
$\\{\mathbf{q}_{1}^{(m)},\mathbf{q}_{2}^{(m)},\mathbf{P}^{(m)}_{\mathrm{com,1}},\mathbf{P}^{(m)}_{\mathrm{com},2},\mathbf{v}_{y}^{(m)}\\}$

In Algorithm 5, we summarize all the steps of the solution to problem
$\mathrm{(P.1)}$, see Figure 3. Leveraging the computational complexity
analysis of the different sub-problems, Algorithm 5 has
$\mathcal{O}\Big{(}M\big{(}M_{1}DN+M_{2}{\log_{2}}(1/{\epsilon_{2}})+M_{3}(3N-1)^{3.5}\big{)}\Big{)}$
worst-case computational time complexity, where $M$ is the number of
iterations required for Algorithm 5 to converge. In practice, convergence is
achieved in just a few iterations. Consequently, due to the low complexity of
Algorithm 5, the $\epsilon_{5}$-optimal step size ${\psi}$ can be derived by
performing an exhaustive search in $[0,1]$ with precision $\epsilon_{5}$. Note
that pso Algorithm 1, being a stochastic approach, can achieve the same
objective function value for sub-problem $\mathrm{(P.1.a)}$ for different
solutions. Therefore, Algorithm 5 might lead to slightly different results for
the overall problem for $\mathrm{(P.1)}$ given the same input, which is a
known property of ao algorithm when involving a stochastic algorithm[45].

TABLE I: System parameters [31, 26, 17, 19, 6].

Parameter  Value Parameter Value [HTML]EFEFEFuav parameters  [HTML]EFEFEFRadar
parameters  Total number of time slots, $N$ 80  Reference target
$x$-coordinate, $x_{t}$ 20 m  Maximum and minimum uav altitudes,
$z_{\mathrm{max}}$ and $z_{\mathrm{min}}$ 100 and 1 m  Time slot and pulse
duration, $\delta_{t}$ and $\tau_{p}$ 1 and $10^{-6}$ s  Maximum and minimum
uav velocities, $v_{\mathrm{max}}$ and $v_{\mathrm{min}}$ 10 and 0.1 m/s
Master look angle and radar beamwidth, $\theta_{1}$ and
$\Theta_{3\mathrm{dB}}$ 45° and 30°  UAV battery capacity, $E_{\mathrm{max}}$
12.22 Wh  Minimum and maximum slave look angle, $\theta_{\rm min}$ and
$\theta_{\rm max}$ 15° and 75°  Minimum baseline distance, $b_{\mathrm{min}}$
2 m  Radar transmit and recieve gains, $G_{t}$ and $G_{r}$ 6 dBi  UAV weight
in Newton, $W_{u}$ 60 N  Radar transmit power and noise figure,
$P_{t,1}=P_{t,2}$ and $F$ 10 dBm and 5 dB  Air density, $\rho$ 1.225 kg/m3
System, azimuth, and atmospheric losses, $L_{\mathrm{sys}}$,
$L_{\mathrm{azm}}$, and $L_{\mathrm{atm}}$ 2, 2, and 0 dB  Fuselage drag
ratio, $d_{0}$ 0.6  Pulse repetition frequency, $\mathrm{PRF}$ 100 Hz  Profile
drag coefficient, $\delta_{u}$ 0.0012  Radar wavelength, $\lambda$ 0.12 m
Rotor radius, $R$ 0.4 m  Radar bandwidth and center frequency,
$B_{\mathrm{Rg}}$ and $f_{0}$ 3 and 2.5 GHz  Rotor disc area, $A_{e}$ 0.503 m2
Number of independent looks, $n_{L}$ 16  Rotor solidity ratio, $s$ 0.05  Radar
system temperature, $T_{\mathrm{sys}}$ 400 K  Blade angular velocity, $\Omega$
300 rad/s  Normalized backscatter coefficient, $\sigma_{0}$ -10 dB  Tip speed
of the rotor blade, $\rho$ 120 m/s  Minimum snr and baseline
decorrelation,$\gamma_{\mathrm{SNR}}^{\mathrm{min}}$ and
$\gamma_{\mathrm{Rg}}^{\mathrm{min}}$ 0.8  Incremental correction factor to
induced power, $k_{u}$ 0.1  Other decorrelation sources, $\gamma_{\rm other}$
0.8  [HTML]EFEFEFCommunication parameters  Minimum hoa and maximum height
error, $h_{\mathrm{amb}}^{\mathrm{min}}$ and $\Delta_{h}^{\rm max}$ 1 and
0.224 m  gs 3D location, $\mathbf{g}=(g_{x},g_{y},g_{z})^{T}$ $(0,-270,5)^{T}$
m  [HTML]EFEFEFAlgorithm parameters  Master communication bandwidth, $B_{c,1}$
1 GHz  pso population and maximum number of iterations, $D$ and $M_{1}$
$2\times 10^{3}$ and $200$ Slave communication bandwidth, $B_{c,2}$ 1 GHz  pso
cognitive and social constants, $c_{1}$ and $c_{2}$ 0.1 and 0.2  Communication
reference channel gains, ${\beta_{c,1}}={\beta_{c,2}}$ 18.75 dB  pso maximum
particle velocity, $v_{\rm PSO}^{\rm max}$ 20 m/s  Number of bits per complex
sample, $n_{B}$ 4  Constant $O$ and number of realizations $\mathcal{N}$ 500
and $10^{3}$ Maximum communication transmit power,
$P_{\mathrm{com}}^{\mathrm{max}}$ 10 dB  Error tolerance,
${\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=\epsilon_{4}}$ and $\epsilon_{5}$
$10^{-4}$ and $10^{-2}$

##  V Simulation Results and Discussion

In this section, we present simulation results for the proposed uav formation
and resource allocation optimization algorithm. Unless specified otherwise,
the system parameters specified in Table I are employed. As explained in the
previous section, since the algorithm proposed to solve $\mathrm{(P.1)}$ is
based on stochastic optimization, different solutions may be obtained for the
same input. Therefore, simulation results are averaged over $10^{3}$
realizations. To assess the performance of the proposed algorithm, we consider
the following benchmark schemes:

  * •

Benchmark scheme 1 (Classical-alternating-optimization): For this scheme, we
employ the classical ao algorithm [45] to solve problem $\mathrm{(P.1)}$,
i.e., without adjusting step size ${\psi}$. This is equivalent to using
Algorithm 5 with ${\psi}=1$.

  * •

Benchmark scheme 2 (Fixed-steady-speed): As in [1], in this scheme, we employ
a fixed and steady radar platform velocity, i.e., $\mathbf{v}_{y}$ is not
optimized but is set to $\mathbf{v}_{y}[n]=4$ m/s$,\forall n$. We use ao to
optimize the communication powers and the uav formation. This is equivalent to
using Algorithm 5 with ${\psi}=0$.

  * •

Benchmark scheme 3 (Fixed-slave-look-angle): Similar to [1], we adopt a
predefined look angle for the slave platform,
$\theta_{2}(\mathbf{q}_{2})=\frac{\pi}{4}$. Problem ($\mathrm{P.1}$) is solved
with a corresponding modified version of Algorithm 5. Here, the value of
$\psi$ leading to the best insar coverage is searched via exhaustive search in
$[0,1]$ with precision $\epsilon_{5}$.

Figure 5:  Interferometric sensing coverage achieved by the proposed solution
and benchmark schemes versus iteration number for different initial points,
denoted by $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$, and given by
$\mathcal{F}_{1}=\Big{\\{}\mathbf{q}^{(1)}_{1}=(-40,60)^{T}$ m,
$\mathbf{q}^{(1)}_{2}=(-45,50)^{T}$ m, $P^{(1)}_{\rm com,1}[n]=P^{(1)}_{\rm
com,2}[n]=7.78$ dB, $v^{(1)}_{y}[n]=4$ m/s, $\forall n\Big{\\}}$ and
$\mathcal{F}_{2}=\Big{\\{}\mathbf{q}^{(1)}_{1}=(-20,40)^{T}$ m,
$\mathbf{q}^{(1)}_{2}=(-30,40)^{T}$ m, $P^{(1)}_{\rm com,1}[n]=P^{(1)}_{\rm
com,2}[n]=7.78$ dB, $v^{(1)}_{y}[n]=4$ m/s, $\forall n\Big{\\}}$. Solid lines
represent mean values while the shaded area corresponds to one standard
deviation, i.e., $1\sigma$. The ${\psi}$ values are optimized based on
exhaustive search.

In Figure 5, we investigate the convergence of the proposed ao algorithm and
the different benchmark schemes by plotting the total insar coverage as a
function of the number of iterations. The figure shows that the classical ao
algorithm, i.e., benchmark scheme 1, outperforms benchmark schemes 2 and 3 due
to the joint optimization of the uav formation and resource allocation. In
fact, benchmark scheme 2 exhibits the lowest performance, underlining the
importance of optimizing the instantaneous velocity of both uav and the slave
look angle to achieve maximal sensing coverage, respectively. Furthermore,
Figure 5 shows that, after convergence, the proposed solution outperforms all
benchmark schemes achieving performance gains of 14.8 %, 41.21 %, and 35.74 %,
compared to benchmark schemes 1, 2, and 3, respectively. In particular, the
superiority of the proposed solution compared to benchmark scheme 1 is
attributed to the tuning of the step size ${\psi}$ to avoid premature
convergence to inferior sub-optimal solutions. This performance gain comes at
the cost of a lower convergence speed, as the proposed solution requires more
iterations to converge. Nevertheless, overall, all benchmark schemes as well
as the proposed solution converge after a relatively small number of
iterations. Lastly, for different initial points $\mathcal{F}_{1}$ and
$\mathcal{F}_{2}$, the proposed ao-based algorithm achieves similar
performance confirming resilience to the initialization.

Figure 6: Interferometric sensing coverage versus step size ${\psi}$ of
Algorithm 5 for different initial points, denoted by $\mathcal{F}_{1}$ and
$\mathcal{F}_{2}$, see Figure 5. Solid lines represent mean values while the
shaded area corresponds to one standard deviation, i.e., $1\sigma$.

In Figure 6, we show the total insar coverage as a function of step size
${\psi}$ defined in (40). The figure confirms that the performance of the
conventional ao algorithm, i.e., for ${\psi}=1$, can be significantly improved
by tuning ${\psi}$ in $[0,1]$. Moreover, the optimal step size depends on the
system parameters and the ao initialization process. For instance, the best
step size for initial point $\mathcal{F}_{1}$ is ${\psi}={0.43}$, while for
$\mathcal{F}_{2}$, it is ${\psi}={0.36}$. Figure 6 reveals that values of
${\psi}$ close to 1 and 0 are not preferable. On the one hand, a value of
${\psi}$ close to 1 results in a fast increase of the velocity in the first
few iterations to improve coverage. However, this leads to low snr for the
master and slave uav, i.e., a high snr decorrelation, see (8). Therefore, the
uav formation in the across-track cannot be properly optimized and both drones
have to be positioned close to the target area to keep short master and slave
slant ranges that compensate for the impact of the high velocity on the snr.
On the other hand, a value of ${\psi}$ close to 0 results in slow velocity
updates, which allows for optimizing the uav formation to improve coverage by
opting for large master and slave slant ranges. However, the velocity can
subsequently not be increased due to the snr decorrelation constraint
$\mathrm{C6}$, resulting in poor performance. The proposed approach overcomes
this issue by adequately selecting ${\psi}$ so that the updates of the
different optimization variables are well-balanced in each iteration to
achieve maximum sensing coverage.

Figure 7: Optimal uav formation represented by the uav baseline and
perpendicular baseline versus the maximum communication transmit power for
${\beta_{c,1}}={\beta_{c,2}}=19.3$ dB.

In Figure 7, we investigate the impact of the maximum communication transmit
power $P_{\rm com}^{\rm max}$ on the optimal uav bistatic formation in terms
of the interferometric baseline $b$ and its perpendicular component
$b_{\bot}$. The figure shows that the uav formation employed for sensing is
heavily affected by the communication transmit power used to offload the sar
data. In particular, for low maximum communication transmit powers, i.e.,
$P_{\rm com}^{\rm max}\leq{9.8}$ dB (blue-colored region), the proposed
optimization algorithm prefers a close formation with short interferometric
baselines as small as 6 m. In contrast, high maximum communication
transmission powers, i.e., $P_{\rm com}^{\rm max}\geq{9.8}$ dB (green-colored
region), result in a large uav formation, with baselines as large as 14 m.
This is because for low maximum communication transmit powers, data offloading
is the performance bottleneck and, thus, Algorithm 5 tends to position both
uav near the target area, resulting in less sensing data. This leads to a
short master slant range $r_{1}$ and, therefore, a shorter perpendicular
baseline is required to satisfy constraint $\mathrm{C8}$, see (11). On the
other hand, when the maximum communication transmit power is sufficiently
high, the proposed solution optimizes the uav formation by flying at higher
altitudes for better coverage, leading to large uav formations. This results
in higher sensing data rates for the slave and master uav. Yet, all data from
both drones can be successfully offloaded because of the high communication
transmission power available.

Figure 8: Interferometric sensing coverage versus the minimum required snr
decorrelation $\gamma_{\rm SNR}^{\rm min}$ for different maximum communication
transmit powers $P_{\mathrm{com}}^{\mathrm{max}}$.

In Figure 8, we show the optimal insar coverage obtained with Algorithm 5 as a
function of the minimum required snr decorrelation for different maximum
communication transmit powers. The figure reveals that the total sensing
coverage decreases as the minimum required snr decorrelation increases, i.e.,
as the sar coherence requirement becomes more stringent. We distinguish
between three regions in Figure 8, namely Regions I, II and III. In Region I,
a high SNR decorrelation is allowed and increasing the maximum communication
transmit power results in better sensing coverage due to the extended range
for communications. Therefore, in Region I, the system performance bottleneck
is the air-to-ground communication. In Region II, stricter sensing
requirements in terms of sar coherence are imposed and the maximum
communication power has a more moderate impact on the coverage compared to
Region I. Finally, in Region III, a high coherence is required and constraint
$\mathrm{C6}$ becomes active, i.e., sensing becomes the system performance
bottleneck. Consequently, the maximum communication transmit power has little
impact on the coverage performance, which highlights the interplay between
communication and sensing in the considered system.

Figure 9:  Average uav velocity as function of the maximum communication
transmit power $P_{\rm com}^{\rm max}$ for different values of the minimum SNR
decorrelation $\gamma_{\rm SNR}^{\rm min}$ and for
$\beta_{c,1}=\beta_{c,2}=20.91$ dB.

In Figure 9, we depict the average uav velocity, given by
${\frac{1}{N-1}}\sum\limits_{n=1}^{N-1}v_{y}[n]$, as function of the maximum
communication transmit power for different minimum SNR decorrelation
thresholds. The figure shows that the average uav velocity increases with the
maximum communication transmit power due to the increased communication range
for both uav, which enables longer trajectories in azimuth direction.
Moreover, the figure shows a strong impact of the sensing requirement on the
drone speed. In fact, to achieve higher insar coherence, a lower speed is
preferred, as it allows for improved uav snr. This further confirms the need
for slowly moving drones for insar imaging applications, especially when
highly reliable sensing is demanded.

##  VI Conclusion

In this paper, we investigated joint formation and resource allocation
optimization for bistatic uav-based 3d insar sensing. To assess the 3d sensing
performance, the insar coverage, interferometric coherence, hoa, and relative
height accuracy were introduced as relevant insar-specific performance
metrics. A non-convex non-smooth optimization problem was formulated and
solved for the maximization of the bistatic insar ground coverage while
enforcing energy, communication, and sensing performance constraints. Our
simulation results confirmed the superior performance of the proposed
algorithm compared to several benchmark schemes and emphasized the relevance
of optimizing the velocity of both radar platforms and the look angle of the
slave uav for maximum insar coverage. We revealed the importance of low uav
speed for high insar performance and showed that large areas can be mapped and
offloaded to ground in real time while ground object heights can be estimated
with high precision. Lastly, we highlighted the significant impact of the uav-
gs communication link on the uav formation, and therefore, sensing
performance.

##  Appendix A Proof of Proposition 1

Constraint $\mathrm{C2}$ makes sub-problem $\mathrm{(P.1.b)}$ a one-
dimensional optimization problem wrt $z_{1}$, as it specifies a linear
relation between the $x$\- and $z$-coordinates of $U_{1}$. In fact, constraint
$\mathrm{C2}$ imposes that
$r_{1}(\mathbf{q}_{1})=\frac{z_{1}}{\cos(\theta_{1})}$, and therefore,
constraints $\mathrm{C1},\mathrm{C3},\mathrm{C5},$ and $\mathrm{C8}$ yield
lower and upper bounds on $z_{1}$ that are collected in variables $Z_{\rm
max}$ and $Z_{\rm min}$, defined in (34) and (33), respectively. In
particular, for constraint $\mathrm{C5}$, the lower and upper bounds on
altitude $z_{1}$ are derived by decomposing the baseline into its parallel and
perpendicular components, i.e.,
$b(\mathbf{q}_{1},\mathbf{q}_{2})=\sqrt{b_{\bot}^{2}(\mathbf{q}_{2})+b_{\parallel}^{2}(\mathbf{q}_{1},\mathbf{q}_{2})}$,
where its parallel components is given by:

$b_{\parallel}(\mathbf{q}_{1},\mathbf{q}_{2})=b(\mathbf{q}_{1},\mathbf{q}_{2})\sin\Big{(}\theta_{1}-\alpha(\mathbf{q}_{1},\mathbf{q}_{2})\Big{)}.$
(41)

Note that, based on constraint $\mathrm{C2}$, it can be shown that the
perpendicular baseline only depends on $\mathbf{q}_{2}$, see Appendix A in
[1]. Hence, the baseline can be rewritten as follows:

$b^{2}(\mathbf{q}_{1},\mathbf{q}_{2})=\big{(}r_{1}(\mathbf{q}_{1})-r_{2}(\mathbf{q}_{2})\cos(\theta_{1}-\theta_{2}(\mathbf{q}_{2}))\big{)}^{2}+b_{\bot}^{2}(\mathbf{q}_{2}).$
(42)

Here, we only focus on the case where
$b_{\mathrm{min}}>b_{\bot}(\mathbf{q}_{2})$, as constraint $\mathrm{C5}$ is
always satisfied otherwise. Therefore, based on (42), the bounds on $z_{1}$
are given by:

$\displaystyle\frac{z_{1}}{\cos(\theta_{1})}=r_{1}(\mathbf{q}_{1})\geq
r_{2}(\mathbf{q}_{2})\cos(\theta_{1}-\theta_{2}(\mathbf{q}_{2}))+\sqrt{b^{2}_{\mathrm{min}}-b^{2}_{\bot}},$
(43) $\displaystyle\frac{z_{1}}{\cos(\theta_{1})}=r_{1}(\mathbf{q}_{1})\leq
r_{2}(\mathbf{q}_{2})\cos(\theta_{1}-\theta_{2}(\mathbf{q}_{2}))-\sqrt{b^{2}_{\mathrm{min}}-b^{2}_{\bot}}.$
(44)

Lastly, constraint $\overline{\mathrm{C11}}$ is derived by expanding
constraint $\mathrm{C11}$ and rearranging its terms. Terms that are increasing
wrt $z_{1}$ are assigned to function ${a_{1}}$, whereas terms that are
decreasing wrt $z_{1}$ are grouped in function ${a_{2}}$.

## References

  * [1] M.-A. Lahmeri _et al._ , “UAV formation optimization for communication-assisted InSAR sensing,” _presented at IEEE Int. Conf. Commun. (ICC)_ , Jun 2024.
  * [2] M. Mozaffari _et al._ , “A tutorial on UAVs for wireless networks: Applications, challenges, and open problems,” _IEEE Commun. Surv. Tutor._ , vol. 21, no. 3, pp. 2334–2360, 2019.
  * [3] A. M. Haimovich, R. S. Blum, and L. J. Cimini, “MIMO radar with widely separated antennas,” _IEEE Signal Process. Mag._ , vol. 25, no. 1, pp. 116–129, 2008.
  * [4] A. Moreira _et al._ , “A tutorial on synthetic aperture radar,” _IEEE Geosci. Remote Sens. Mag._ , vol. 1, no. 1, pp. 6–43, 2013.
  * [5] R. Bamler and P. Hartl, “Synthetic aperture radar interferometry,” _Inverse Problems_ , vol. 14, no. 4, 1998.
  * [6] M. Martone _et al._ , “Coherence evaluation of TanDEM-X interferometric data,” _ISPRS J. Photogramm. Remote Sens._ , vol. 73, pp. 21–29, 2012.
  * [7] K. Meng _et al._ , “UAV-enabled integrated sensing and communication: Opportunities and challenges,” _IEEE Wirel. Commun._ , pp. 1–9, 2023.
  * [8] A. Khalili _et al._ , “Energy-aware resource allocation and trajectory design for UAV-enabled ISAC,” in _Proc. IEEE Global Commun. Conf._ , Dec 2023, pp. 4193–4198.
  * [9] K. Meng _et al._ , “UAV trajectory and beamforming optimization for integrated periodic sensing and communication,” _IEEE Wireless Commun. Letters_ , vol. 11, no. 6, pp. 1211–1215, 2022.
  * [10] Z. Lyu, G. Zhu, and J. Xu, “Joint maneuver and beamforming design for UAV-enabled integrated sensing and communication,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 4, pp. 2424–2440, 2023.
  * [11] Z. Sun _et al._ , “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, 2016.
  * [12] ——, “Trajectory optimization for maneuvering platform bistatic SAR with geosynchronous illuminator,” _IEEE Trans. Geosci. Remote Sens._ , vol. 62, pp. 1–15, 2024.
  * [13] ——, “Performance analysis and system implementation for energy-efficient passive UAV radar imaging system,” _IEEE Trans. Veh. Technol._ , vol. 72, no. 8, pp. 9938–9955, 2023.
  * [14] F. Xu _et al._ , “Heuristic path planning method for multistatic UAV-borne SAR imaging system,” _IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens._ , vol. 14, pp. 8522–8536, 2021.
  * [15] 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, 2022.
  * [16] M.-A. Lahmeri, W. R. Ghanem, C. Bonfert, and R. Schober, “Robust trajectory and resource optimization for communication-assisted UAV SAR sensing,” _IEEE Open J. Commun. Soc._ , pp. 1–1, 2024, early access.
  * [17] M.-A. Lahmeri _et al._ , “Trajectory and resource optimization for UAV synthetic aperture radar,” in _Proc. IEEE Global Commun. Conf._ , 2022, pp. 897–903.
  * [18] H. Ren _et al._ , “Swarm UAV SAR for 3-D imaging: System analysis and sensing matrix design,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–16, 2022.
  * [19] G. Krieger _et al._ , “TanDEM-X: A satellite formation for high-resolution SAR interferometry,” _IEEE Trans. Geosci. Remote Sens._ , vol. 45, no. 11, pp. 3317–3341, 2007.
  * [20] ——, “Interferometric synthetic aperture radar (SAR) missions employing formation flying,” _Proc. IEEE_ , vol. 98, no. 5, pp. 816–843, 2010.
  * [21] X. Chen _et al._ , “Performance of joint sensing-communication cooperative sensing UAV network,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 12, pp. 15 545–15 556, 2020.
  * [22] Z. Xiao and Y. Zeng, “Waveform design and performance analysis for full-duplex integrated sensing and communication,” _IEEE J. Sel. Areas Commun._ , vol. 40, no. 6, pp. 1823–1837, 2022.
  * [23] Z. Xu _et al._ , “Effect analysis and spectral weighting optimization of sidelobe reduction on SAR image understanding,” _IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens._ , vol. 12, no. 9, pp. 3434–3444, 2019.
  * [24] K. Meng, X. He, Q. Wu, and D. Li, “Multi-UAV collaborative sensing and communication: Joint task allocation and power optimization,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 6, pp. 4232–4246, 2023.
  * [25] D. L. Schuler _et al._ , “Terrain topography measurement using multipass polarimetric synthetic aperture radar data,” _Rad. Sci._ , vol. 35, no. 3, pp. 813–832, 2000.
  * [26] V. Mustieles-Perez _et al._ , “Towards UAV-based ultra-wideband multi-baseline SAR interferometry,” in _Proc. 20th European Rad. Conf_ , 2023, pp. 233–236.
  * [27] M. Zink _et al._ , “The TanDEM-X mission: Overview and status,” in _Proc. IEEE Int. Geosci. Remote Sens. Symp._ , 2007, pp. 3944–3947.
  * [28] I. Joughin, D. Winebrenner, and D. Percival, “Probability density functions for multilook polarimetric signatures,” _IEEE Trans. Geosci. Remote Sens._ , vol. 32, no. 3, pp. 562–574, 1994.
  * [29] S. Yueh _et al._ , “K-distribution and polarimetric terrain radar clutter,” _J. Electromagn. Waves Appl._ , vol. 3, no. 8, pp. 747–768, 1989\.
  * [30] M. Younis, “Digital beam-forming for high resolution wide swath real and synthetic aperture radar,” Ph.D. dissertation, Karlsruhe, Univ., 2004.
  * [31] 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.
  * [32] J. Schneider and S. Kirkpatrick, _Stochastic Optimization_. Springer Science & Business Media, 2007.
  * [33] D. Wang, D. Tan, and L. Liu, “Particle swarm optimization algorithm: an overview,” _Soft Computing_ , vol. 22, pp. 387–408, 2018.
  * [34] J. Kennedy and R. Eberhart, “Particle swarm optimization,” in _Proc. Int. Conf. Neural Netw._ , vol. 4, 1995, pp. 1942–1948.
  * [35] K. E. Parsopoulos _et al._ , “Particle swarm optimization method for constrained optimization problems,” _Intelligent Technologies: Theory and Application: New Trends in Intelligent Technologies_ , vol. 76, no. 1, pp. 214–220, 2002.
  * [36] T. M. Shami _et al._ , “Particle swarm optimization: A comprehensive survey,” _IEEE Access_ , vol. 10, pp. 10 031–10 061, 2022.
  * [37] J. Robinson and Y. Rahmat-Samii, “Particle swarm optimization in electromagnetics,” _IEEE Trans. Antennas Propag._ , vol. 52, no. 2, pp. 397–407, 2004.
  * [38] Q. He and L. Wang, “An effective co-evolutionary particle swarm optimization for constrained engineering design problems,” _Eng. Appl. Artif. Intell._ , vol. 20, no. 1, pp. 89–99, 2007.
  * [39] H. Garg, “A hybrid PSO-GA algorithm for constrained optimization problems,” _Appl. Math. Comput._ , vol. 274, pp. 292–305, 2016.
  * [40] Y. J. A. Zhang _et al._ , “Monotonic optimization in communication and networking systems,” _Found. Trends Netw._ , vol. 7, no. 1, pp. 1–75, Oct. 2013.
  * [41] H. Tuy, “Monotonic optimization: Problems and solution approaches,” _SIAM J. Optim._ , vol. 11, no. 2, pp. 464–494, Nov. 2000.
  * [42] S. Diamond and S. Boyd, “CVXPY: A python-embedded modeling language for convex optimization,” _J. Mach. Learn. Res._ , vol. 17, no. 83, pp. 1–5, 2016.
  * [43] Q. T. Dinh and M. Diehl, “Local convergence of sequential convex programming for nonconvex optimization,” _Recent Adv. Optim. Appl. Eng._ , pp. 93–102, 2010.
  * [44] G. Zhang, Q. Wu, M. Cui, and R. Zhang, “Securing UAV communications via joint trajectory and power control,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 2, pp. 1376–1389, 2019.
  * [45] J. C. Bezdek and R. J. Hathaway, “Some notes on alternating optimization,” in _Proc. AFSS Int. Conf. Fuzzy Sys._ Springer, 2002, pp. 288–300.