# Verification of some functional inequalities via polynomial optimization

Giovanni Fantuzzi  Department of Aeronautics, Imperial College London, London
SW7 2AZ (e-mail: giovanni.fantuzzi10@imperial.ac.uk).

###### Abstract

Motivated by the application of Lyapunov methods to partial differential
equations (PDEs), we study functional inequalities of the form
$f(I_{1}(u),\ldots,I_{k}(u))\geq 0$ where $f$ is a polynomial, $u$ is any
function satisfying prescribed constraints, and $I_{1}(u),\ldots,I_{k}(u)$ are
integral functionals whose integrands are polynomial in $u$, its derivatives,
and the integration variable. We show that such functional inequalities can be
strengthened into sufficient polynomial inequalities, which in principle can
be checked via semidefinite programming using standard techniques for
polynomial optimization. These sufficient conditions can be used also to
optimize functionals with affine dependence on tunable parameters whilst
ensuring their nonnegativity. Our approach relies on a measure-theoretic
lifting of the original functional inequality, which extends both a recent
moment relaxation strategy for PDE analysis and a dual approach to
inequalities for integral functionals.

###### keywords:

Functional inequalities, moment relaxations, polynomial optimization, Lyapunov
methods, partial differential equations, distributed parameter systems

††thanks: This work has been submitted to IFAC for possible publication. The
author was funded by an Imperial College Research Fellowship.

##  1 Introduction

Lyapunov methods are a powerful tool to analyse dynamical systems governed by
differential equations, wherein key system properties are inferred by
constructing one or more functions of the system’s state (and, possibly, of
time) subject to inequality constraints. For example, classical Lyapunov
functions that are positive definite and decay along trajectories prove the
stability of equilibrium states (Lyapunov, 1992; Datko, 1970). Other types of
Lyapunov-like functions, generically called auxiliary functions in this work,
can bound the effect of external disturbances (Willems, 1972; Dashkovskiy and
Mironchenko, 2013; Ahmadi et al., 2016), certify safety (Prajna et al., 2007;
Ahmadi et al., 2017; Miller et al., 2021), approximate reachable sets and
basins of attraction (Tan and Packard, 2006; Korda et al., 2013; Henrion and
Korda, 2014; Valmorbida and Anderson, 2017), estimate extreme behaviour
(Chernyshenko et al., 2014; Fantuzzi et al., 2016; Fantuzzi and Goluskin,
2020; Goluskin and Fantuzzi, 2019; Tobasco et al., 2018; Korda et al., 2021;
Goluskin, 2020), and solve optimal control problems (Lasserre et al., 2008;
Henrion et al., 2008; Korda et al., 2018).

For systems governed by ordinary differential equations (ODEs) with polynomial
right-hand side, polynomial auxiliary functions can often be constructed via
semidefinite programming. This is because the constraints on polynomial
auxiliary functions reduce to polynomial inequalities imposed on the entire
state space or on a subset thereof. While usually intractable, such
inequalities can be strengthened by insisting that all polynomials required to
be nonnegative admit a decomposition into sums of squares (SOS), which can be
searched for by solving semidefinite programs (SDPs) (Laurent, 2009; Parrilo,
2013; Lasserre, 2015). Given enough computational resources to solve these
SDPs, and barring issues with poor numerical conditioning, Lyapunov methods
can therefore be applied to polynomial ODEs in a systematic way.

The use of Lyapunov methods to study dynamical systems governed by partial
differential equations (PDEs), instead, remains a challenge for two reasons.
One is that there is considerable freedom in the structure of auxiliary
functions (more precisely, functionals) for PDEs, and it is not clear how to
generalise the notion of polynomial auxiliary functions for ODEs. The second
reason is that verifying if a candidate auxiliary functional satisfies all
relevant constraints—such as the Lyapunov conditions for nonlinear
stability—requires checking functional inequalities over infinite-dimensional
function spaces. This is usually hard to do, even with computer assistance.

Valmorbida et al. (2016, 2015) and Ahmadi et al. (2016, 2017, 2019)
demonstrated that SOS polynomials and SDPs can be used to verify the
nonnegativity of integral functionals with polynomial integrands in one or two
spatial dimensions. This is achieved without discretization of the PDE state,
but rather by requiring the polynomial integrand to be nonnegative pointwise
after augmenting it by terms that integrate to zero. This approach,
generalized to arbitrary spatial dimensions by Chernyavsky et al. (2021),
enables one to search for auxiliary functions in the form of integrals of
polynomials of the PDE state. A dual formulation of the same strategy, based
on so-called occupation measures, was independently introduced by Korda et al.
(2018) to relax optimal PDE control problems into SDPs.

Unfortunately, auxiliary functionals for PDEs that can be represented as
integrals of polynomials may not be sufficiently general. For example,
Lyapunov analysis of fluid flows (Goulart and Chernyshenko, 2012; Huang et
al., 2015; Goluskin and Fantuzzi, 2019; Fuentes et al., 2019) required
functionals of the form

$\mathcal{V}(u)=V(\mathcal{I}_{1}(u),\ldots,\mathcal{I}_{k}(u)),$ (1)

where $V$ is a $k$-variate polynomial and
$\mathcal{I}_{1}(u),\ldots,\mathcal{I}_{k}(u)$ are fixed integral functionals
of the fluid’s velocity $u$. Such auxiliary functionals were optimized by
projecting the governing PDEs onto a finite-dimensional ODE and estimating the
projection error.

Motivated by the desire to avoid application-specific estimates, here we study
the fundamental problem of verifying if functionals in the form (1) are
nonnegative for all $u$ satisfying prescribed constraints such as boundary
conditions. Leveraging the occupation measure framework of Korda et al.
(2018), we show that the broad class of functional inequalities obtained when
(1) is a _polynomial of integrals_ can be strengthened into sufficient
polynomial inequalities over finite-dimensional sets specified by linear
constraints and linear matrix inequalities. These purely algebraic problems
can be further strengthened into SDPs using standard techniques for polynomial
optimization, so one can attempt to verify hard functional inequalities using
algorithms for convex programming. In principle, therefore, one can implement
Lyapunov methods for PDEs with no need for rigorous ODE projections. Doing
this in practice, of course, requires overcoming well-known challenges related
to the numerical solution of SDPs, which may be too large or too ill-
conditioned for currently available general-purpose solvers. We do not address
such implementation issues here and leave them to future work. Nevertheless,
we demonstrate that our approach to verifying functional inequalities yields
sharp results on two relatively simple examples.

##  2 A class of functional inequalities

Let $\Omega\subset\mathbb{R}^{n}$ be an open bounded domain with smooth
boundary $\partial\Omega$. Given an integer $p$ with $1<p<\infty$, let
$W^{1,p}(\Omega;\mathbb{R}^{m})$ be the Sobolev space of $p$-integrable
functions $u:\Omega\to\mathbb{R}^{m}$ with $p$-integrable gradients. The
$m\times n$ Jacobian matrix $\nabla u$ has entries $(\nabla
u)_{ij}=\frac{\partial u_{i}}{\partial x_{j}}$.

Given an integer $d$, let $\mathcal{S}$ be a subset of the multi-index set

$\displaystyle\\!\mathcal{E}_{d,p}:=\\{(\alpha,\beta,\gamma)\in\mathbb{N}^{n}\times\mathbb{N}^{m}\times\mathbb{N}^{m\times
n}\\!:$ $\displaystyle\left|\alpha\right|\leq d,$
$\displaystyle\left|\beta\right|+\left|\gamma\right|\leq p\\}.$ (2)

Consider a functional $\mathcal{F}$ on $W^{1,p}(\Omega;\mathbb{R}^{m})$
defined as

$\mathcal{F}(u):=f\bigg{(}\left\\{\int_{\Omega}x^{\alpha}u(x)^{\beta}\nabla
u(x)^{\gamma}\,\,{\rm
d}{x}\right\\}_{(\alpha,\beta,\gamma)\in\mathcal{S}}\bigg{)},$ (3)

where $f$ is a function in $\left|\mathcal{S}\right|$ variables. Here and
throughout the paper, $x^{\alpha}=x_{1}^{\alpha_{1}}\cdots x_{n}^{\alpha_{n}}$
is the multivariate monomial with variables $x=(x_{1},\ldots,x_{n})$, exponent
$\alpha\in\mathbb{N}^{n}$, and degree
$\left|\alpha\right|=\alpha_{1}+\cdots+\alpha_{n}$. The multivariate monomials
$u^{\beta}$ and $\nabla u^{\gamma}$ are defined similarly. The bound
$\left|\beta\right|+\left|\gamma\right|\leq p$ in (2) ensures that the
integrals in (3) are well-defined.

We are interested in the following fundamental problems.

######  Problem 1 (Feasibility)

Check if $\mathcal{F}(u)$ in (3) is nonnegative for all functions $u$ in the
subset $\mathcal{U}$ of $W^{1,p}(\Omega;\mathbb{R}^{m})$ defined by the
differential and boundary conditions

$\displaystyle a\left({x},{u}({x}),\nabla{u}({x})\right)$
$\displaystyle=0\quad\text{on }\Omega,$ (4a) $\displaystyle
b\left({x},{u}({x})\right)$ $\displaystyle=0\quad\text{on }\partial\Omega,$
(4b)

where $a:\overline{\Omega}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times
n}\to\mathbb{R}$ and $b:\partial\Omega\times\mathbb{R}^{m}\to\mathbb{R}$ are
given continuous functions that do not grow faster than a degree-$p$
polynomial in the second and third arguments.

######  Problem 2 (Optimization)

Suppose $f$ in (3) depends affinely on a vector of parameters
$\lambda\in\mathbb{R}^{\ell}$ and let $\mathcal{U}$ be as in Problem 1. Given
a convex function $b:\mathbb{R}^{\ell}\to\mathbb{R}$, solve the convex
optimization problem

$\min_{\lambda\in\mathbb{R}^{\ell}}\,b(\lambda):\;\mathcal{F}(u)\geq
0\quad\forall u\in\mathcal{U}.$ (5)

Our approach to these problems is to derive sufficient but more tractable
conditions for the functional inequality $\mathcal{F}(u)\geq 0$. We do so in
three steps. First, we lift the functional inequality into an equivalent
inequality on so-called _occupation_ and _boundary measures_ (Korda et al.,
2018). We then strengthen the latter by enforcing it over a larger convex set
of measures, obtaining a measure-theoretic condition whose validity suffices
to conclude that $\mathcal{F}(u)\geq 0$ on $\mathcal{U}$. Finally, for
problems with polynomial data (cf. §4) we replace this measure-theoretic
condition with a sufficient constrained polynomial inequality that, at least
in principle, can be tested computationally. This polynomial inequality
depends affinely on $\lambda$ when so does $\mathcal{F}$, so it can be used to
strengthen the constraint in (5) and compute a feasible (and, sometimes, near-
optimal) $\lambda$.

###### Remark 3

Our results immediately extend to functions $f$ in (3) that depend also on
boundary integrals. Moreover, it is straightforward to consider function sets
$\mathcal{U}$ whose definitions, in addition to the equalities in (4), include
differential and boundary inequalities, equality or inequality constraints on
integrals over the domain or the boundary, and multiple constraints of each
type. Finally, only minor modifications are needed to consider domains
$\Omega$ whose boundary is only piecewise-smooth and Lipschitz.

###### Remark 4

If the function $f$ in (3) is linear, then $\mathcal{F}$ is an integral
functional and may be rewritten as

$\mathcal{F}(u)=\int_{\Omega}\phi({x},{u}(x),\nabla{u}(x))\,\,{\rm d}{x}$

for some polynomial $\phi$. In such cases, a duality argument (Fantuzzi, 2019)
shows that our sufficient conditions for the functional inequality in Problem
1 are equivalent to those proposed by Valmorbida et al. (2016, 2015), Ahmadi
et al. (2016, 2017, 2019), and Chernyavsky et al. (2021).

##  3 Lifting via occupation measures

This section describes the first two steps of the strategy outlined above.
Occupation and boundary measures are briefly reviewed in §3.1, while
sufficient measure-theoretic conditions for the inequality $\mathcal{F}(u)\geq
0$ on $\mathcal{U}$ are derived in §3.2. To lighten the notation, we will
often write $\langle f,\mu\rangle$ to indicate the integral of a function $f$
from a set $\Theta$ into $\mathbb{R}$ against a measure $\mu$ supported on
$\Theta$, i.e.,

$\langle f,\mu\rangle:=\int_{\Theta}f(\vartheta)\,{\rm d}\mu.$

###  3.1 Occupation and boundary measures

Given any function $u\in W^{1,p}(\Omega;\mathbb{R}^{m})$, let $\mu_{u}$ be the
pushforward of the Lebesgue measure on $\Omega$ through the map
$x\mapsto(x,u(x),\nabla u(x))$. By construction, $\mu_{u}$ is supported on the
product set

$D:=\overline{\Omega}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times n}$ (6)

and satisfies

$\int_{D}\eta(x,y,Z)\,{\rm d}\mu_{u}=\int_{\Omega}\eta(x,u(x),\nabla
u(x))\,\,{\rm d}{x}$ (7)

for all continuous functions $\eta:D\to\mathbb{R}$ that grow no faster than a
degree-$p$ polynomial in the last two arguments (so the right-hand side is
well-defined). Following Korda et al. (2018), we call $\mu_{u}$ the
_occupation measure_ generated by $u$. In complete analogy, the _boundary
measure_ $\nu_{u}$ generated by $u\in W^{1,p}(\Omega;\mathbb{R}^{m})$ is the
pushforward of the surface measure on $\partial\Omega$ through the map
$x\mapsto(x,u(x))$. By construction, $\nu_{u}$ is supported on the product set

$B:=\partial\Omega\times\mathbb{R}^{m}$ (8)

and satisfies

$\int_{B}\zeta(x,y)\,{\rm d}\nu_{u}=\int_{\partial\Omega}\zeta(x,u(x))\,{\rm
d}S$ (9)

for all continuous functions $\zeta:B\to\mathbb{R}$ that grow no faster than a
degree-$p$ polynomial in the second argument.

The occupation and boundary measures $\mu_{u}$ and $\nu_{u}$ are related via
identities implied by the divergence theorem. Let $\Phi$ be the linear space
of continuously differentiable $n$-valued functions
$\varphi:\overline{\Omega}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$ such that,
given any $u\in W^{1,p}(\Omega;\mathbb{R}^{m})$, the function
$x\mapsto\varphi(x,u(x))$ is differentiable using the chain rule and has an
integrable derivative. To every $\varphi\in\Phi$ we associate the “total
divergence” function

$\mathcal{D}\varphi(x,y,Z)=\sum_{i=1}^{n}\bigg{[}\frac{\partial\varphi_{i}}{\partial
x_{i}}(x,y)+\sum_{j=1}^{m}Z_{ji}\frac{\partial\varphi_{i}}{\partial
y_{j}}(x,y)\bigg{]}$

from $\overline{\Omega}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times n}$ into
$\mathbb{R}$. This definition ensures that $\mathcal{D}\varphi(x,u(x),\nabla
u(x))$ is the divergence of $\varphi(x,u(x))$ when the latter is calculated
using the chain rule. With this notation, it is straightforward to combine the
divergence theorem with identities (7) and (9) to conclude the following
result. (See also Theorem 1 in Korda et al., 2018.)

###### Lemma 5

Let $\mu_{u}$ and $\nu_{u}$ be the occupation and boundary measures generated
by a function $u\in W^{1,p}(\Omega;\mathbb{R}^{m})$, and let $\hat{n}(x)$ be
the outwards unit normal vector to $\partial\Omega$. Then,

$\langle\mathcal{D}\varphi,\mu_{u}\rangle=\langle\varphi\cdot\hat{n},\nu_{u}\rangle\quad\forall\varphi\in\Phi.$
(10)

###  3.2 Measure-theoretic strengthening of Problem 1

Consider Problem 1, which asks to verify if $\mathcal{F}(u)\geq 0$ for all
functions $u$ in the set $\mathcal{U}$ defined by the differential and
boundary constraints in (4). Since identity (7) implies that

$\mathcal{F}(u)=f(\\{\langle
x^{\alpha}y^{\beta}Z^{\gamma},\mu_{u}\rangle\\}_{(\alpha,\beta,\gamma)\in\mathcal{S}}),$

it is natural to extend $\mathcal{F}$ into a functional on a suitable space of
measures, and then reformulate the functional inequality $\mathcal{F}(u)\geq
0$ as a constraint on occupation measures.

To do this, let $\mathcal{M}_{+}(D)$ and $\mathcal{M}_{+}(B)$ be the cones of
Borel measures supported on the sets $D$ and $B$ defined in (6) and (8),
respectively. The functional

$\mathcal{G}(\mu,\nu)=f(\\{\langle
x^{\alpha}y^{\beta}Z^{\gamma},\mu\rangle\\}_{(\alpha,\beta,\gamma)\in\mathcal{S}})$
(11)

from $\mathcal{M}_{+}(D)\times\mathcal{M}_{+}(B)$ into $\mathbb{R}$ clearly
satisfies $\mathcal{G}(\mu_{u},\nu_{u})=\mathcal{F}(u)$ whenever
$(\mu_{u},\nu_{u})$ is a pair of occupation and boundary measures generated by
a function $u\in W^{1,p}(\Omega;\mathbb{R}^{m})$. Problem 1 may therefore be
restated as follows.

###### Problem 6

Check if the functional $\mathcal{G}(\mu,\nu)$ in (11) is nonnegative when
$(\mu,\nu)\in\mathcal{M}_{+}(D)\times\mathcal{M}_{+}(B)$ is a pair of
occupation and boundary measures generated by $u\in\mathcal{U}$.

###### Remark 7

The apparently redundant dependence of $\mathcal{G}$ on $\nu$ is useful in
what follows to exploit the relations between occupation and boundary measures
given by Lemma 5. Moreover, it enables us to easily extend our approach to
functionals $\mathcal{F}(u)$ that depend also on the surface integrals
$\int_{\partial\Omega}x^{\alpha}u(x)^{\beta}\,{\rm d}S$ with
$\left|\alpha\right|\leq d$ and $\left|\beta\right|\leq p$, which can be
rewritten as $\langle x^{\alpha}y^{\beta},\nu_{u}\rangle$ using (9).

Problem 6 is just as hard as verifying that $\mathcal{F}(u)\geq 0$ on the
function set $\mathcal{U}$. However, our measure-theoretic reformulation
suggests that to check this functional inequality it is enough to verify that
$\mathcal{G}$ is nonnegative on a subset
$\mathcal{A}\subseteq\mathcal{M}_{+}(D)\times\mathcal{M}_{+}(B)$ that includes
all occupation and boundary measures generated by $u\in\mathcal{U}$. For
suitably chosen $\mathcal{A}$, this sufficient measure-theoretic condition can
be more tractable that the original functional inequality.

One particularly convenient strategy is to define $\mathcal{A}$ using linear
constraints necessarily satisfied by occupation and boundary measures
generated by $u\in\mathcal{U}$. Lemma 5 already provides an uncountable number
of such constraints. Additional ones can be derived from identities (7) and
(9), and from the differential and boundary conditions (4a) and (4b) imposed
on $u\in\mathcal{U}$. Specifically, let $a$ and $b$ be the functions appearing
in (4). Further, let $\Psi_{a}$ and $\Psi_{b}$ be the sets of continuous
functions $\psi:D\to\mathbb{R}$ and $\rho:B\to\mathbb{R}$ such that the
functions $(\psi a)(x,y,Z)$ and $(\rho b)(x,y)$ grow no faster than a
degree-$p$ polynomial in $y$ and $Z$. Then,

###### Lemma 8

Every pair of occupation and boundary measures generated by $u\in\mathcal{U}$
belongs to the set $\mathcal{A}$ of pairs
$(\mu,\nu)\in\mathcal{M}_{+}(D)\times\mathcal{M}_{+}(B)$ satisfying

$\displaystyle\langle\mathcal{D}\varphi,\mu\rangle$
$\displaystyle=\langle\varphi\cdot\hat{n},\nu\rangle$
$\displaystyle\forall\varphi\in\Phi,$ (12a)
$\displaystyle\langle\eta,\mu\rangle$
$\displaystyle=\textstyle\int_{\Omega}\eta(x)\,\,{\rm d}{x}$
$\displaystyle\forall\eta\in C(\overline{\Omega}),$ (12b)
$\displaystyle\langle\zeta,\nu\rangle$
$\displaystyle=\textstyle\int_{\partial\Omega}\zeta(x)\,{\rm d}S$
$\displaystyle\forall\zeta\in C(\partial\Omega),$ (12c)
$\displaystyle\langle\psi a,\mu\rangle$ $\displaystyle=0$
$\displaystyle\forall\psi\in\Psi_{a},$ (12d) $\displaystyle\langle\rho
b,\nu\rangle$ $\displaystyle=0$ $\displaystyle\forall\rho\in\Psi_{b}.$ (12e)

If $\mu$ and $\nu$ are occupation and boundary measures generated by
$u\in\mathcal{U}$, condition (12a) is just a restatement of Lemma 5.
Conditions (12b) and (12c) follow directly from identities (7) and (9).
Condition (12d) holds because

$\langle\psi a,\mu\rangle=\int_{D}\psi(x,y,Z)a(x,y,Z)\,{\rm d}\mu\\\
=\int_{\Omega}\psi(x,u,\nabla u)\,a(x,u,\nabla u)\,\,{\rm d}{x}=0,$

where the last two equalities are a consequence of (7) and (4a), respectively.
Condition (12e) is proven analogously using (9) and (4b).  $\square$

Combining Lemma 8 with the preceding discussion yields the following result.

###### Proposition 9

Let $\mathcal{G}(\mu,\nu)$ be as in (11), and let $\mathcal{A}$ be the set of
measures $(\mu,\nu)\in\mathcal{M}_{+}(D)\times\mathcal{M}_{+}(B)$ satisfying
(12a)–(12e). If

$\mathcal{G}(\mu,\nu)\geq 0\quad\forall(\mu,\nu)\in\mathcal{A},$ (13)

then the functional $\mathcal{F}(u)$ in (3) is nonnegative for all functions
$u$ in the set $\mathcal{U}$ defined in Problem 1.

###### Remark 10

If the functions $u\in\mathcal{U}$ satisfy Dirichlet boundary conditions,
meaning that $b(x,u)=u-h(x)$ in (4b) for some fixed function $h$, then their
boundary measures $\nu_{u}$ coincide with the known measure $\nu_{h}$. In such
cases, one can replace $\nu$ with $\nu_{h}$ in Lemma 8 and Proposition 9.
Then, condition (12c) can be dropped, while the term
$\langle\varphi\cdot\hat{n},\nu\rangle$ in (12a) becomes the real number
$\langle\varphi\cdot\hat{n},\nu_{h}\rangle=\int_{\partial\Omega}\varphi(x,h(x))\cdot\hat{n}(x)\,{\rm
d}S$, which in principle can be computed explicitly for every choice of
$\varphi\in\Phi$.

###  3.3 Exploiting symmetries

Proposition 9 can be refined if the functional inequality in Problem 1 enjoys
symmetries. Specifically, if the domain $\Omega$, the integrals entering (3),
and the constraints in (4) are all invariant under a group of linear
transformations, then it suffices to verify (13) on the group-invariant subset
of $\mathcal{A}$.

To make these idea precise, let $\mathbb{G}$ be a group of orthonormal linear
transformations on $\mathbb{R}^{n}\times\mathbb{R}^{m}$, meaning a family
(finite, countable or uncountable) of pairs of matrices
$(A,B)\in\mathbb{R}^{n\times n}\times\mathbb{R}^{m\times m}$ such that:

  1. (G1)

$A^{\mathsf{T}}=A^{-1}$ and $B^{\mathsf{T}}=B^{-1}$;

  2. (G2)

If $(A,B)\in\mathbb{G}$, then $(A^{-1},B^{-1})\in\mathbb{G}$;

  3. (G3)

If $(A_{1},B_{1}),\,(A_{2},B_{2})\in\mathbb{G}$, then
$(A_{1}A_{2},B_{1}B_{2})\in\mathbb{G}$.

Given $(A,B)\in\mathbb{G}$, define linear operators $G_{A,B}$ and $H_{A,B}$ on
$\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times n}$ and
$\mathbb{R}^{n}\times\mathbb{R}^{m}$, respectively, via

$\displaystyle G_{A,B}(x,y,Z)$
$\displaystyle=\left(Ax,\,By,\,BZA^{\mathsf{T}}\right)$ $\displaystyle
H_{A,B}(x,y)$ $\displaystyle=\left(Ax,\,By\right).$

We say that Problem 1 is invariant under $\mathbb{G}$ if:

  1. (A1)

$x\in\Omega\iff Ax\in\Omega$ for all $(A,B)\in\mathbb{G}$;

  2. (A2)

$a\circ G_{A,B}=a$ for all $(A,B)\in\mathbb{G}$;

  3. (A3)

$b\circ H_{A,B}=b$ for all $(A,B)\in\mathbb{G}$;

  4. (A4)

$(Ax)^{\alpha}(By)^{\beta}(BZA^{\mathsf{T}})^{\gamma}=x^{\alpha}y^{\beta}Z^{\gamma}$
for all $(\alpha,\beta,\gamma)\in\mathcal{S}$ and all $(A,B)\in\mathbb{G}$.

In such cases, Proposition 9 can be replaced with

###### Proposition 11

Let $\mathcal{G}(\mu,\nu)$ be as in (11), and let $\mathcal{A}_{\mathbb{G}}$
be the set of measures
$(\mu,\nu)\in\mathcal{M}_{+}(D)\times\mathcal{M}_{+}(B)$ that satisfy
(12a)–(12e) as well as the invariance conditions

$\displaystyle\langle\eta\circ G_{A,B},\mu\rangle$
$\displaystyle=\langle\eta,\mu\rangle$ $\displaystyle\forall\eta\in
C(D),\,(A,B)\in\mathbb{G},$ (15a) $\displaystyle\langle\zeta\circ
H_{A,B},\nu\rangle$ $\displaystyle=\langle\zeta,\nu\rangle$
$\displaystyle\forall\zeta\in C(B),\,(A,B)\in\mathbb{G}.$ (15b)

If assumptions 1–4 hold and

$\mathcal{G}(\mu,\nu)\geq 0\quad\forall(\mu,\nu)\in\mathcal{A}_{\mathbb{G}},$
(16)

then the functional $\mathcal{F}(u)$ in (3) is nonnegative for all functions
$u$ in the set $\mathcal{U}$ defined in Problem 1.

This result is a direct consequence of the following lemma.

###### Lemma 12

Under conditions 1–4, for every $(\mu,\nu)\in\mathcal{A}$ there exists
$(\overline{\mu},\overline{\nu})\in\mathcal{A}_{\mathbb{G}}$ such that
$\mathcal{G}(\overline{\mu},\overline{\nu})=\mathcal{G}(\mu,\nu)$.

Let $\left|\mathbb{G}\right|$ be the cardinality of $\mathbb{G}$. Fix
$(\mu,\nu)\in\mathcal{A}$ and define the group-averaged measures

$\displaystyle\overline{\mu}$
$\displaystyle:=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\\!\\!G_{A,B}\sharp\mu,$
$\displaystyle\overline{\nu}$
$\displaystyle:=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\\!\\!H_{A,B}\sharp\nu,$

where $G_{A,B}\sharp\mu$ and $H_{A,B}\sharp\nu$ are the pushforwards of $\mu$
and $\nu$ by $G_{A,B}$ and $H_{A,B}$. The identity
$\mathcal{G}(\overline{\mu},\overline{\nu})=\mathcal{G}(\mu,\nu)$ follows upon
observing that

$\displaystyle\langle x^{\alpha}y^{\beta}Z^{\gamma},\overline{\mu}\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle
x^{\alpha}y^{\beta}Z^{\gamma},G_{A,B}\sharp\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle(Ax)^{\alpha}(By)^{\beta}(BZA^{\mathsf{T}})^{\gamma},\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle
x^{\alpha}y^{\beta}Z^{\gamma},\mu\rangle$ $\displaystyle=\langle
x^{\alpha}y^{\beta}Z^{\gamma},\mu\rangle$

for all $(\alpha,\beta,\gamma)\in\mathcal{S}$. The first equality is a direct
consequence of the linearity of integration against measures; the second one
follows from the definition of the pushforward measure; the third one exploits
assumption 4.

There remains to show that $\overline{\mu}$ and $\overline{\nu}$ satisfy
(15a), (15b), and (12a)–(12e). For (15a), observe that

$\displaystyle\langle\eta\circ G_{A,B},\overline{\mu}\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A^{\prime},B^{\prime})\in\mathbb{G}}\langle\eta\circ
G_{A,B},G_{A^{\prime},B^{\prime}}\sharp\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A^{\prime},B^{\prime})\in\mathbb{G}}\langle\eta\circ
G_{A,B}\circ G_{A^{\prime},B^{\prime}},\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A^{\prime\prime},B^{\prime\prime})\in\mathbb{G}}\langle\eta\circ
G_{A^{\prime\prime},B^{\prime\prime}},\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A^{\prime\prime},B^{\prime\prime})\in\mathbb{G}}\langle\eta,G_{A^{\prime\prime},B^{\prime\prime}}\sharp\mu\rangle$
$\displaystyle=\langle\eta,\overline{\mu}\rangle.$

The first and last equalities exploit the linearity of integration against
measures; the second and fourth ones follow from the definition of the
pushforward measure; the third one combines conditions 2 and 3 with the
definition of $G_{A,B}$. Analogous steps show that $\overline{\nu}$ satisfies
(15b).

Next, we show that $\overline{\mu}$ and $\overline{\nu}$ satisfy (12a). Fix an
arbitrary function $\varphi\in\Phi$ and, for every $(A,B)\in\mathbb{G}$,
define another function $\varphi_{A,B}\in\Phi$ via

$\varphi_{A,B}(x,y):=A^{\mathsf{T}}\varphi(Ax,By).$

Thanks to condition 1, it is an exercise in vector calculus to check that
$\mathcal{D}\varphi_{A,B}=\mathcal{D}\varphi\circ G_{A,B}$. Consequently,

$\displaystyle\langle\mathcal{D}\varphi,\overline{\mu}\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle\mathcal{D}\varphi,G_{A,B}\sharp\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle\mathcal{D}\varphi\circ
G_{A,B},\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle\mathcal{D}\varphi_{A,B},\mu\rangle$
$\displaystyle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle\varphi_{A,B}\cdot\hat{n},\nu\rangle.$

The last equality is true because $(\mu,\nu)\in\mathcal{A}$ by assumption, so
$\mu$ and $\nu$ satisfy (12a). Now, assumption 1 implies that the unit normal
vector $\hat{n}(x)$ to $\partial\Omega$ satisfies $\hat{n}(Ax)=A\hat{n}(x)$,
so $\varphi_{A,B}(x,y)\cdot\hat{n}(x)=(\varphi\cdot\hat{n})\circ H_{A,B}$.
Consequently,

$\langle\mathcal{D}\varphi,\overline{\mu}\rangle=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle(\varphi\cdot\hat{n})\circ
H_{A,B},\nu\rangle\\\
=\frac{1}{\left|\mathbb{G}\right|}\sum_{(A,B)\in\mathbb{G}}\langle\varphi\cdot\hat{n},H_{A,B}\sharp\nu\rangle=\langle\varphi\cdot
n,\overline{\nu}\rangle.$

Since $\varphi\in\Phi$ was arbitrary, the pair
$(\overline{\mu},\overline{\nu})$ satisfies (12a) as claimed. The remaining
conditions, (12b)–(12e), are verified with similar steps (omitted for
brevity).  $\square$

##  4 Reduction to an SDP

The measure-theoretic inequality (13) and its symmetry-exploiting analogue
(16), although simpler than the original functional inequality
$\mathcal{F}(u)\geq 0$ from Problem 1, remain difficult to check in general.
Often, however, they can be further strengthened into finite-dimensional
sufficient conditions that can be tested computationally. One case in which
this can be done systematically is that of problems with polynomial data,
meaning that:

  1. (P1)

The function $f$ in (3) is polynomial;

  2. (P2)

The functions $a$ and $b$ in (4) are polynomials;

  3. (P3)

The domain $\Omega$ is a basic semialgebraic set, i.e., it is defined by
polynomial equalities and inequalities.

We shall also assume that:

  1. (P4)

There exist polynomials $g$ and $h$ such that

$\displaystyle\Omega$ $\displaystyle=\\{x\in\mathbb{R}^{n}:g(x)\geq 0\\},$
(18a) $\displaystyle\partial\Omega$
$\displaystyle=\\{x\in\mathbb{R}^{n}:g(x)=0\\}.$ (18b)

  2. (P5)

The unit normal vector $\hat{n}(x)=\nabla g(x)\left|\nabla g(x)\right|^{-1}$
to $\partial\Omega$ is polynomial.

The last two assumptions are made purely for simplicity. It is straightforward
to extend our discussion to domains with semialgebraic definitions more
complicated than (18). When $\hat{n}(x)$ is not polynomial, instead, one can
redefine the boundary measures in §3.1 to be the pushforward of the weighted
surface measure $\left|\nabla g(x)\right|^{-1}{\rm d}S$. Then, all results in
§3 hold with $\hat{n}(x)$ replaced by $\nabla g(x)$ and with an extra factor
of $\left|\nabla g(x)\right|^{-1}$ in the surface integrals in (9) and (12c).

###  4.1 A sufficient polynomial inequality

Under assumptions 1–5, the measure-theoretic inequalities (13) and (16) can be
strengthened into polynomial inequalities over sets defined by linear
equations and linear matrix inequalities (LMIs). To do this, observe that
$\mathcal{G}$ depends on the moments of $\mu$ with exponents
$\mathcal{S}\subset\mathcal{E}_{d,p}$, which have degree
$\left|\alpha\right|\leq d$ and $\left|\gamma\right|+\left|p\right|\leq p$
(cf. §2). Recall also that the polynomials $a(x,y,Z)$ and $b(x,y)$ in (4)
cannot have degree larger than $p$ in $y$ and $Z$, and there is no loss of
generality in assuming that their degree in $x$ is $d$ or less. We may
similarly assume that the polynomial $g$ in (18) has degree no larger than
$d$.

Now, fix any integer $\omega\geq d/2$ (called the relaxation order) and define
the moment vectors

$\displaystyle\xi$ $\displaystyle:=\\{\langle
x^{\alpha}y^{\beta}Z^{\gamma},\mu\rangle:\;\left|\alpha\right|\leq
2\omega,\;\left|\beta\right|+\left|\gamma\right|\leq p\\},$
$\displaystyle\theta$ $\displaystyle:=\\{\langle
x^{\alpha}y^{\beta},\nu\rangle:\;\left|\alpha\right|\leq
2\omega+\deg(\hat{n}),\;\left|\beta\right|\leq p\\}.$

The functional $\mathcal{G}$ in (11) is clearly a polynomial of $\xi$. If
$2\omega=d$ we have exactly $\mathcal{G}(\mu,\nu)=f(\xi)$, where $f$ is the
polynomial in (3). With a slight abuse of notation we write the same when
$2\omega>d$, which amounts to viewing $f$ as a higher-dimensional polynomial.
Thus, to prove inequalities (13) and (16) it suffices to show that $f$ is
nonnegative for all moment vectors $\xi$. In turn, this is true if $f(\xi)\geq
0$ for all $\xi$ in a set defined by LMIs and linear equations necessarily
satisfied by all moment vectors.

The LMIs involve the well-known _moment matrices_ and _localizing matrices_
(Laurent, 2009; Lasserre, 2015). Precisely, given any polynomial
$h(x)=\sum_{\left|\alpha\right|\leq d_{h}}h_{\alpha}x^{\alpha}$ of degree
$d_{h}\leq d$, let $M_{h}(\xi)$ be a matrix with rows and columns indexed by
the multi-indices $(\alpha,\beta,\gamma)$ with
$\left|\alpha\right|\leq\omega-\lceil\frac{1}{2}d_{h}\rceil$ and
$\left|\beta\right|+\left|\gamma\right|\leq\lfloor\frac{1}{2}p\rfloor$, and
whose entries are

$[M_{h}(\xi)]_{(\alpha,\beta,\gamma),(\alpha^{\prime},\beta^{\prime},\gamma^{\prime})}=\sum_{\left|\alpha^{\prime\prime}\right|\leq
d_{h}}h_{\alpha^{\prime\prime}}\xi_{\alpha+\alpha^{\prime}+\alpha^{\prime\prime},\beta+\beta^{\prime},\gamma+\gamma^{\prime}}$

Let $M_{h}(\theta)$ be defined in an analogous way, except that its rows and
columns are indexed by multi-indices $(\alpha,\beta)$ with
$\left|\alpha\right|\leq\omega-\lceil\frac{1}{2}d_{h}\rceil$ and
$\left|\beta\right|\leq\lfloor\frac{1}{2}p\rfloor$. Then, under assumptions
1–5, the moment vectors $\xi$ and $\theta$ are known to satisfy the block-
diagonal LMI

$M(\xi,\theta):=\begin{pmatrix}M_{1}(\xi)\\\ &M_{g}(\xi)\\\
&&M_{1}(\theta)\end{pmatrix}\succeq 0.$ (20)

The linear equality constraints on $\xi$ and $\theta$, instead, are obtained
from (12a)–(12e) and (in the presence of symmetries) from (15a) and (15b).
This is done by restricting the test functions in these conditions to be
monomials of suitably bounded degree. Let $\Phi_{\omega}\subset\Phi$ be a
monomial basis for the space of $n$-valued polynomials $\varphi(x,y)$ of
degree $2\omega$ in $x$ and degree $p$ in $y$:

$\Phi_{\omega}:=\big{\\{}x^{\alpha}y^{\beta}\hat{e}_{k}:\;\left|\alpha\right|\leq
2\omega,\,\left|\beta\right|\leq p,\,k=1,\ldots,n\big{\\}},$

where $\hat{e}_{k}$ is the $n$-dimensional unit vector in the $k$-th
coordinate direction. Further, define the monomial sets

$\displaystyle\Psi_{a,\omega}$
$\displaystyle:=\big{\\{}x^{\alpha}y^{\beta}Z^{\gamma}:\,\left|\beta\right|+\left|\gamma\right|+\deg_{y,Z}(a)\leq
p,$ $\displaystyle\hskip 96.0pt\left|\alpha\right|+\deg_{x}(a)\leq
2\omega\big{\\}},$ $\displaystyle\Psi_{b,\omega}$
$\displaystyle:=\big{\\{}x^{\alpha}y^{\beta}:\,\left|\beta\right|+\deg_{y}(b)\leq
p,$ $\displaystyle\hskip 52.0pt\left|\alpha\right|+\deg_{x}(b)\leq
2\omega+\deg(\hat{n})\big{\\}}.$

Then, conditions (12a)–(12e) imply

$\displaystyle\langle\mathcal{D}\varphi,\mu\rangle$
$\displaystyle=\langle\varphi\cdot\hat{n},\nu\rangle$
$\displaystyle\forall\varphi\in\Phi_{\omega},$ (22a) $\displaystyle\langle
x^{\alpha},\mu\rangle$ $\displaystyle=\textstyle\int_{\Omega}x^{\alpha}\,{\rm
d}{x}$ $\displaystyle\forall\alpha:\left|\alpha\right|\leq 2\omega,$ (22b)
$\displaystyle\langle x^{\alpha},\nu\rangle$
$\displaystyle=\textstyle\int_{\partial\Omega}x^{\alpha}\,{\rm d}S$
$\displaystyle\forall\alpha:\left|\alpha\right|\leq 2\omega+\deg(\hat{n}),$
(22c) $\displaystyle\langle\psi a,\mu\rangle$ $\displaystyle=0$
$\displaystyle\forall\psi\in\Psi_{a,\omega},$ (22d) $\displaystyle\langle\rho
b,\nu\rangle$ $\displaystyle=0$ $\displaystyle\forall\rho\in\Psi_{b,\omega}.$
(22e)

Thanks to the imposed degree bounds, these identities are linear equations for
the moment vectors $\xi$ and $\theta$. For problems with symmetries, further
linear equations can be obtained by setting
$\eta=x^{\alpha}y^{\beta}Z^{\gamma}$ with $\left|\alpha\right|\leq 2\omega$
and $\left|\beta\right|+\left|\gamma\right|\leq p$ in (15a), and by taking
$\zeta=x^{\alpha}y^{\beta}$ with $\left|\alpha\right|\leq 2\omega$ and
$\left|\beta\right|\leq p$ in and (15b). For notational convenience, we
combine all equality constraints into the equation

$e(\xi,\theta):=A\xi+B\theta+c=0.$ (23)

The matrices $A$ and $B$ and the vector $c$ can be constructed explicitly in
applications.

The next result is immediate.

###### Proposition 13

Inequalities (13) and (16), and hence the functional inequality in Problem 1,
hold if

$f(\xi)\geq 0\quad\forall(\xi,\theta)\text{ s.t.~{}\eqref{e:lmi}
and~{}\eqref{e:equalities}}.$ (24)

###### Remark 14

The sufficient condition for the functional inequality in Problem 1 provided
by Proposition 13 becomes less restrictive when the relaxation order $\omega$
(hence, the size of the vectors $\xi$ and $\theta$) is increased. However, in
general we do not expect inequality (24) to become a necessary condition for
Problem 1 as $\omega\to\infty$.

###  4.2 Implementation via semidefinite programming

Inequality (24) can be verified computationally via semidefinite programming
using (a variation of) known techniques based on sum-of-squares (SOS)
polynomials and SOS polynomial matrices (Scherer and Hol, 2006).

We begin by stating some preliminary results about positive semidefinite (PSD)
matrices. For any integers $m$ and $n$, let $\mathbb{S}^{n}$ be the space of
$n\times n$ symmetric matrices and let $\mathbb{S}^{n}_{+}$ be the PSD cone.
Given $P\in\mathbb{S}^{n}$ and $Q\in\mathbb{S}^{mn}$ partitioned into $m^{2}$
blocks $Q_{ij}$ of size $n\times n$, define

$\mathcal{T}(P,Q):=\begin{pmatrix}\operatorname{trace}(PQ_{11})&\cdots&\operatorname{trace}(PQ_{1n})\\\
\vdots&\ddots&\vdots\\\
\operatorname{trace}(PQ_{n1})&\cdots&\operatorname{trace}(PQ_{nn})\\\
\end{pmatrix}\in\mathbb{S}^{m}.$

Note that $\mathcal{T}(P,Q)$ is simply the trace inner product of $P$ and $Q$
if $m=1$. When $m=n^{k}$ for some integer $k\geq 1$, we define the matrices
$\mathcal{T}_{i}(P,Q)\in\mathbb{S}^{n^{k-i}}$ for $i=1,\ldots,k$ via

$\displaystyle\mathcal{T}_{1}(P,Q)$ $\displaystyle=\mathcal{T}(P,Q),$
$\displaystyle\mathcal{T}_{i}(P,Q)$
$\displaystyle=\mathcal{T}[P,\,\mathcal{T}_{i-1}(P,Q)],\quad i=2,\ldots,k.$

Observe that $\mathcal{T}_{i}(P,Q)$ depends linearly on $Q$ for every $i$, and
that $\mathcal{T}_{k}(P,Q)$ is a scalar.

###### Lemma 15

If $P\in\mathbb{S}^{n}_{+}$ and $Q\in\mathbb{S}^{n^{k}}_{+}$, then
$\mathcal{T}_{i}(P,Q)$ is PSD for each $i=1,\ldots,k$. In particular,
$\mathcal{T}_{k}(P,Q)\geq 0$.

The case $i=1$ was proven by Scherer and Hol (2006, §2.2). The general case
follows by induction.  $\square$

Lemma 15 immediately implies the next result, which may be seen as an
application of the so-called S-procedure.

###### Proposition 16

Let $\rho_{\omega}$ be the number of rows in the LMI (20) and
$\sigma_{\omega}$ be the number of constraints in (23). The polynomial
inequality (24), hence the functional inequality in Problem 1, holds if there
exist a polynomial vector $q(\xi,\theta)$ with $\sigma_{\omega}$ entries and
PSD polynomial matrices $\\{Q_{k}(\xi,\theta)\\}_{k=1,\ldots,K}$, each of size
$(\rho_{\omega})^{k}\times(\rho_{\omega})^{k}$, such that

$f(\xi)\geq q(\xi,\theta)\cdot
e(\xi,\theta)+\sum_{k=1}^{K}\mathcal{T}_{k}(M(\xi,\theta),Q_{k}(\xi,\theta)).$
(25)

It is well known (see, e.g., Lasserre, 2015; Laurent, 2009; Scherer and Hol,
2006; Parrilo, 2013) that if inequality (25) and the positive semidefinitess
constraints on the polynomial matrices $Q_{k}$ are strengthened into SOS
constraints, then the search for these matrices and for the polynomial vector
$q$ can be reformulated as an SDP once the number $K$ and the polynomial
degrees are fixed. The same is true when $f$ depends affinely on parameters
$\lambda$, as in Problem 2, because in such cases (25) is jointly affine in
$\lambda$, $q$, and the matrices $Q_{k}$. Consequently, one can use
semidefinite programming to (attempt to) verify the functional inequality in
Problem 1 and to compute feasible $\lambda$ for the convex optimization
problem (5).

###### Remark 17

For $K=1$, inequality (25) reduces to the sufficient conditions for polynomial
inequalities considered by Scherer and Hol (2006). This is enough when the set
of $(\xi,\theta)$ satisfying (20) and (23) is compact and satisfies the so-
called Archimedean condition (Scherer and Hol, 2006, Theorem 1), but not
otherwise: as demonstrated in §5.2, in general it is necessary to use $K>1$.
Beyond compactness, we are not aware of conditions ensuring that (25) is
feasible if $K$ and the degrees of $q$ and each $Q_{k}$ are sufficiently
large.

##  5 Examples

We conclude by showing that Proposition 16 works well on two examples with
known analytical solutions.

###  5.1 A Jensen–Poincaré inequality

Consider the problem of finding the minimum $\lambda$ such that

$\lambda\int_{-1}^{1}u_{x}^{2}\,{\rm d}{x}-\left(\int_{-1}^{1}u\,{\rm
d}{x}\right)^{2}\geq 0$ (26)

for all functions $u:(-1,1)\to\mathbb{R}$ satisfying the boundary conditions
$u(\pm 1)=0$. This problem fits the framework of §2 with $\Omega=(-1,1)$,
$n=m=1$, $p=2$, $a(x,u,\nabla u)\equiv 0$ and $b(x,u)=u$. Analytical solution
using the calculus of variation shows that (26) holds if
$\lambda\geq\frac{2}{3}$. Here, we verify this via the sufficient condition in
Proposition 16. We use $\omega=1$ and take a semialgebraic definition of
$\Omega$ as in (18) with $g(x)=1-x^{2}$.

Using the vector $\xi=\\{\xi_{\alpha,\beta,\gamma}:\,\alpha\leq
2,\beta+\gamma\leq 2\\}$, the polynomial inequality corresponding to (26) is

$f(\xi)=\lambda\xi_{0,0,2}-\xi_{0,1,0}^{2}\geq 0.$ (27)

To determine the set over which this inequality is required to hold, recall
from Remark 10 that we can fix $\nu=\nu_{0}$ to be the boundary measure of the
zero function, so $\langle x^{\alpha},\nu\rangle=1^{\alpha}-(-1)^{\alpha}$ for
all $\alpha$ and $\langle x^{\alpha}y^{\beta},\nu\rangle=0$ for all $\alpha$
and all $\beta>0$. With this simplification, identities (22a)–(22e) give
equality constraints on the entries of $\xi$ alone. Specifically, (22a) yields

$\displaystyle 2-\xi_{0,0,0}$ $\displaystyle=0,$ $\displaystyle\xi_{1,0,0}$
$\displaystyle=0,$ $\displaystyle\tfrac{2}{3}-\xi_{2,0,0}$ $\displaystyle=0.$
(28)

Condition (22b), instead, gives

$\displaystyle\xi_{0,0,1}$ $\displaystyle=0,$ $\displaystyle\xi_{0,1,1}$
$\displaystyle=0,$ (29a) $\displaystyle 2-\xi_{0,0,0}$ $\displaystyle=0,$
$\displaystyle\xi_{0,1,0}+\xi_{1,0,1}$ $\displaystyle=0,$ (29b)
$\displaystyle\xi_{0,2,0}+2\xi_{1,1,1}$ $\displaystyle=0,$
$\displaystyle\xi_{1,0,0}$ $\displaystyle=0,$ (29c) $\displaystyle
2\xi_{1,1,0}+\xi_{2,0,1}$ $\displaystyle=0,$
$\displaystyle\xi_{1,2,0}+\xi_{2,1,1}$ $\displaystyle=0.$ (29d)

Identities (22c)–(22e) give no constraints because $\nu$ is fixed and $a=0$.
However, since our example is invariant under the group
$\mathbb{G}=\\{(1,1),\,(-1,1)\\}$ generating the variable transformation
$(x,u,u_{x})\mapsto(-x,u,-u_{x})$ we can use (15a) to obtain a final set of
constraints:

$\displaystyle\xi_{0,0,1}$ $\displaystyle=0,$ $\displaystyle\xi_{0,1,1}$
$\displaystyle=0,$ $\displaystyle\xi_{1,0,0}$ $\displaystyle=0,$
$\displaystyle\xi_{1,0,2}$ $\displaystyle=0,$ (30a) $\displaystyle\xi_{1,1,0}$
$\displaystyle=0,$ $\displaystyle\xi_{1,2,0}$ $\displaystyle=0,$
$\displaystyle\xi_{2,0,1}$ $\displaystyle=0,$ $\displaystyle\xi_{2,1,1}$
$\displaystyle=0.$ (30b)

The LMI (20) for this example does not depend on $\theta$ and it is block-
diagonal with two blocks,

$M(\xi,\theta)=\begin{pmatrix}M_{1}(\xi)&0\\\ 0&M_{g}(\xi)\end{pmatrix}\succeq
0.$ (31)

Using (28)–(30) to eliminate variables that have a known value or can be
expressed in terms of other ones, we find

$\displaystyle
M_{1}(\xi)={\small\begin{pmatrix}2&\xi_{0,1,0}&-\xi_{0,1,0}&0&0&0\\\
\xi_{0,1,0}&-2\xi_{1,1,1}&\xi_{1,1,1}&0&0&0\\\
-\xi_{0,1,0}&\xi_{1,1,1}&\xi_{2,0,2}&0&0&0\\\
0&0&0&\xi_{0,0,2}&-\xi_{0,1,0}&\xi_{1,1,1}\\\
0&0&0&-\xi_{0,1,0}&2/3&\xi_{2,1,0}\\\
0&0&0&\xi_{1,1,1}&\xi_{2,1,0}&\xi_{2,2,0}\end{pmatrix}}$ and $\displaystyle
M_{g}(\xi)={\small\begin{pmatrix}4/3&\xi_{0,1,0}-\xi_{2,1,0}&0\\\
\xi_{0,1,0}-\xi_{2,1,0}&-\xi_{1,1,1}-\xi_{2,2,0}&0\\\
0&0&\xi_{0,0,2}-\xi_{2,0,2}\end{pmatrix}}.$

Next, we apply Proposition 16 to verify (27) for all vectors $\xi$ satisfying
(28)–(31). Take $q(\xi,\theta)=0$, $K=1$, and

$Q_{1}(\xi,\theta)=\kappa^{2}vv^{\mathsf{T}}+\tfrac{3}{2}w(\xi)w(\xi)^{\mathsf{T}},$

where $\kappa$ is an arbitrary constant and

$\displaystyle v^{\mathsf{T}}$
$\displaystyle=\left(0,\,0,\,0,\,1,\,0,\,0,\,0,\,0,\,0\right)$ $\displaystyle
w(\xi)^{\mathsf{T}}$
$\displaystyle=\left(0,\,0,\,0,\,2/3,\,\xi_{0,1,0},\,0,\,0,\,0,\,0\right).$

With these choices, inequality (25) reduces to

$\left(\lambda-\tfrac{2}{3}-\kappa^{2}\right)\xi_{0,0,2}\geq 0$

and it holds if $\lambda=\frac{2}{3}+\kappa^{2}$. Since $\kappa$ is arbitrary,
we conclude that the original functional inequality (26) is satisfied for all
$\lambda\geq\frac{2}{3}$, as claimed at the start of the example.

###  5.2 An inhomogeneous inequality

Let $\|v\|_{2}=\int_{-1}^{1}v^{2}\,{\rm d}{x}$ be the usual $L^{2}$ norm of a
square-integrable function $v:(-1,1)\to\mathbb{R}$. Given any positive integer
$d$, consider finding the minimum scalar $\lambda$ such that

$\lambda\|u_{x}\|_{2}^{2d}\|u\|_{2}^{2d}-2\|u\|_{2}^{2d}+1\geq 0$ (32)

for all $u:(-1,1)\to\mathbb{R}$ satisfying $u(\pm 1)=0$. It follows from the
Poincaré inequality $\|u_{x}\|_{2}^{2}\geq(\pi/2)^{2}\|u\|_{2}^{2}$ that (32)
holds if $\lambda\geq(2/\pi)^{2d}=:\lambda_{d}^{*}$.

Table 1 lists numerical upper bounds on $\lambda_{d}^{*}$ obtained for $d=1$,
$2$, $3$ and $4$, computed using MOSEK (MOSEK ApS, 2019) to solve the SDPs
corresponding to the SOS strengthening of the sufficient condition (25) from
Proposition 16. This condition reads

$f(\xi):=\lambda\xi_{0,0,2}^{d}\xi_{0,2,0}^{d}-2\xi_{0,2,0}^{d}+1\geq 0$ (33)

for all vectors $\xi=\\{\xi_{\alpha,\beta,\gamma}:\alpha\leq
2\omega,\,\beta+\gamma\leq 2\\}$ satisfying suitable equality constraints and
an LMI $M(\xi)\succeq 0$, which can be explicitly constructed but are too long
to be reported. As in our previous example, there is no dependence on $\theta$
by virtue of the Dirichlet boundary conditions on $u$. The tabulated bounds
were obtained for increasing values of the relaxation order $\omega$ with
$K=2$ and polynomial matrices $Q_{1}$ and $Q_{2}$ of degree $2d-2$, after
enforcing all equality constraints on $\xi$ explicitly. Computations with
$K=1$ give identical results if $d$ is even but are infeasible if $d$ is odd,
showing that setting $K>1$ is necessary in general.

The bounds in Table 1 improve as the relaxation order $\omega$ is raised, and
although they appear conservative they provably converge to the analytical
optimum $\lambda^{*}_{d}$ as $\omega\to\infty$. Loosely speaking, this is
because analysis by Chernyavsky et al. (2021, Theorem 4.3), combined with a
duality argument (Fantuzzi, 2019), guarantees that our sufficient conditions
capture with arbitrary accuracy the Poincaré inequality needed to verify (32).
Precisely, for every $\kappa>(\pi/2)^{2}$ there exists a sufficiently large
$\omega$ and a PSD matrix $P$ of the same size as the matrix $M(\xi)$ for that
$\omega$ value such that, for all $\xi$ that satisfy the equality constraints,

$\xi_{0,0,2}=\kappa\xi_{0,2,0}+\operatorname{trace}[M(\xi)P].$

Taking the $d$-th power of this identity, multiplying both sides by
$\xi_{0,2,0}^{d}$, and rearranging gives

$\xi_{0,0,2}^{d}\xi_{0,2,0}^{d}-\kappa^{d}\xi_{0,2,0}^{2d}=\sum_{j=1}^{d}\textstyle\binom{d}{j}\displaystyle\kappa^{d-j}\xi_{0,2,0}^{2d-j}s_{j}(\xi),$
(34)

where
$s_{j}(\xi)=\left(\operatorname{trace}[M(\xi)P]\right)^{j}=(\xi_{0,0,2}-\kappa\xi_{0,2,0})^{j}$.
Now, it can be shown that $\xi_{0,2,0}$ is a diagonal entry of $M(\xi)$, so
there exists a PSD diagonal matrix $D$ such that
$\operatorname{trace}(DM(\xi))=\xi_{0,2,0}$. Using this matrix, it is possible
to rewrite the right-hand side of (34) as $\mathcal{T}_{2}[M(\xi),R(\xi)]$
with

$R(\xi):=\sum_{\begin{subarray}{c}j=1\\\ j\text{
odd}\end{subarray}}^{d}\xi_{0,2,0}^{2d-j-1}s_{j-1}(\xi)\left(D\otimes
P\right)\\\\[-12.91663pt] +\sum_{\begin{subarray}{c}j=2\\\ j\text{
even}\end{subarray}}^{d}\xi_{0,2,0}^{2d-j}s_{j-2}(\xi)\left(P\otimes
P\right),$

which can be easily shown to be an SOS matrix. Then, if $f(\xi)$ is as in (33)
and we fix $K=2$, $q(\xi)=0$, $Q_{1}(\xi)=0$, and $Q_{2}(\xi)=\lambda R(\xi)$,
the sufficient condition (25) from Proposition 16 for our example reduces to

$\displaystyle 0\leq
f(\xi)-\mathcal{T}_{2}[M(\xi),Q_{2}(\xi)]=\lambda\kappa^{d}\xi_{0,2,0}^{2d}-2\xi_{0,2,0}^{d}+1.$

This inequality, hence the original functional inequality (32), holds if and
only if $\lambda\geq\kappa^{-d}$. This construction works for every
$\kappa>(\frac{\pi}{2})^{2}$ provided that $\omega$ is large enough, so the
functional inequality (32) can be verified with the sufficient condition in
Proposition 16 for $\lambda$ as close to
$(\frac{2}{\pi})^{2d}=\lambda_{d}^{*}$ as desired.

Table 1: Upper bounds on the minimum $\lambda$ for which inequality (32) holds. The analytical minimum $\lambda_{d}^{*}$ is reported for comparison. $\omega$ | $d=1$ | $d=2$ | $d=3$ | $d=4$
---|---|---|---|---
1 | 4.000000 | 16.00000 | 64.00000 | 256.0000
2 | 1.333333 | 1.777775 | 2.370369 | 3.160493
3 | 0.854470 | 0.730117 | 0.623862 | 0.533069
4 | 0.677588 | 0.459121 | 0.311073 | 0.210776
5 | 0.590550 | 0.348746 | 0.205849 | 0.121443
$\lambda_{d}^{*}$ | 0.405285 | 0.164256 | 0.066570 | 0.026980

This work benefited from discussions with David Goluskin, Jason Bramburger,
Alexandr Chernyavsky and Sergei Chernyshenko.

## References

  * Ahmadi et al. (2019) Ahmadi, M., Valmorbida, G., Gayme, D., and Papachristodoulou, A. (2019).  A framework for input-output analysis of wall-bounded shear flows.  _J. Fluid Mech._ , 873, 742–785.
  * Ahmadi et al. (2016) Ahmadi, M., Valmorbida, G., and Papachristodoulou, A. (2016).  Dissipation inequalities for the analysis of a class of PDEs.  _Automatica_ , 66, 163–171.
  * Ahmadi et al. (2017) Ahmadi, M., Valmorbida, G., and Papachristodoulou, A. (2017).  Safety verification for distributed parameter systems using barrier functionals.  _Syst. Control Lett._ , 108, 33–39.
  * Chernyavsky et al. (2021) Chernyavsky, A., Bramburger, J.J., Fantuzzi, G., and Goluskin, D. (2021).  Convex relaxations of integral variational problems: pointwise dual relaxation and sum-of-squares optimization.  arXiv:2110.03079.
  * Chernyshenko et al. (2014) Chernyshenko, S.I., Goulart, P.J., Huang, D., and Papachristodoulou, A. (2014).  Polynomial sum of squares in fluid dynamics: a review with a look ahead.  _Philos. Trans. Roy. Soc. A_ , 372(2020), 20130350.
  * Dashkovskiy and Mironchenko (2013) Dashkovskiy, S. and Mironchenko, A. (2013).  Input-to-state stability of infinite-dimensional control systems.  _Math. Control Signals Systems_ , 25(1), 1–35.
  * Datko (1970) Datko, R. (1970).  Extending a theorem of A. M. Liapunov to Hilbert Space.  _J. Math. Anal. Appl._ , 32(3), 610–616.
  * Fantuzzi (2019) Fantuzzi, G. (2019).  Duality of convex relaxations for constrained variational problems.  arXiv:1906.12136.
  * Fantuzzi and Goluskin (2020) Fantuzzi, G. and Goluskin, D. (2020).  Bounding extreme events in nonlinear dynamics using convex optimization.  _SIAM J. Appl. Dyn. Syst._ , 19(3), 1823–1864.
  * Fantuzzi et al. (2016) Fantuzzi, G., Goluskin, D., Huang, D., and Chernyshenko, S.I. (2016).  Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization.  _SIAM J. Appl. Dyn. Syst._ , 15(4), 1962–1988.
  * Fuentes et al. (2019) Fuentes, F., Goluskin, D., and Chernyshenko, S. (2019).  Global stability of fluid flows despite transient growth of energy.  arXiv:1911.09079.
  * Goluskin (2020) Goluskin, D. (2020).  Bounding extrema over global attractors using polynomial optimisation.  _Nonlinearity_ , 33(9), 4878–4899.
  * Goluskin and Fantuzzi (2019) Goluskin, D. and Fantuzzi, G. (2019).  Bounds on mean energy in the Kuramoto–Sivashinsky equation computed using semidefinite programming.  _Nonlinearity_ , 32(5), 1705–1730.
  * Goulart and Chernyshenko (2012) Goulart, P.J. and Chernyshenko, S.I. (2012).  Global stability analysis of fluid flows using sum-of-squares.  _Phys. D_ , 241(6), 692–704.
  * Henrion and Korda (2014) Henrion, D. and Korda, M. (2014).  Convex computation of the region of attraction of polynomial control systems.  _IEEE Trans. Automat. Control_ , 59(2), 297–312.
  * Henrion et al. (2008) Henrion, D., Lasserre, J.B., and Savorgnan, C. (2008).  Nonlinear optimal control synthesis via occupation measures.  In _Proc. IEEE Conf. Decis. Control_ , 4749–4754.
  * Huang et al. (2015) Huang, D., Chernyshenko, S.I., Goulart, P.J., Lasagna, D., Tutty, O.R., and Fuentes, F. (2015).  Sum-of-squares of polynomials approach to nonlinear stability of fluid flows: an example of application.  _Proc. Roy. Soc. A_ , 471(2183), 20150622.
  * Korda et al. (2013) Korda, M., Henrion, D., and Jones, C.N. (2013).  Inner approximations of the region of attraction for polynomial dynamical systems.  _IFAC Proc. Volumes_ , 43(23), 534–539.
  * Korda et al. (2018) Korda, M., Henrion, D., and Lasserre, J.B. (2018).  Moments and convex optimization for analysis and control of nonlinear partial differential equations.  arXiv:1804.07565.
  * Korda et al. (2021) Korda, M., Henrion, D., and Mezić, I. (2021).  Convex computation of extremal invariant measures of nonlinear dynamical systems and Markov processes.  _J. Nonlinear Sci._ , 31, 14(1–26).
  * Lasserre (2015) Lasserre, J.B. (2015).  _An introduction to polynomial and semi-algebraic optimization_.  Cambridge University Press.
  * Lasserre et al. (2008) Lasserre, J.B., Henrion, D., Prieur, C., and Trélat, E. (2008).  Nonlinear optimal control via occupation measures and LMI-relaxations.  _SIAM J. Control Optim._ , 47(4), 1643–1666.
  * Laurent (2009) Laurent, M. (2009).  Sums of squares, moment matrices and optimization over polynomials.  In M. Putinar and S. Sullivant (eds.), _Emerging Applications of Algebraic Geometry_ , The IMA Volumes in Mathematics and its Applications, vol 149, 157–270. Springer New York.
  * Lyapunov (1992) Lyapunov, A.M. (1992).  Stability of motion: General problem.  _Internat. J. Control_ , 55(3), 539–589.  Translated by A. T. Fuller from a French translation of Lyapunov’s 1892 dissertation.
  * Miller et al. (2021) Miller, J., Henrion, D., and Sznaier, M. (2021).  Peak Estimation Recovery and Safety Analysis.  _IEEE Control Syst. Lett._ , 5(6), 1982–1987.
  * MOSEK ApS (2019) MOSEK ApS (2019).  _The MOSEK optimization toolbox for MATLAB manual. Version 9.2._ URL https://docs.mosek.com/9.2/toolbox/index.html.
  * Parrilo (2013) Parrilo, P.A. (2013).  Polynomial optimization , sums of squares and applications.  In G. Blekherman, P.A. Parrilo, and R.R. Thomas (eds.), _Semidefinite optimization and convex algebraic geometry_ , chapter 3, 47–157. SIAM.
  * Prajna et al. (2007) Prajna, S., Jadbabaie, A., and Pappas, G.J. (2007).  A framework for worst-case and stochastic safety verification using barrier certificates.  _IEEE Trans. Automat. Control_ , 52(8), 1415–1428.
  * Scherer and Hol (2006) Scherer, C.W. and Hol, C.W.J. (2006).  Matrix sum-of-squares relaxations for robust semi-definite programs.  _Math. Program. B_ , 107, 189–211.
  * Tan and Packard (2006) Tan, W. and Packard, A. (2006).  Stability region analysis using sum of squares programming.  In _Proc. Amer. Control Conf._ , 2297–2302.
  * Tobasco et al. (2018) Tobasco, I., Goluskin, D., and Doering, C.R. (2018).  Optimal bounds and extremal trajectories for time averages in nonlinear dynamical systems.  _Phys. Lett. A_ , 382(6), 382–386.
  * Valmorbida et al. (2015) Valmorbida, G., Ahmadi, M., and Papachristodoulou, A. (2015).  Convex solutions to integral inequalities in two-dimensional domains.  In _Proc. 54 th IEEE Conf. Decis. Control_, 7268–7273.
  * Valmorbida et al. (2016) Valmorbida, G., Ahmadi, M., and Papachristodoulou, A. (2016).  Stability analysis for a class of partial differential equations via semidefinite programming.  _IEEE Trans. Automat. Control_ , 61(6), 1649–1654.
  * Valmorbida and Anderson (2017) Valmorbida, G. and Anderson, J. (2017).  Region of attraction estimation using invariant sets and rational Lyapunov functions.  _Automatica_ , 75, 37–45.
  * Willems (1972) Willems, J.C. (1972).  Dissipative dynamical systems part I: General theory.  _Arch. Ration. Mech. Anal._ , 45(5), 311–314.


# On the Optimal Beamwidth of UAV-Assisted Networks Operating at Millimeter
Waves

Manishika Rawat1, Marco Giordani1, Brejesh Lall2, Abdelaali Chaoub3, and
Michele Zorzi1
E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>1Department of Information Engineering, University of
Padova, Padova, Italy
2Department of Electrical Engineering, Indian Institute of Technology Delhi,
New Delhi, India
3National Institute of Posts and Telecommunications (INPT), Rabat, Morocco

###### Abstract

The millimeter-wave (mm-wave) bands enable very large antenna arrays that can
generate narrow beams for beamforming and spatial multiplexing. However,
directionality introduces beam misalignment and leads to reduced energy
efficiency. Thus, employing the narrowest possible beam in a cell may not
necessarily imply maximum coverage. The objective of this work is to determine
the optimal sector beamwidth for a cellular architecture served by an unmanned
aerial vehicle (UAV) acting as a base station (BS). The users in a cell are
assumed to be distributed according to a Poisson Point Process (PPP) with a
given user density. We consider hybrid beamforming at the UAV, such that
multiple concurrent beams serve all the sectors simultaneously. An
optimization problem is formulated to maximize the sum rate over a given area
while limiting the total power available to each sector. We observe that, for
a given transmit power, the optimal sector beamwidth increases as the user
density in a cell decreases, and varies based on the height of the UAV. Thus,
we provide guidelines towards the optimal beamforming configurations for users
in rural areas.

###### Index Terms:

UAV-BS, millimeter-wave, optimal sector beamwidth, rural connectivity.

This paper has been accepted for publication at the IEEE Wireless
Communications and Networking Conference Workshops (WCNC WKSHPS), 2023.
Copyright may change without notice.

##  I Introduction

According to the World Bank, 43 % of the world population lives in rural areas
[1]. However, these regions remain mostly unserved while the world prepares to
roll out the fifth generation (5G) of mobile networks. According to a report
published by the International Telecommunication Union in 2021, the share of
Internet users in urban areas is twice the number in rural areas [2]. The
primary cause behind this is the lack of communication infrastructure in
remote areas. The next generation of wireless networks (6G) emerges as a
solution to this challenge [3]: for example, 6G is focusing on non-terrestrial
networks (NTN) using Unmanned Aerial Vehicles (UAVs), High Altitude Platforms
(HAPs), and satellites to promote ubiquitous and high-capacity global
connectivity [4]. Notably, these modules can serve as aerial base stations
(BS) or to assist the terrestrial BS in providing on-demand, cost-effective
coverage in unserved and poorly served areas [5].

UAVs, in particular, have been proposed to bridge the digital divide and
provide on-demand networks for applications such as disaster management,
medical camps, network exploration, and surveillance [6, 7]. Deploying a UAV
as a BS (UAV-BS) is quick and affordable compared to a terrestrial network
infrastructure, and when operating in the millimeter-wave (mm-wave) bands can
promote cost-effective ubiquitous coverage, high throughput, and low latency
even in rural areas [8].

However, how to deploy UAV-BSs is a challenging design issue, and has been
discussed in detail in [9, 10]. At the same time, allocating the (limited)
resources to the ground users is critical. Several resource allocation
problems have been formalized in the literature to meet different quality of
service (QoS) requirements such as coverage, fairness, and energy efficiency
[11, 12, 13, 14, 15, 16, 17].

Figure 1: Geometry of the scenario where a UAV-BS serves a cellular area of
radius $R$.

An energy-efficient UAV communication system was proposed in [11] by
optimizing the UAV trajectory and jointly considering the communication
throughput and the energy consumption. In [12], the authors introduced an
approximate beam pattern and provided a solution for the UAV deployment and
beam gain allocation to maximize the capacity over a given area. In [13],
multi-UAV communication and non-orthogonal multiple access (NOMA) have been
combined for the purpose of constructing high-capacity Internet of Things
(IoT) uplink transmission systems. The channel assignment, the uplink transmit
power of IoT nodes, and the flying heights of UAVs have been jointly optimized
to maximize the system capacity. In [14], NOMA transmission with UAV-BS is
implemented to provide coverage over a dense user region. A beam scanning
approach has been proposed to identify the optimal area within the user
region, and hence maximize the achievable sum rate. In [15], the locations of
transceivers in the downlink and uplink were modeled using a Poisson Point
Process (PPP) or Poisson Cluster Process to derive closed-form expressions of
the coverage probability. In [16], an optimal resource allocation problem has
been investigated for downlink coverage. It considered concurrent transmission
to all sectors, and an asymptotically-optimal solution has been proposed to
solve a mixed integer non-linear programming problem. The authors in [17] have
proposed an intelligent UAV-BS placement and power allocation algorithm to
maximize the sum rate in a region.

The optimal beamwidth for UAV-assisted multi-user systems has been studied in
[18] considering the main lobe of the directional antenna serving a cell. The
ground terminals were partitioned into disjoint clusters, which were
sequentially served by the UAV, and the joint UAV height and beamwidth have
been investigated. The altitude, beamwidth, and location of the UAV and the
bandwidth allocated to each user were jointly optimized in [19] for uplink UAV
communication to minimize the sum uplink power. An algorithm was proposed to
obtain a suboptimal solution by assigning different bandwidths to ground
terminals. In [20], the UAV location and antenna beamwidth were jointly
optimized for quasi-stationary and mobile UAVs. The impact of the beamwidth on
UAV location and trajectory has been investigated for minimizing serving time
and increasing throughput. All of these works, however, assume that the
vertical main lobe of the UAV antenna covers all users.

In this paper we determine the optimal sector beamwidth for a UAV-assisted
cellular network implementing hybrid beamforming, where users are distributed
according to a PPP with a given user density. An optimization problem is
developed for efficient resource allocation to maximize the sum rate in a
sector while ensuring fairness. We observe that the optimal sector beamwidth
depends on the user density, the height of the UAV-BS, and the propagation
scenario. For a sub-urban scenario, with cell radius of $100$ m and a UAV-BS
deployed at a height of $100$ m, the optimal beamwidth decreases from
$10^{\circ}$ to $5^{\circ}$ as the user density increases from $0.0005$ to
$0.002$ UEs/$\text{\,}{\mathrm{m}}^{2}$. On the other hand, for a fixed user
density of $0.0005$ UEs/$\text{\,}{\mathrm{m}}^{2}$, the optimal beamwidth
initially decreases from $12^{\circ}$ to $10^{\circ}$ and then increases to
$15^{\circ}$ as the UAV height increases from $10$ m to $200$ m. We further
observe that the number of sectors required to optimally serve a given number
of users in an urban region is much larger than that in a rural region for the
same QoS requirements.

The paper has been organized as follows. In Sec. II we present our system
model; in Sec. III we describe the optimization problem to maximize the sum
rate in a cell within given constraints, and present the algorithm we used to
compute the sum rate as a function of the sector beamwidth; in Sec. IV we show
our simulation results; conclusions are given in Sec. V.

Notations: For a random variable $X$, $X\sim\text{Poisson}(x)$ denotes that
$X$ is Poisson distributed with rate parameter $x$,
$X\sim\text{U}(x_{1},x_{2})$ indicates that $X$ is uniformly distributed in
the range $(x_{1},x_{2}]$, and $X\sim\text{Tr}(x_{1},x_{2})$ denotes that $X$
has triangular distribution in the range $(x_{1},x_{2}]$ [21].

##  II System Model

The system model involves a UAV-BS deployed at a height $h$ and serving a
cellular area of radius $R$, as shown in Fig. 1. The cell is divided into $S$
sectors, each of beamwidth $\theta$. The UAV is mounted with a uniform planar
array (UPA) antenna of $N$ elements such that $\theta\approx 2/N$ radians
[22]. The users, represented by crosses in Fig. 1, are distributed in the cell
according to a PPP with density $\lambda$. The location of the $k^{th}$ user
in a cell is given by $(d_{k},\phi_{k})$, where $d_{k}$ is the horizontal
distance between the $k^{th}$ user and the UAV-BS and $\phi_{k}$ is the phase
of the $k^{th}$ user measured counterclockwise. Here,
$\phi\sim\text{U}(0,2\pi)$ and $d\sim\text{Tr}(0,R)$. The average number of
users in a cell is $M=\pi R^{2}\lambda$. We assume fixed transmit power
$P_{t}$ at the UAV, and orthogonal frequency division multiplexing (OFDM) to
serve multiple users on a single beam. Therefore, the total bandwidth $B$ is
split into $N_{c}$ subcarriers to serve multiple users in a sector. We operate
at mm-wave frequency in an effort to maximize the communication capacity, and
assume hybrid beamforming at the UAV such that all sectors are served by
concurrent beams [23].

The data rate $r_{n,k}$ for the $n^{th}$ subcarrier allocated to the $k^{th}$
user is given by

$\displaystyle r_{n,k}=(B/N_{c})\log_{2}(1+P_{n,k}\gamma_{n,k}),$ (1)

where $\gamma_{n,k}$ is given by [24], i.e.,

$\displaystyle\gamma_{n,k}(l_{k})=$
$\displaystyle\frac{|\mathbbm{h}_{LoS}^{n,k}|^{2}PL_{LoS}^{n,k}(l_{k})P_{r}(l_{k})+|\mathbbm{h}_{NLoS}^{n,k}|^{2}PL_{NLoS}^{n,k}(l_{k})(1-P_{r}(l_{k}))}{N_{0}B/N_{c}GG_{r}}.$
(2)

In Eq. (2), $l_{k}=\sqrt{h^{2}+d_{k}^{2}}$ is the distance between the UAV and
the $k^{th}$ user, and $\mathbbm{h}_{LoS}^{n,k}$ and
$\mathbbm{h}_{NLoS}^{n,k}$ are the path gains for the line-of-sight (LoS) and
non-line-of-sight (NLoS) paths between the UAV and the $k^{th}$ user.
Specifically, $k\in\\{1,2,\dots,K_{s}\\}$, where $K_{s}$ represents the number
of users in the $s^{th}$ sector, $s\in\\{1,2,\dots,S\\}$, and
$n\in\\{1,2,\dots,N_{c}\\}$. $N_{0}$ is the noise power spectral density, and
$G$ and $G_{r}$ represent the beamforming gains of the transmitting and
receiving antennas, respectively. $P_{r}(l_{k})$ is the LoS probability
between the UAV and the $k^{th}$ user, and is given by

$\displaystyle
P_{r}(l_{k})=\frac{1}{1+e^{\alpha_{1}\psi^{3}+\alpha_{2}\psi^{2}+\alpha_{3}\psi+\alpha_{4}}},$
(3)

where $\psi=\sin^{-1}\left(h/l_{k}\right)$ is the elevation angle of the UAV
with respect to the $k^{th}$ user, and $\alpha_{1}$, $\alpha_{2}$,
$\alpha_{3}$, and $\alpha_{4}$ are parameters of the LoS probability
distribution defined for sub-urban, urban, dense-urban, and high-rise-urban
environments in [25, Table 2]. $PL_{LoS}$ and $PL_{NLoS}$ represent the path
losses for LoS and NLoS links; for mm-wave communication [26], they are
expressed in dB as

$\displaystyle PL_{LoS}(l_{k})=61.4+20\log_{10}(l_{k})+\mathcal{N}(0,33.64);$
(4) $\displaystyle
PL_{NLoS}(l_{k})=72.0+29.2\log_{10}(l_{k})+\mathcal{N}(0,75.69).$ (5)

where, $\mathcal{N}(\mu,\sigma^{2})$ represents a normal random variable with
mean $\mu$ and variance $\sigma^{2}$. In the following section we will develop
an optimization problem to maximize the sum rate in a cell while limiting the
total power available in a sector.

##  III Sum Rate Maximization Problem

The sum rate maximization problem for a sector can be expressed as:

$\displaystyle\\!\max_{\pi_{n,k},P_{n,k}}$
$\displaystyle\sum_{k=1}^{K_{s}}\sum_{n=1}^{N_{c}}r_{n,k}\pi_{n,k}$ (6a) s.t.
$\displaystyle C_{1}:$ $\displaystyle\sum_{k=1}^{K_{s}}\pi_{n,k}\leq
1\quad\forall n,$ (6b) $\displaystyle C_{2}:$
$\displaystyle\sum_{n=1}^{N_{c}}r_{n,k}\pi_{n,k}\geq R_{0}\quad\forall k,$
(6c) $\displaystyle C_{3}:$
$\displaystyle\sum_{k=1}^{K_{s}}\sum_{n=1}^{N_{c}}P_{n,k}\leq\frac{P_{t}\theta}{360},$
(6d) $\displaystyle C_{4}:$ $\displaystyle\pi_{n,k}\in\\{0,1\\}\quad\forall
n,k,$ (6e) $\displaystyle C_{5}:$ $\displaystyle P_{n,k}\geq 0\quad\forall
n,k.$ (6f)

In (6a), we introduce a binary variable $\pi_{n,k}$ to ensure that at least
one subcarrier is allocated to one user in a sector in the beam serving time.
Here, $\theta$ is the beamwidth of a sector, and $R_{0}$ is the bit rate.
$P_{n,k}$ is the power transmitted over the $n^{th}$ subcarrier allocated to
the $k^{th}$ user. Constraints $C_{1}$ and $C_{4}$ ensure that the same
subcarrier is not allocated to different users, and thus $\pi_{n,k}$ can take
only binary values such that the sum across the users is less than or equal to
$1$. $C_{2}$ ensures a minimum data rate to each user. As the cell is divided
into sectors, the total power available to a sector would be proportional to
the sector beamwidth. Therefore, $C_{3}$ limits the total power available to a
sector. This limits the maximum data rate available to users as the sector
beamwidth decreases or the number of sectors in a cell increases. $C_{5}$
ensures that the power allocated to each user is non-negative.

The optimization problem in (6a) is a mixed integer non-linear programming
problem (MINLP). Due to the non-convexity, the global optimal solution cannot
be achieved. The time taken to run this module increases proportionally with
the number of users and sectors. However, we can simplify the problem using
the mixed integer (MI) property of $\pi_{n,k}$. The simplified optimization
problem can be expressed as

$\displaystyle\\!\max_{\pi_{n,k},P_{n,k}}$
$\displaystyle\sum_{k=1}^{K_{s}}\sum_{n=1}^{N_{c}}r_{n,k}$ (7a) s.t.
$\displaystyle C_{1},C_{4},C_{5},$ $\displaystyle C_{6}:$
$\displaystyle\sum_{n=1}^{N_{c}}r_{n,k}\geq R_{0}\quad\forall k,$ (7b)
$\displaystyle C_{7}:$
$\displaystyle\sum_{k=1}^{K_{s}}\sum_{n=1}^{N_{c}}P_{n,k}\leq\pi_{n,k}\frac{P_{t}\theta}{360},$
(7c) $\displaystyle C_{8}:$ $\displaystyle
r_{n,k}\leq\pi_{n,k}R_{max}\quad\forall n,k.$ (7d)

initialize $N_{c}$, $P_{t}$, $B$, $R_{0}$, $G_{0}$, $G_{r}$, $N_{0}$, $R$,
$h$, $\lambda$, $\theta\in div(360)$,

for  _each value of $\theta$_ do

$S=360/\theta$ ;

/*number of sectors

_ _ for  _$\eta=1,2,\dots,\mathbb{N}$ ;_

__ /*MC simulations

do

$\chi\sim\text{Poisson}(\pi R^{2}\lambda)$, $\phi\sim\text{U}(0,2\pi)$,
$d\sim\text{Tr}(0,R)$;

/*user distribution parameters

for  _$s=1,2,\dots,S$_ do

Based on the value of $\phi_{k}$ and $\theta$, determine the users that lie in
the $s^{th}$ sector;

Compute $g_{n,k}(N_{c},K_{s})$ for each sector;

Run optimization module given in (7) to allocate $N_{c}$ subcarriers to
$K_{s}$ users;

Compute $r_{n,k}(N_{c},K_{s})$ for each sector;

$Op_{s}(s)=\sum_{n=1}^{N_{c}}\sum_{k=1}^{K_{s}}r_{n,k}(N_{c},K_{s})$ ;

/*sum rate per sector

end for

$Op_{\eta}(\eta)=\sum_{s=1}^{S}Op_{s}$

$Ap_{\eta}(\eta)=Op_{\eta}(\eta)/\chi$

end for

$Op=\sum_{\eta=1}^{\mathbb{N}}Op_{\eta}/\mathbb{N}$;

/*sum rate per cell

$Ap=\sum_{\eta=1}^{\mathbb{N}}Ap_{\eta}/\mathbb{N}$;

/*average rate

end for

Algorithm 1 To compute the sum rate per cell as a function of $\theta$ for a
given value of $\lambda$

We have omitted $\pi_{n,k}$ from (7a) and (7b) to reduce the redundancy in
problem formulation. The condition that $r_{n,k}=0$ when $\pi_{n,k}=0$ is
enforced by $C_{8}$. $C_{5}$ and $C_{7}$ ensure that $P_{n,k}$ drops to zero
when $\pi_{n,k}=0$. Thus, the algorithm needs to run only for non-zero entries
in $\pi_{n,k}$, which reduces the problem complexity. $R_{max}$ in $C_{8}$ is
the maximum data rate that can be achieved. It can be used to limit the number
of iterations of the solver to save the processing time.

The resource allocation problem in (7) can be explained by taking the example
of Fig. 1. Here, $S=8$ and we set $N_{c}=4$ so that both $\pi_{n,k}$ and
$P_{n,k}$ for the second sector measured counterclockwise in Fig. 1 are of
size $4\times 2$. Thus, the four subcarriers in a sector can be allocated to
two users. Because of $C_{1}$, one of the possible solutions would be
$\pi_{n,k}=\begin{pmatrix}1&0\\\ 0&1\\\ 0&1\\\ 0&1\end{pmatrix}$. Accordingly,
$P_{n,k}=\begin{pmatrix}P_{1,2}&0\\\ 0&P_{2,3}\\\ 0&P_{3,3}\\\
0&P_{4,3}\end{pmatrix}$ such that $\sum_{n=1}^{4}\sum_{k=2}^{3}P_{n,k}\leq
P_{t}/8$.

Algorithm 1 specifies the steps to compute the sum rate per cell and average
rate per user as a function of $\theta$ for a given user density. First, we
generate users distributed as a PPP in a cell of radius $R$ for a given value
of $\theta$ and $\lambda$. The users in the same sector are then categorized
for analysis. The optimization function maximizes the sum rate per sector,
which is added to produce the sum rate per cell. The average rate per user is
computed by dividing the sum rate per cell by the number of users. Here,
$\mathbb{N}$ represents the number of Monte Carlo (MC) simulations and
$div(x)$ represents the divisors of $x$.

##  IV Results and Discussions

In this section, we evaluate the impact of user density and UAV height on the
optimal sector beamwidth of a cell. We work with $B=1$ GHz, $P_{t}=10$ W,
$N_{c}=30$, $R_{0}=1$ Gbps, $R_{max}=50$ Gbps, $f_{c}=28$ GHz, and
$N_{0}=-174$ dBm/Hz. The channel is assumed to be Rician with a distribution
parameter of $8$ dB [27]. We use the solving constraint integer program (scip)
in the General Algebraic Modeling System (GAMS), a high-level modeling system
for mathematical optimization to solve the MINLP problem in (6a). We consider
(i) a rural/sub-urban scenario with a lower user density
$\lambda=\\{0.0005,0.0008,0.001,0.002\\}$ UEs/$\text{\,}{\mathrm{m}}^{2}$ in a
cell of radius $R=100$ m (Sec. IV-A), and (ii) an urban scenario with a higher
user density $\lambda=\\{0.05,0.08,0.1\\}$ UEs/$\text{\,}{\mathrm{m}}^{2}$ in
a cell of radius $R=10$ m (Sec. IV-B). The impact of the UAV height is
explored in Sec. IV-C. We assume $G=N$ and $G_{r}=1$ for the analysis [22].
The results have been obtained for $500$ MC simulations.

(a) Sum rate per cell.

(b) Average rate per user.

(c) Jain’s fairness index.

Figure 2: Sum rate per cell, average rate per user, and Jain’s fairness index
in a rural/sub-urban scenario with $R=100$ m and $h=100$ m, as a function of
$\theta$, for different user densities.

###  IV-A Rural/Sub-Urban Scenario

Fig. 2(a) plots the sum rate obtained by solving the optimization problem in
(7) as a function of $\theta$. The result is obtained for $R=100$ m, $h=100$
m, and different values of $\lambda$. First, we observe that as the user
density increases, the number of users in the cell increases, thus increasing
the sum rate per cell, which validates the accuracy of our results. Initially,
as $\theta$ increases the sum rate also increases given that the power
allocated to each sector increases proportionally. However, after a certain
value of $\theta$, the sum rate decreases exponentially. This can be
attributed to the lower spatial reuse as $\theta$ increases, i.e., as the
number of sectors decreases. Therefore, we conclude that the smallest
beamwidth does not necessarily ensure the maximum sum rate. The sector
beamwidth at which the sum rate is maximized is represented by $\theta_{opt}$.
According to Fig. 2(a), $\theta_{opt}=10^{\circ}$ for $\lambda=0.0005$. As
$\lambda$ increases from $0.0008$ to $0.002$, $\theta_{opt}$ decreases from
$9^{\circ}$ to $5^{\circ}$, respectively. This implies that in a rural
environment, it is desirable to operate through wider beams as the density of
users decreases.

Fig. 2(b) plots the average rate per user as a function of $\theta$ for
different values of $\lambda$. The use of the mm-wave technology ensures that
the rate is always above 1 Gbps, and in line with the requirements of most
5G/6G applications. As expected, as the number of users in the cell increases,
the average rate per user decreases. Moreover, as $\theta$ increases the
average rate also increases, and follows the same trend of the sum rate as
explained in the previous paragraph. However, $\theta_{opt}$ decreases as the
user density increases. This is because, as the number of users in each sector
grows, the amount of power allocated to each user decreases. Consequently, in
order to maximize the sum rate, the objective function drives $\theta$ towards
smaller values. Therefore, we conclude that for a densely populated rural
region, a higher number of sectors would be optimal.

(a) Sum rate per cell.

(b) Average rate per user.

Figure 3: Sum rate per cell and average rate per user in an urban scenario
with $R=10$ m and $h=20$ m as a function of $\theta$, for different user
densities.

In order to demonstrate the impact of the fairness among the users in a cell,
Fig. 2(c) plots Jain’s Fairness index as a function of $\theta$ for different
user densities. It is computed as
$J=\dfrac{\left(\sum_{i}R_{i}\right)^{2}}{M\sum_{i}R_{i}^{2}}$, where $R_{i}$
represents the rate per user [28]. The value of Jain’s fairness index for
$\lambda=\\{0.0005,0.0008,0.001,0.002\\}$ at $\theta_{opt}$ is
$\\{0.96,0.9504,0.9461,0.9345\\}$, respectively. We observe that $J$ increases
with $\theta$, and good fairness is achieved at the optimal sector beamwidth.
This can be attributed to constraint $C_{6}$ in (7), where a minimum data rate
is guaranteed to each user. Thus, we conclude that the proposed resource
allocation ensures fairness among all the users in a cell.

###  IV-B Urban Scenario

Fig. 3(a) plots the sum rate as a function of $\theta$ for various user
densities in a cell of radius $R=10$ m. The UAV is deployed at $h=20$ m. As
observed in Sec. IV-A, the sum rate initially increases and then decreases
exponentially after $\theta$ reaches its optimal value. However,
$\theta_{opt}$ obtained here is generally smaller than in the rural/sub-urban
scenario. At $\lambda=0.05$, $\theta_{opt}=8^{\circ}$ which gives $M=15.7$ and
$S=45$ to serve the users optimally. As the user density increases from
$\lambda=0.08$ to $0.1$, $\theta_{opt}$ decreases from $6^{\circ}$ to
$5^{\circ}$. This implies that the number of sectors required to serve all the
users increase from $S=60$ to $S=72$, respectively. A similar trend for the
value of $\theta_{opt}$ is observed in the plot of the average rate per user
in Fig. 3(b). Notice that these results were obtained up to a beamwidth of
$120^{\circ}$, with the cell having only three sectors to provide a more
comprehensive picture.

Figure 4: Average rate per user vs. $h$ and $\theta$ with $R=100$ m and
$\lambda=0.0005$.

###  IV-C Impact of the UAV Height

Fig. 4 plots the average rate per user as a function of $\theta$ and $h$ in a
rural/sub-urban cell with radius $R=100$ m and $\lambda=0.0005$
UEs/$\text{\,}{\mathrm{m}}^{2}$. The LoS probability ($P_{r}$) as a function
of $h$ defined in Eq. (3) is shown in the inset. We observe that $P_{r}$ is
small for $h\leq 50$ m. As a result, the average rate per user for $h=50$ m is
higher than for $h=10$ m. As $h$ continues to increase, the LoS probability
increases, but so does the path loss. The impact of the path loss becomes the
dominant factor when $h>100$ m, even though the LoS probability also
increases. The effect is also visible in the value of $\theta_{opt}$, which
decreases initially from $12^{\circ}$ to $10^{\circ}$ as $h$ increases from
$50$ m to $100$ m, and then increases to $\theta_{opt}=15^{\circ}$ at $h=200$
m. This is because, for $50\leq h\leq 100$ m, the path loss is relatively
small, thus the optimization problem drives $\theta$ to a smaller value to
increase the antenna gain and the sum rate. However, to compensate for the
very high path loss at $h>100$ m, the optimization problem attempts to
increase the power allocated to the users by increasing $\theta_{opt}$.
Therefore, we conclude that the optimal sector beamwidth is a function of $h$.

In order to compare the optimal number of sectors required to serve a given
number of users in the rural/sub-urban and urban scenarios, we consider two
configurations with the same LoS probability, i.e., $R=100$ m and $h=50$ m,
and $R=10$ m and $h=20$ m, respectively. For the first configuration with
$\lambda=0.0005$ and $M=15.7$, the optimal number of sectors at $\theta_{opt}$
is $S=30$. This is much lower than $S=45$ obtained in the second configuration
with the same number of users in Sec. IV-B. Consequently, the number of
elements required in a UPA antenna decreases from $N=15$ in an urban
environment to around $N=12$ in a sub-urban scenario. This implies a
considerable reduction in the cost and complexity of the antenna system for
serving a rural/sub-urban area.

##  V Conclusion

In this work, we investigated the optimal sector beamwidth for a cell served
by a UAV acting as a BS. Assuming that users are distributed according to a
PPP, we observed that there is an optimal beamwidth to maximize the sum rate
over a region. For a given transmit power, this optimal beamwidth is a
function of the user density and the height of the UAV. Based on simulations
we observed that, for the same LoS probability, $S=45$ sectors are required in
an urban region vs. $S=30$ sectors in a rural region to optimally serve an
average number of $15.7$ users in a cell with the same QoS requirements. Also,
we showed that a remote area with lower density of users can be optimally
served by a UAV-BS flying at lower altitude. This implies lower complexity in
terms of antenna size and transmit power, and less spatial reuse and
interference among the sectors, respectively, compared to the urban
environment. Our conclusions encourage the use of UAV-BSs to connect remote
and poorly served areas.

## References

  * [1] The World Bank, “Percentage of rural population,” Accessed: 09/20/2022. [Online]. Available: https://data.worldbank.org/indicator/SP.RUR.TOTL.ZS
  * [2] ITU, “Measuring digital development: Facts and figures,” _Technical Report_ , 2021, Accessed: 09/20/2022. [Online]. Available: https://www.itu.int/en/ITU-D/Statistics/Pages/facts/
  * [3] M. Giordani, M. Polese, M. Mezzavilla, S. Rangan, and M. Zorzi, “Toward 6G networks: Use cases and technologies,” _IEEE Commun. Mag._ , vol. 58, no. 3, pp. 55–61, Mar. 2020.
  * [4] M. Giordani and M. Zorzi, “Non-terrestrial networks in the 6G era: Challenges and opportunities,” _IEEE Network_ , vol. 35, no. 2, pp. 244–251, Dec. 2020.
  * [5] A. Chaoub, M. Giordani, B. Lall, V. Bhatia, A. Kliks, L. Mendes, K. Rabie, H. Saarnisaari, A. Singhal, N. Zhang _et al._ , “6G for bridging the digital divide: Wireless connectivity to remote areas,” _IEEE Wireless Commun._ , vol. 29, no. 1, pp. 160–168, July 2021.
  * [6] D. Wang, M. Giordani, M.-S. Alouini, and M. Zorzi, “The potential of multi-layered hierarchical non-terrestrial networks for 6G,” _IEEE Veh. Technol. Mag._ , vol. 16, no. 3, pp. 99–107, July 2020.
  * [7] A. U. Rahman, F. Fourati, K.-H. Ngo, A. Jindal, and M.-S. Alouini, “Network graph generation through adaptive clustering and infection dynamics: A step toward global connectivity,” _IEEE Commun. Lett._ , vol. 26, no. 4, pp. 783–787, Jan. 2022.
  * [8] M. Matracia, M. A. Kishk, and M.-S. Alouini, “Coverage analysis for UAV-assisted cellular networks in rural areas,” _IEEE Open J. Veh. Technol._ , vol. 2, pp. 194–206, Apr. 2021.
  * [9] I. A. Elnabty, Y. Fahmy, and M. Kafafy, “A survey on UAV placement optimization for UAV-assisted communication in 5G and beyond networks,” _Phys. Commun._ , vol. 51, p. 101564, Apr. 2022.
  * [10] M. Alzenad, A. El-Keyi, and H. Yanikomeroglu, “3D placement of an unmanned aerial vehicle base station for maximum coverage of users with different QoS requirements,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 1, pp. 38–41, Sept. 2017.
  * [11] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 6, pp. 3747–3760, Mar. 2017.
  * [12] Z. Xiao, H. Dong, L. Bai, D. O. Wu, and X.-G. Xia, “Unmanned aerial vehicle base station (UAV-BS) deployment with millimeter-wave beamforming,” _IEEE Internet of Things J._ , vol. 7, no. 2, pp. 1336–1349, Nov. 2019.
  * [13] R. Duan, J. Wang, C. Jiang, H. Yao, Y. Ren, and Y. Qian, “Resource allocation for multi-UAV aided IoT NOMA uplink transmission systems,” _IEEE Internet of Things J._ , vol. 6, no. 4, pp. 7025–7037, Apr. 2019.
  * [14] N. Rupasinghe, Y. Yapıcı, I. Güvenç, and Y. Kakishima, “Non-orthogonal multiple access for mm-wave drone networks with limited feedback,” _IEEE Trans. Commun._ , vol. 67, no. 1, pp. 762–777, Aug. 2018\.
  * [15] W. Yi, Y. Liu, E. Bodanese, A. Nallanathan, and G. K. Karagiannidis, “A unified spatial framework for UAV-aided mmwave networks,” _IEEE Trans. Commun._ , vol. 67, no. 12, pp. 8801–8817, Oct. 2019.
  * [16] S. Kumar, S. Suman, and S. De, “Dynamic resource allocation in UAV-enabled mmwave communication networks,” _IEEE Internet Things J._ , vol. 8, no. 12, pp. 9920–9933, Sept. 2021.
  * [17] Z. Kaleem, W. Khalid, A. Muqaibel, A. A. Nasir, C. Yuen, and G. K. Karagiannidis, “Learning-aided UAV 3D placement and power allocation for sum-capacity enhancement under varying altitudes,” _IEEE Commun. Lett._ , vol. 26, no. 7, pp. 1633–1637, May 2022.
  * [18] H. He, S. Zhang, Y. Zeng, and R. Zhang, “Joint altitude and beamwidth optimization for UAV-enabled multiuser communications,” _IEEE Commun. Lett._ , vol. 22, no. 2, pp. 344–347, Feb. 2018.
  * [19] Z. Yang, C. Pan, M. Shikh-Bahaei, W. Xu, M. Chen, M. Elkashlan, and A. Nallanathan, “Joint altitude, beamwidth, location, and bandwidth optimization for UAV-enabled communications,” _IEEE Commun. Lett._ , vol. 22, no. 8, pp. 1716–1719, Aug. 2018.
  * [20] N. Tang, H. Tang, B. Li, and X. Yuan, “Joint maneuver and beamwidth optimization for UAV-enabled multicasting,” _IEEE Access_ , vol. 7, pp. 149 503–149 514, Oct. 2019.
  * [21] M. Rawat, B. Lall, and S. Srirangarajan, “Statistical characterization of the angle of arrival in an underwater channel due to scattering from the sea surface modeled as sum of sinusoids,” _IEEE Trans. Veh. Technol._ , vol. 71, no. 2, pp. 1562–1574, Nov. 2021.
  * [22] C. A. Balanis, _Antenna Theory: Analysis and Design_ , 4th ed. John Wiley & Sons, Inc, 2015.
  * [23] Y. Wang, M. Giordani, X. Wen, and M. Zorzi, “On the beamforming design of millimeter wave UAV networks: Power vs. capacity trade-offs,” _Computer Networks_ , vol. 205, p. 108746, Mar. 2022.
  * [24] N. Varshney and S. De, “Optimum downlink beamwidth estimation in mm-wave communications,” _IEEE Trans. Commun._ , vol. 69, no. 1, pp. 544–557, Sept. 2020.
  * [25] I. Mohammed, I. B. Collings, and S. V. Hanly, “Line of sight probability prediction for UAV communication,” in _IEEE ICC Workshops, Montreal, QC, Canada_ , June 2021, pp. 1–5.
  * [26] M. R. Akdeniz, Y. Liu, M. K. Samimi, S. Sun, S. Rangan, T. S. Rappaport, and E. Erkip, “Millimeter wave channel modeling and cellular capacity evaluation,” _IEEE J. Sel. Areas Commun._ , vol. 32, no. 6, pp. 1164–1179, June 2014.
  * [27] T. Zhou, C. Tao, L. Liu, and Z. Tan, “Ricean k-factor measurements and analysis for wideband radio channels in high-speed railway u-shape cutting scenarios,” in _Proc. 79th IEEE Veh. Technol. Conf. (VTC Spring), Seoul, South Korea_ , May 2014, pp. 1–5.
  * [28] J. V. Evangelista, Z. Sattar, G. Kaddoum, and A. Chaaban, “Fairness and sum-rate maximization via joint subcarrier and power allocation in uplink SCMA transmission,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 12, pp. 5855–5867, Sept. 2019.