# Resource Reservation in Backhaul and Radio Access Network with Uncertain
User Demands

Navid Reyhanian, Hamid Farmanbar, and Zhi-Quan Luo N. Reyhanian is with the
Department of Electrical and Computer Engineering, University of Minnesota,
Minneapolis, MN, 55455 USA (e-mail: navid@umn.edu).H. Farmanbar is with Huawei
Canada Research Center, Ottawa, Canada (e-mail:
hamid.farmanbar@huawei.com).Z.-Q. Luo is with Shenzhen Research Institute of
Big Data, The Chinese University of Hong Kong, Shenzhen, China (e-mail:
luozq@cuhk.edu.cn).This paper was presented in part at the $21^{\text{st}}$
IEEE International Workshop on Signal Processing Advances in Wireless
Communications (SPAWC), Atlanta, GA, USA, May 26–29, 2020 [1].

###### Abstract

Resource reservation is an essential step to enable wireless data networks to
support a wide range of user demands. In this paper, we consider the problem
of joint resource reservation in the backhaul and Radio Access Network (RAN)
based on the statistics of user demands and channel states, and also network
availability. The goal is to maximize the sum of expected traffic flow rates,
subject to link and access point budget constraints, while minimizing the
expected outage of downlinks. The formulated problem turns out to be non-
convex and difficult to solve to global optimality. We propose an efficient
Block Coordinate Descent (BCD) algorithm to approximately solve the problem.
The proposed BCD algorithm optimizes the link capacity reservation in the
backhaul using a novel multi-path routing algorithm that decomposes the
problem down to link-level and parallelizes the computation across backhaul
links, while the reservation of transmission resources in RAN is carried out
via a novel scalable and distributed algorithm based on Block Successive
Upper-bound Minimization (BSUM). We prove that the proposed BCD algorithm
converges to a Karush–Kuhn–Tucker (KKT) solution. Simulation results verify
the efficiency and the efficacy of our BCD approach against two heuristic
algorithms.

###### Index Terms:

Resource reservation, multi-path routing, traffic maximization, outage
minimization, parallel computation.

##  I Introduction

Resource reservation is an important step in network planning and management
due to its significant effects on the user quality of service. For wireless
data networks operating in random and dynamic environments, finding resource
reservation protocols that remain robust under uncertain user demands is
challenging. Resource reservation, which balances network performance and its
hardware costs, involves traffic forecasting and resource allocation for the
predicted traffic [2, 3, 4]. Resource reservation in the backhaul and Radio
Access Network (RAN) should satisfy a wide range of applicable traffic
demands. In particular, both the link capacity in the backhaul and
transmission resources in RAN should be sliced and reserved for users such
that upon the arrival of a new demand, the network is able to support it.

Resource reservation for the uncertain demand was first studied by Gomory and
Hu in [5], which reserved link capacities using a single commodity routing
problem with a finite number of sources. For communication networks, where
both link budget and node budget are to be reserved, different approaches are
proposed for resource reservation. In traffic oblivious approaches, to make
reservations and slice the network resources, user demand and its statistics
are not considered in the problem formulation [6, 7, 8]. The drawback of
traffic oblivious approaches is that they limit the ability of a network to
adapt to any given demand. To reserve link capacities in flow networks, a
collection of predicted demand scenarios are considered in [9, 10]. The
proposed algorithms in [9, 10] reserve link capacities such that the predicted
demand scenarios are supported as much as possible. The accuracy of the
reservations in [9, 10] is based on the number of predicted scenarios.
However, as the number of scenarios increases, the complexity of solving the
problem increases. Short term user demands are predicted by Long Short-Term
Memory (LSTM) neural networks in [11, 12, 13]. Recurring resource reservations
based on the short-term traffic variations incur reconfiguration costs,
service interruptions, and overhead in networks [14]. The mean of user demands
is used in [15] to balance the workload among a set of data centers in a
network that consists of the backhaul and RAN such that the utilization of
resources is maximized. The joint reservation of computational and radio
resources is studied in [16], where different ranges are considered for
uncertain user demands. A linear program is formulated in [16] to support the
uncertain user demands, which vary in given ranges, as much as the network
allows. In [17], the transmission resource reservation in RAN is considered
where the minimum requirements of users are known and deterministic. The
authors of [17] proposed a matching-based algorithm to solve an optimization
problem with the goal of minimizing the consumption of network resources while
meeting the requirements of users.

Optimal routing is studied widely for many settings, e.g., [18, 19, 20], while
optimal resource allocation in RAN has also been studied for different
wireless channels, e.g., [12, 21, 22, 23, 24, 25, 26, 17]. The joint routing
in the backhaul network and resource allocation RAN is studied in a number of
more recent papers [27, 28, 29, 30, 31, 32, 33, 34]. In [27] and [31], the
user demand requirements are deterministic and known. On the other hand, in
[28, 30, 34, 32, 33], the traffic of users is maximized as much as the network
is able to support, regardless of user demand statistics. To find a robust
resource reservation, network resources should be reserved based on demand
statistics. In [28, 29, 30] and [34, 32, 33], the wireless channel capacity is
a deterministic function of input power. Moreover, the convexity of the
problem is assumed in [28, 29, 30, 33]. Neither of these assumptions holds in
practice, where the wireless channel capacity is random and its distribution
is a function of supplied transmission resources [35, 36, 37].

In addition to different proposed formulations for resource allocations and
network planning with certain and uncertain user demands in existing
literature, several algorithms have been used to solve the resulting
optimization problems. Among them, the Alternating Direction Method of
Multipliers (ADMM) has been used widely [10, 38, 27, 39, 40]. ADMM enables
flow decoupling in the network optimization process. The efficiency of ADMM
depends on the number of auxiliary link variables introduced to make the
optimization subproblems separable. For networks with a large number of links,
ADMM can be slow, i.e., requiring a large number of iterations. A dual
decomposition method for path-based routing is used in [41], where a gradient
ascent approach has been proposed to solve the dual problem. Since in most
problems the dual function is non-smooth, the gradient ascent approach has to
take small steps, resulting in slow convergence. A distributed approach for
large-scale revenue management problems in airline networks is proposed by
Kemmer et al. in [42]. The single-path dynamic programming approach in [42]
has shown great success in practice despite the absence of convergence or
solution enhancement guarantees.

In this paper, we propose a resource management scheme for end-to-end resource
reservation, i.e., from data centers to users, based on user demand and
downlink achievable rate statistics for a data network consisting of the
backhaul and RAN. We consider a multi-path routing in our formulation, where a
user can be served by several Access Points (APs) through multiple paths from
a data center. We formulate the problem of jointly reserving the transmission
resources in RAN and link capacities in the backhaul based on user demand and
downlink achievable rate statistics so as to maximize the total expected
supportable user traffic, while minimizing the expected outage of downlinks.
Since the formulated problem is non-convex and hard to solve, we propose an
efficient Block Coordinate Descent (BCD) algorithm, which is convergent to a
Karush–Kuhn–Tucker (KKT) solution of the resource reservation problem.

In the proposed BCD approach, one block of variables determines the link
capacity reservation in the backhaul and the other block of variables
specifies the transmission resource reservation in RAN. We alternately
optimize the two blocks of variables in the BCD algorithm. Fixing the
transmission resources in RAN, we update the link capacity reservation in the
backhaul via a novel multi-path routing algorithm. Inspired by the resource
level decomposition ideas in [42], the proposed multi-path routing decomposes
the problem down to link-level and parallelizes the computation across
backhaul links. Based on the convergence theory for Block Successive Upper-
bound Minimization (BSUM) methods in [43], we prove that the proposed multi-
path routing is convergent to the global minima of an arbitrary convex cost
function with Lipschitz continuous gradient. The required computation time for
each iteration of the proposed multi-path routing is equal to that for one
link regardless of the network size. After updating the link capacity
reservations, we update the transmission resource reservation in RAN. Since
the resource reservation problem in RAN is possibly non-convex, we propose a
distributed algorithm based on the BSUM techniques to iteratively solve a
sequence of convex approximations of the original problem. We prove that the
proposed BCD algorithm converges to a KKT solution. To verify the performance
of the proposed algorithm, two heuristic algorithms are also developed and
used as benchmarks to evaluate the efficiency and the efficacy of the proposed
approach via simulations.

The rest of this paper is organized as follows. The system model and problem
formulation are given in Section II. Section III describes a general scalable
and distributed algorithm for the multi-path flow routing. In Section IV, we
propose a BCD algorithm for the network resource reservation problem. The
simulation results are given in Section V, and concluding remarks are given in
Section VI.

##  II System Model and Problem Formulation

Consider a typical scenario whereby user data is transmitted via backhaul
network links from data centers to APs in RAN, which in turn relay the data to
the desired users as depicted in Fig. 1. Suppose $\mathcal{B}$ denotes the set
of APs and $\mathcal{K}$ denotes the set of mobile users. The set of directed
wired links of the backhaul is denoted by $\mathcal{L}$. A path connects a
data center and an AP through a sequence of wired links in the backhaul and
finally goes through one downlink to reach the end user. The downlinks between
APs and users are predetermined according to channel quality, interference
levels, and path loss.

Figure 1: A network comprised of APs and backhaul parts.

We consider each user demands one commodity and there are $K=|\mathcal{K}|$
datastreams in the backhaul network. The proposed scheme can be easily
extended to the scenario that each user demands multiple commodities. To serve
each user, several candidate paths are selected between the origin and
destination, and traffic reservation for the corresponding commodities is
implemented over those paths. The candidate paths can go through different
APs, and the joint transmission of APs to a user (coordinated multi-point
mode) is considered in this paper. Only the last hop on each path is wireless.

Each path is denoted by $p$, and the set of all paths is represented by
$\mathcal{P}$. The set of paths that carry user $k$ data is denoted by
$\mathcal{P}_{k}$. The backhaul network links comprising path $p$ for serving
user $k$ are represented by the set $\mathcal{L}_{k}^{p}$. Similarly, the
network nodes on path $p\in\mathcal{P}_{k}$ are denoted by the set
$\mathcal{U}_{k}^{p}$. The demand of user $k$ is a random variable represented
by $d_{k}$. It follows a certain Probability Density Function (PDF) denoted by
$f_{k}(d_{k})$. The corresponding Cumulative Density Function (CDF) is
represented by $F_{k}(d_{k})$. Let $r_{k}$ denote the traffic rate reserved
for user $k$. The actual traffic flow of user $k$ supported by the network is
a random variable given by

$\min(d_{k},r_{k})=\left\\{\rule{0.0pt}{22.76219pt}\right.\begin{array}[]{@{}cl}r_{k},&\text{if}\>\>\>r_{k}\leq
d_{k},\\\\[8.53581pt] d_{k},&\text{otherwise}.\end{array}$

We calculate the expected supportable traffic rate for user $k$ as follows:

$\displaystyle\mathbb{E}\left(\min(d_{k},r_{k})\right)=\int_{0}^{r_{k}}y_{k}.f_{k}(y_{k})dy_{k}+r_{k}\int_{r_{k}}^{\infty}f_{k}(y_{k})dy_{k}.$

Since the network is not able to support the demand when it exceeds the
reserved rate, we have the minimum in the above expectation. In the first
integral, the random demand of user $k$ falls below the reserved rate. In the
second integral, the random demand exceeds $r_{k}$.

Since a user receives their data from multiple APs, transmission resources
should be reserved in multiple APs for the paths available to the user. The
resource reservation in the backhaul and RAN is limited by two physical
constraints:

  * •

The aggregate reserved traffic rate for paths that share a link must not
exceed the link capacity. Therefore, we have the following constraint:

$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}\leq
C_{l},\hskip 14.22636pt\forall l\in\mathcal{L},$ (1)

where $r_{k}^{p}$ is the reserved traffic rate for path $p$ (for serving user
$k$). Moreover, the capacity of link $l$ is denoted by $C_{l}$. Flows on
different paths available to one user are treated as separate flows. Thus, we
have the inner summation in the above constraint.

  * •

The total reserved transmission resources for different paths must not exceed
the AP capacity. Hence, we have

$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}t_{k}^{p}\leq
C_{b},\hskip 14.22636pt\forall b\in\mathcal{B},$ (2)

where $t_{k}^{p}$ is the reserved transmission resources in AP $b$ to transmit
incoming data from path $p\in\mathcal{P}_{k}$ to user $k$. Moreover, the
capacity of AP $b$ is denoted by $C_{b}$.

In addition to the above physical constraints, our multi-path model enforces
another constraint. Since each datastream originating from a data center
splits into a number of sub-flows, we have the following constraint:

  * •

The aggregate reserved traffic rate for the different paths which carry data
to one user is equal to the reserved rate for that user. Hence, we have the
following constraint:

$\displaystyle\sum_{p\in\mathcal{P}_{k}}r_{k}^{p}=r_{k},\hskip
14.22636pt\forall k.$ (3)

In the considered model, we do not make any assumption about the type of the
transmission resource. It can be bandwidth, transmission power, or time-slot
fraction. Based on the allocated resources, the distribution of the achievable
rate of a downlink follows a particular PDF. As only the last hop on each path
is wireless, path $p$ uniquely identifies the downlink of the last hop. The
achievable rate (i.e., instantaneous capacity) of the downlink of path $p$ is
random and follows an arbitrary distribution with a PDF represented by
$z_{k}^{p}(v_{k}^{p},t_{k}^{p})$ and a CDF denoted by
$Z_{k}^{p}(v_{k}^{p},t_{k}^{p})$. The PDF is a function of two variables: the
achievable rate of the downlink, denoted by $v_{k}^{p}$, and the allocated
transmission resource, denoted by $t_{k}^{p}$. When the achievable rate of a
downlink falls below the reserved rate $r_{k}^{p}$, some outage is experienced
and its amount is $r_{k}^{p}-v_{k}^{p}$, given that the amount of allocated
transmission resources to the downlink is $t_{k}^{p}$. The probability that
this amount of outage takes place is $z_{k}^{p}(v_{k}^{p},t_{k}^{p})$. In
light of the above arguments, the expected outage of the downlink of path $p$
is obtained as follows:

$\displaystyle\int_{0}^{r_{k}^{p}}z_{k}^{p}(v_{k}^{p},t_{k}^{p})\>(r_{k}^{p}-v_{k}^{p})dv_{k}^{p}.$
(4)

Since the achievable rate is a continuous random variable, we have the above
integral.

In this paper, we aim to maximize the expected traffic of users as much as the
network is able to support, while minimizing the expected outage of downlinks.
We formulate the following optimization problem to find resource reservations
in the backhaul and RAN:

$\displaystyle\underset{\mathbf{r,t}}{\text{max}}$
$\displaystyle\sum_{k=1}^{K}\Big{[}\mathbb{E}[\min(r_{k},d_{k})]-\theta_{k}\sum_{p\in\mathcal{P}_{k}}\int_{0}^{r_{k}^{p}}z_{k}^{p}(v_{k}^{p},t_{k}^{p})\>(r_{k}^{p}-v_{k}^{p})dv_{k}^{p}\Big{]}$
s.t.
$\displaystyle\eqref{eq:linkcap},\eqref{eq:nodecap},\eqref{eq:split},r_{k},r_{k}^{p},t_{k}^{p}\geq
0,\hskip 14.22636ptp\in\mathcal{P}_{k},\forall k,$ (5)

where $\theta_{k}:\theta_{k}\geq 0$ is a coefficient chosen by the system
designer that adjusts the priorities of maximizing the expected supportable
traffic of user $k$ and the minimization of the aggregate outage of downlinks,
which serve user $k$. The two blocks of variables in the above problem are
$\mathbf{r}=\\{r_{k},r_{k}^{p}\\}_{p\in\mathcal{P}_{k},k=1:K}$ and
$\mathbf{t}=\\{t_{k}^{p}\\}_{p\in\mathcal{P}_{k},k=1:K}$.

######  Remark 1.

Suppose that multiple paths available to user $k$ share a downlink (the last
hop). The aggregate outage of downlinks for serving user $k$ is calculated as
follows:

$\displaystyle\sum_{w\in\mathcal{W}_{k}}\int_{0}^{\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}r_{k}^{p}}z_{k}^{w}(v_{k}^{w},t_{k}^{w})\times(\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}r_{k}^{p}-v_{k}^{w})dv_{k}^{w},$ (6)

where $\mathcal{W}_{k}$ is the set of downlinks, each denoted by $w$, for
serving user $k$. When multiple paths available to user $k$ share a downlink,
the above outage is placed in the objective function of (5) instead of its
second term, which includes (4).

The maximization problem (5) is not easy to solve to global optimality. The
objective function of (5) is in general not necessarily jointly concave in
$\mathbf{r}$ and $\mathbf{t}$ for an arbitrary PDF
$z_{k}^{p}(v_{k}^{p},t_{k}^{p})$. The reason is that
$\int_{0}^{r_{k}^{p}}\partial^{2}z_{k}^{p}(v_{k}^{p},t_{k}^{p})/(\partial
t_{k}^{p})^{2}\>(r_{k}^{p}-v_{k}^{p})dv_{k}^{p}$ is not always non-negative.

######  Proposition 1.

Given $\mathbf{t}$, the optimization in (5) becomes concave in $\mathbf{r}$.

###### Proof.

Fixing $\mathbf{t}$, the objective function is separable in $k$. We find the
Hessian with respect to $r_{k}$ and $\\{r_{k}^{p}\\}_{p\in\mathcal{P}_{k}}$
for those objective function terms which are associated with user $k$ as
follows:

$\displaystyle\mathbf{H}_{k}=\begin{pmatrix}-f_{k}(r_{k})&0&\dots&0\\\
0&-z_{k}^{1}(r_{k}^{1},t_{k}^{1})&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\dots&-z_{k}^{|\mathcal{P}_{k}|}(r_{k}^{|\mathcal{P}_{k}|},t_{k}^{|\mathcal{P}_{k}|})\end{pmatrix}.$

The overall Hessian matrix is

$\displaystyle\mathbf{H}=\begin{pmatrix}\mathbf{H}_{1}&&\mathbf{0}\\\
&\ddots&\\\ \mathbf{0}&&\mathbf{H}_{K}\end{pmatrix}.$

It is observed that the above matrix is negative semidefinite. Since the
constraints of problem (5) are all affine, it follows that the maximization
(5) is concave with fixed t. ∎

Separable constraints on r and t in (5) motivate the BCD algorithm. It is
straightforward to show that with (6) instead of (4) in the objective
function, the optimization in (5) remains concave in r.

##  III Distributed Multi-Path Routing in the Backhaul

This section is concerned with solving (5) when $\mathbf{t}$ is kept fixed,
and (5) is converted to the minimization format after multiplying the
objective function by $-1$. In particular, we study a general multi-path
routing to minimize any convex cost function with a Lipschitz continuous
gradient. We develop an algorithm that is dual-based and decomposes the
problem down to link-level and parallelizes computations across links of the
network. The required computation time for each iteration of the proposed
multi-path routing algorithm is equal to that for one link regardless of the
network size. This interesting property makes the proposed algorithm
appropriate for the online optimization of large networks.

For each datastream in the network, several candidate paths are selected. We
assume that each flow can be split into multiple sub-flows. To formulate the
multi-path routing problem, we first assume that the cost function is
separable in variables, i.e.,
$\psi(\mathbf{r})=\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\psi_{k}^{p}(r_{k}^{p})$,
where each $\psi_{k}^{p}(r_{k}^{p})$ is strictly convex.

The optimization problem for the multi-path flow routing can be written as
follows:

$\displaystyle\underset{\mathbf{r}}{\min}$
$\displaystyle\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\psi_{k}^{p}(r_{k}^{p})$
(7) s.t. $\displaystyle\eqref{eq:linkcap},r_{k}^{p}\geq 0,\hskip
14.22636ptp\in\mathcal{P}_{k},\forall k.$

Since typically the number of variables is greater than the number of
constraints in the above optimization, solving the problem is easier in the
dual domain. The Lagrangian function for the above problem is

$\displaystyle
L_{c}(\mathbf{r},\boldsymbol{\mu},\boldsymbol{\phi})=\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\psi_{k}^{p}(r_{k}^{p})+\sum_{l\in\mathcal{L}}\mu_{l}(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l})$
$\displaystyle-\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\phi_{k}^{p}r_{k}^{p},$
(8)

where $\mu_{l}:\mu_{l}\geq 0$ is the Lagrange multiplier for the capacity
constraint of link $l$, and $\phi_{k}^{p}:\phi_{k}^{p}\geq 0$ is the Lagrange
multiplier for constraint $r_{k}^{p}\geq 0$. Furthermore,
$\boldsymbol{\mu}=\\{\mu_{l}\\}_{l\in\mathcal{L}}$ and
$\boldsymbol{\phi}=\\{\phi_{k}^{p}\\}_{p\in\mathcal{P}_{k},k=1:K}$. We find
the dual problem of (7) as follows:

$\displaystyle\underset{\boldsymbol{\mu},\boldsymbol{\phi}}{\text{max}}$
$\displaystyle\min_{\mathbf{r}}L_{c}(\mathbf{r},\boldsymbol{\mu},\boldsymbol{\phi})$
(9) s.t.
$\displaystyle\boldsymbol{\mu}\geq\mathbf{0},\boldsymbol{\phi}\geq\mathbf{0}.$

For many cost functions, no closed-form solution for
$\mathbf{r}=\arg\min_{\mathbf{r}}L_{c}(\mathbf{r},\boldsymbol{\mu},\boldsymbol{\phi})$
exists. Therefore, commonly, the above problem is solved via a primal-dual
method such as ADMM [10, 38, 39, 40]. However, the auxiliary link variables
introduced to make the per-flow subproblems of optimization in (7) separable
can slow down ADMM in practice.

Resource level decomposition for large-scale single-path applications was
first proposed in [42] to solve the revenue management problems in airline
networks. The proposed decomposition in [42] does not involve any auxiliary
variables. In spite of the absence of convergence or solution enhancement
guarantees, the resource level decomposition has been rather successful in
practice. We leverage resource level decomposition ideas to develop a
distributed algorithm to solve the general multi-path routing problem (7) in a
parallel fashion such that the traffic passing on each link can be obtained
independently from the other links. Unlike the dynamic programming approach in
[42], an optimization-based approach is proposed here to solve subproblems. In
each iteration, the proposed dual algorithm decomposes the problem in (9) and
solves the subproblems globally and in parallel. The optimized $\mu_{l}$ in
the $j^{\text{th}}$ iteration of the proposed algorithm is denoted by
$\mu_{l}^{j}$. Here, we explain the decomposition. The dualized link capacity
constraints
$\sum_{l\in\mathcal{L}}\mu_{l}(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l})$
in the Lagrangian (8) are separable across links. Each link $l$ receives
$\mu_{l}(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l})$.
In each iteration, based on
$\boldsymbol{\mu}^{j-1}=\\{\mu_{l}^{j-1}\\}_{l\in\mathcal{L}}$ in the previous
iteration, we decompose the non-separable terms in the Lagrangian (8), which
include $r_{k}^{p}$, across links on path $p$. Each link $l$ of path $p$
receives a portion of

$\displaystyle\alpha_{k,l}^{p,j}=\mu_{l}^{j-1}/\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\mu_{l^{\prime}}^{j-1},$
(10)

In the $j^{\text{th}}$ iteration, the decomposed per-link Lagrangian function
is as follows:

$\displaystyle
L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})=\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}\alpha_{k,l}^{p,j}\psi_{k}^{p}(r_{k}^{p})$
$\displaystyle+\mu_{l}(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l})-\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}\alpha_{k,l}^{p,j}\phi_{k}^{p}r_{k}^{p},$
(11)

where
$\mathbf{r}_{l}=\\{r_{k}^{p}\\}_{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p},k=1:K}$
and
$\boldsymbol{\phi}_{l}=\\{\phi_{k}^{p}\\}_{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p},k=1:K}$.
We notice that based on (10),
$\\{\alpha_{k,l}^{p,j}\\}_{l\in\mathcal{L}_{k}^{p}}$ in (11) is calculated
using $\boldsymbol{\mu}^{j-1}$. Based on the above decomposition, we obtain

$\displaystyle
L_{c}(\mathbf{r},\boldsymbol{\mu},\boldsymbol{\phi})=\sum_{l\in\mathcal{L}}L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1}).$
(12)

Instead of solving the problem in (9), we solve

$\displaystyle\underset{\\{\mu_{l},\boldsymbol{\phi}_{l}\\}_{l\in\mathcal{L}}}{\text{max}}$
$\displaystyle\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})$
(13) s.t. $\displaystyle\mu_{l}\geq
0,\boldsymbol{\phi}_{l}\geq\mathbf{0},\hskip 14.22636ptl\in\mathcal{L},$

iteratively and then update $\alpha_{k,l}^{p,j+1}$ for iteration $j+1$. The
above problem is decomposable in $\\{\mu_{l},\boldsymbol{\phi}_{l}\\}$ and can
be solved in parallel for all links. Due to strong duality [44, p. 226–p.
227], each subproblem of (13) is equivalent to the following per-link problem
in the primal domain:

$\displaystyle\underset{\mathbf{r}_{l}}{\min}$
$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}\alpha_{k,l}^{p,j}\psi_{k}^{p}(r_{k}^{p})$
(14) s.t.
$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}\leq
C_{l},$ $\displaystyle\alpha_{k,l}^{p,j}r_{k}^{p}\geq 0,\hskip
14.22636ptp\in\mathcal{P}_{k},\forall k,l\in\mathcal{L}_{k}^{p}.$

The optimal $r_{k}^{p}$ and $\mu_{l}$ can be obtained using the first-order
optimality condition for the per-link subproblem in (14). Here, we list KKT
conditions as follows:

$\displaystyle\frac{\partial
L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})}{\partial
r_{k}^{p}}=\alpha_{k,l}^{p,j}\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\mu_{l}-\alpha_{k,l}^{p,j}\phi_{k}^{p}=0,$ (15a)
$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}\leq
C_{l},$ (15b)
$\displaystyle\mu_{l}\>(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l})=0,\>\>\mu_{l}\geq
0,$ (15c)
$\displaystyle\alpha_{k,l}^{p,j}r_{k}^{p}\>\phi_{k}^{p}=0,\>\>r_{k}^{p}\geq
0,\>\>\phi_{k}^{p}\geq 0.$ (15d)

First, we consider that $r_{k}^{p}>0$ and $\phi_{k}^{p}=0$. Due to the strict
convexity of $\psi_{k}^{p}(r_{k}^{p})$,
$\partial\psi^{p}_{k}(r_{k}^{p})/\partial r^{p}_{k}$ is strictly increasing.
Thus, given $\mu_{l}$, there is a unique $r_{k}^{p}$ to solve
$\alpha_{k,l}^{p,j}\partial\psi_{k}^{p}(r_{k}^{p})/\partial
r_{k}^{p}+\mu_{l}=0$. Since $\partial\psi^{p}_{k}(r_{k}^{p})/\partial
r^{p}_{k}$ is strictly increasing, we implement a bisection search on
$r_{k}^{p}$ in the non-negative orthant $r_{k}^{p}\geq 0$ to find $r_{k}^{p}$
from $\alpha_{k,l}^{p,j}\partial\psi_{k}^{p}(r_{k}^{p})/\partial
r_{k}^{p}+\mu_{l}=0$. If the obtained $r_{k}^{p}$ is positive, we keep
$\phi_{k}^{p}=0$. Otherwise, we set $r_{k}^{p}=0$ and find
$\phi_{k}^{p}=\partial\psi_{k}^{p}(r_{k}^{p})/\partial
r_{k}^{p}|_{r_{k}^{p}=0}+\mu_{l}/\alpha_{k,l}^{p,j}$. For a given $\mu_{l}$,
we obtain each $r_{k}^{p}$ variable associated with link $l$, i.e.,
$r_{k}^{p}:p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p},k=1:K$. The dual
approach for solving the optimization in (14) works as follows: implement a
bisection search on the Lagrange multiplier $\mu_{l}$ in the positive orthant
and numerically find each $r_{k}^{p}$ variable from (15a) and (15d) for each
$\mu_{l}$ until we have
$\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}=C_{l}$.
If there is no such positive $\mu_{l}$, we drop the first constraint from
optimization (14) and solve (14) by setting the gradient of the cost function
to zero. Then, we project the solution to the positive orthant. Due to the
strict convexity of each subproblem, the optimal primal variables are unique.
The optimal $\mu_{l}$ for each per-link subproblem is also unique. We justify
this claim.

0\. Initialization $s_{1}=0$, $s_{2}=\text{large number}$, $q_{1}=0$,
$q_{2}=0$, $q_{3}=0$;

repeat

1\. $s_{3}=(s_{1}+s_{2})/2$;

2\. Implement a bisection search to solve (15a) with $\phi_{k}^{p}=0$ and find
$r_{k}^{p}:r_{k}^{p}\geq 0$, where $\mu_{l}=s_{1}$;

3\. if  _there is no positive solution for $r_{k}^{p}$_ then

$r_{k}^{p}=0$ and $\phi_{k}^{p}=\partial\psi_{k}^{p}(r_{k}^{p})/\partial
r_{k}^{p}|_{r_{k}^{p}=0}+\mu_{l}/\alpha_{k,l}^{p,j}$;

4\.
$q_{1}=\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l}$;

5\. Implement a bisection search to solve (15a) with $\phi_{k}^{p}=0$ and find
$r_{k}^{p}:r_{k}^{p}\geq 0$, where $\mu_{l}=s_{2}$;

6\. if  _there is no positive solution for $r_{k}^{p}$_ then

$r_{k}^{p}=0$ and $\phi_{k}^{p}=\partial\psi_{k}^{p}(r_{k}^{p})/\partial
r_{k}^{p}|_{r_{k}^{p}=0}+\mu_{l}/\alpha_{k,l}^{p,j}$;

7\.
$q_{2}=\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l}$;

8\. Implement a bisection search to solve (15a) with $\phi_{k}^{p}=0$ and find
$r_{k}^{p}:r_{k}^{p}\geq 0$, where $\mu_{l}=s_{3}$;

9\. if  _there is no positive solution for $r_{k}^{p}$_ then

$r_{k}^{p}=0$ and $\phi_{k}^{p}=\partial\psi_{k}^{p}(r_{k}^{p})/\partial
r_{k}^{p}|_{r_{k}^{p}=0}+\mu_{l}/\alpha_{k,l}^{p,j}$;

10\.
$q_{3}=\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l}$;

11\. if  _$q_{1}.q_{3} <0$_ then

$s_{2}=s_{3}$;

12\. if  _$q_{2}.q_{3} <0$_ then

$s_{1}=s_{3}$;

13\. if  _$q_{1} <0$, $q_{2}<0$, $q_{3}<0$ _ then

13.1. $\mu_{l}=0$;

13.2. Solve $\partial\psi_{k}^{p}(r_{k}^{p})/\partial r_{k}^{p}=0$ to find
$r_{k}^{p}$;

13.3. Project the obtained $r_{k}^{p}$ variable to the positive orthant;

13.4. $s_{2}=s_{1}$;

until  _$s_{2}-s_{1}$ is small enough_;

Algorithm 1 Dual algorithm to solve the per-link optimization in (14)

If the link capacity constraint is not tight, then due to (15c), $\mu_{l}$ has
to be zero. If the link capacity is tight, then at least one
$r_{k}^{p}:p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p},k=1:K$ is non-zero and
$\phi_{k}^{p}=0$. Due to a) the strict convexity of $\psi^{p}_{k}(r_{k}^{p})$
and the monotone variation of $\psi^{p}_{k}(r_{k}^{p})/\partial r^{p}_{k}$;
and b) the uniqueness of the optimal $r_{k}^{p}$, the obtained $\mu_{l}$ from
(15a) is unique. We justify the bisection search on $\mu_{l}$ as follows: if
the unique optimal $\mu_{l}$ is positive, from (15c), we observe that we must
have
$\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}=C_{l}$,
where each $r_{k}^{p}$ is found from (15a) and (15d). Such positive $\mu_{l}$
can be uniquely found using a bisection search due to the strictly monotone
variation of
$\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p},$ with
$\mu_{l}$ (strict convexity of $\psi_{k}^{p}(r_{k}^{p})$ as explained above).
If the optimal $\mu_{l}$ is zero, then (15b) and (15c) are already satisfied
and it is enough to find the unique non-negative minimizer of each
$\psi_{k}^{p}(r_{k}^{p})$ from (15a) and (15d). In light of the above
arguments, two nested bisection methods are required to solve (14): the inner
bisection works on $r_{k}^{p}$ and the outer one works on $\mu_{l}$. The
summary of the proposed bisection approach to solve the per-link optimization
in (14) is given in Algorithm 1.

0\. Initialization Assign some small positive number to each $\mu_{l}^{0}$,
$j=0$;

repeat

for  _all links_ do

1\. Find
$\alpha_{k,l}^{p,j+1}=\mu_{l}^{j}/\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\mu_{l^{\prime}}^{j}$;

if  _$\mu_{l}^{j} >0$_ then

2\. Apply Algorithm 1 to find $\mu_{l}^{j+1}$;

3\. $j=j+1$;

until  _$\boldsymbol{\mu}^{j}$ converge_;

for  _all $\\{r_{k}^{p}\\}_{p\in\mathcal{P}_{k},k=1:K}$ variables_ do

4\. Use the latest computed $r_{k}^{p}$ by Algorithm 1 from a per-link
subproblem, where $l\in\mathcal{L}_{k}^{p}$ and $\mu_{l}^{j}>0$;

Algorithm 2 Multi-path routing algorithm to solve the optimization in (7)

Suppose that the optimization in (14) is iteratively solved in parallel for
all links of the network. For a link with a large capacity, the link capacity
constraint is not tight and Algorithm 1 finds $\mu_{l}^{j}=0$ and we have
$\alpha_{k,l}^{p,j+1}=0$. For those links, we do not need to continue
computation as the KKT conditions listed in (15a)–(15d) remain satisfied. In
the following iterations, we ignore those links and consider links with
$\mu_{l}^{j}>0$. We alternate between solving the optimization in (14) in
parallel for all links and updating $\alpha_{k,l}^{p,j+1}$ until all
$\\{\mu_{l}^{j}\\}_{l\in\mathcal{L}}$ variables converge, i.e.,
$\left\lVert\boldsymbol{\mu}^{j}-\boldsymbol{\mu}^{j-1}\right\rVert_{2}<\epsilon$.
Once $\boldsymbol{\mu}^{j}$ converges, for each $r_{k}^{p}$ variable, we use
the computed $r_{k}^{p}$ in the last iteration of Algorithm 1 from a
subproblem with $\mu_{l}^{j}>0,l\in\mathcal{L}_{k}^{p}$. A brief description
of the proposed dual algorithm for solving the optimization in (7) is given in
Algorithm 2.

After Algorithm 2 converges, we use the obtained $r_{k}^{p}$ from a per-link
problem with tight link capacity constraint, i.e., $l:\mu_{l}^{j}>0$, for the
other links on that path for which Algorithm 2 finds $\mu_{l}^{j}=0$. The key
property of Algorithm 2 is that after convergence, the obtained $r_{k}^{p}$ on
different links of one path are identical.

######  Proposition 2.

Upon convergence of Algorithm 2, the flow rates across links on each path are
identical.

###### Proof.

Algorithm 2 finds $r_{k}^{p}$, from the per-link subproblem for link
$l\in\mathcal{L}_{k}^{p},\mu_{l}^{j}>0$, using the following equation:

$\displaystyle\frac{\partial
L_{l}(\mathbf{r}_{l},\mu_{l}^{j},\boldsymbol{\phi}_{l}^{j},\boldsymbol{\mu}^{j-1})}{\partial
r_{k}^{p}}=\alpha_{k,l}^{p,j}\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\mu_{l}^{j}-\alpha_{k,l}^{p}\phi_{k}^{p,j}$
$\displaystyle=\frac{\mu_{l}^{j-1}}{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\mu_{l^{\prime}}^{j-1}}\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\mu_{l}^{j}-\frac{\mu_{l}^{j-1}}{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\mu_{l^{\prime}}^{j-1}}\phi_{k}^{p,j}=0.$
(16)

Suppose Algorithm 2 has converged in the $j^{\text{\text{th}}}$ iteration; we
have
$\left\lVert\boldsymbol{\mu}^{j}-\boldsymbol{\mu}^{j-1}\right\rVert_{2}<\epsilon$.
Then, we have
$\boldsymbol{\mu}^{j}-\boldsymbol{\mu}^{j-1}=\boldsymbol{\vartheta}$, where
$\left\lVert\boldsymbol{\vartheta}\right\rVert_{2}<\epsilon$. We have

$\displaystyle\frac{1}{\alpha_{k,l}^{p,j}}=\frac{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\mu_{l^{\prime}}^{j-1}}{\mu_{l}^{j-1}}=\frac{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}(\mu_{l^{\prime}}^{j}-\vartheta_{l^{\prime}})}{\mu_{l}^{j}-\vartheta_{l}}.$

We multiply (16) by $1/\alpha_{k,l}^{p,j}$ and we have the following:

$\displaystyle\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\frac{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}(\mu_{l^{\prime}}^{j}-\vartheta_{l^{\prime}})}{\mu_{l}^{j}-\vartheta_{l}}\mu_{l}^{j}-\phi_{k}^{p,j}$
$\displaystyle=\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}(\mu_{l^{\prime}}^{j}-\vartheta_{l^{\prime}})(1+\frac{\vartheta_{l}}{\mu_{l}^{j}-\vartheta_{l}})-\phi_{k}^{p,j}$
$\displaystyle=\underbrace{\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\sum_{l\in\mathcal{L}_{k}^{p}}\mu_{l}^{j}-\phi_{k}^{p,j}}_{\partial
L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})/\partial
r_{k}^{p}}-\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\vartheta_{l^{\prime}}(1+\frac{\vartheta_{l}}{\mu_{l}^{j}-\vartheta_{l}})$
$\displaystyle+\frac{\vartheta_{l}(\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\mu_{l^{\prime}}^{j})}{\mu_{l}^{j}-\vartheta_{l}}=\frac{\partial
L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})}{\partial
r_{k}^{p}}-\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\vartheta_{l^{\prime}}(1+\frac{\vartheta_{l}}{\mu_{l}^{j}-\vartheta_{l}})$
$\displaystyle+\frac{\vartheta_{l}(\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\mu_{l^{\prime}}^{j})}{\mu_{l}^{j}-\vartheta_{l}}=0.$
(17)

When $\epsilon$ tends to zero, then
$\boldsymbol{\vartheta}\rightarrow\mathbf{0}$ and from (17) we find $\partial
L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})/\partial
r_{k}^{p}=0$. Moreover, we observe that $\partial
L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})/\partial
r_{k}^{p}$ is independent of the link index on path $p$. This means that
$\\{r_{k}^{p}\\}_{p\in\mathcal{P}_{k}}$ variables obtained by solving the link
subproblems are identical for all links along each path $p$ for which
$\mu_{l}^{j}>0$. They are also equal to the minimizer of Lagrangian function
$L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})$ in (8). ∎

######  Theorem 1.

If $\psi(\mathbf{r})$ is strictly convex and separable, then the primal and
dual iterates of Algorithm 2 will converge to the optimal primal and dual
solutions of (7).

###### Proof.

Notice that based on the definition of $\alpha_{k,l}^{p,j}$, given identical
feasible variables $r_{k}^{p}$, $\hat{\mu}_{l}$ and
$\hat{\boldsymbol{\phi}_{l}}$ to both Lagrangian functions in (8) and (11),
from (12), we have
$\sum_{l\in\mathcal{L}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})=L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$.
First, we show that
$\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})$
is a lower-bound for
$\min_{\mathbf{r}}L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$.
Since the minimum of
$L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})$
is less than or equal to the other values of
$L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})$,
we have

$\displaystyle\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})\leq
L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1}).$

Thus, we obtain

$\displaystyle\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})$
$\displaystyle\leq\sum_{l\in\mathcal{L}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})=L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}}),$

where the equality is due to (12). In
$L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$, we choose
$\mathbf{r}$ to be the minimizer of
$L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$. Thus, we
obtain

$\displaystyle\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})\leq\min_{\mathbf{r}}L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}}).$
(18)

From (18), we observe that solving the problem in (13) iteratively is a
successive lower-bound maximization (upper-bound minimization if we rewrite
problems (9) and (13) as minimizations).

We justify the claim that the primal and dual solutions obtained from solving
(13) successively converge to the primal and dual solutions of (7). We build
our proof based on the convergence theory for BSUM given in [43]. We show, in
the same order given in the Appendix, that the lower-bound satisfies all four
convergence conditions given in [43, Assumption 2]:

  1. 1.

At feasible points $\hat{\boldsymbol{\mu}}\geq\mathbf{0}$ and
$\hat{\boldsymbol{\phi}}\geq\mathbf{0}$, we show that
$\min_{\mathbf{r}}L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})=\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$.
From KKT conditions for each subproblem, we obtain

$\displaystyle\frac{L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})}{\partial
r_{k}^{p}}=\frac{\hat{\mu}_{l}}{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\hat{\mu}_{l^{\prime}}}\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\hat{\mu}_{l}$
$\displaystyle-\frac{\hat{\mu}_{l}}{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\hat{\mu}_{l^{\prime}}}\hat{\phi}_{k}^{p}=0.$
(19)

Assuming $\hat{\mu}_{l}>0$, after multiplication by
$\frac{\sum_{l^{\prime}\in\mathcal{L}_{k}^{p}}\hat{\mu}_{l^{\prime}}}{\hat{\mu}_{l}}$,
we obtain

$\displaystyle\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\sum_{l\in\mathcal{L}_{k}^{p}}\hat{\mu}_{l}-\hat{\phi}_{k}^{p}=\frac{\partial
L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})}{\partial
r_{k}^{p}}=0.$ (20)

When $\psi(\mathbf{r})$ is strictly convex, there is a unique minimizer for
each $L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$ and
$L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$.
We observe from (19) and (20) that, at feasible points
$\hat{\boldsymbol{\mu}}\geq\mathbf{0}$ and
$\hat{\boldsymbol{\phi}}\geq\mathbf{0}$, the minimizer of
$L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$ is equal to
that of
$L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$.
From (12), applying identical variables $r_{k}^{p}$, $\hat{\mu}_{l}$ and
$\hat{\phi}_{l}$ to both
$L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$ and
$L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$,
we have
$L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}},\hat{\boldsymbol{\mu}})=\sum_{l\in\mathcal{L}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$.
We choose each $r_{k}^{p}$ to be the minimizer, and we find
$\min_{\mathbf{r}}L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}},\hat{\boldsymbol{\mu}})=\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$.

  2. 2.

From (18), we observe that
$\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})$
is a lower-bound.

  3. 3.

We deploy [45, Proposition 7.1.1] to find the derivative of
$\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$
with respect to $\mu_{l}$. There are three satisfied conditions that ensure
the existence of the derivative: a) the feasible set of (14) is compact; b)
$L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$
is continuous in $\mu_{l}$; and c) for each $\hat{\mu}_{l}$, the equation
$\partial
L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})/\partial
r_{k}^{p}=0$ has a unique solution for $r_{k}^{p}$ due to the strict convexity
of $\psi(\mathbf{r})$. Given identical $\hat{\mu}_{l}$ and $\hat{\phi}_{l}$ to
both Lagrangian functions (8) and (11), the derivative of
$\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$
with respect to $\mu_{l}$ is

$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l},$
$\displaystyle\hskip 59.75095pt\text{where}\hskip
2.84544ptr_{k}^{p}=\arg\min_{r_{k}^{p}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}}).$

The derivative of
$\min_{\mathbf{r}}L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$
with respect to $\mu_{l}$ is

$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},l\in\mathcal{L}_{k}^{p}\\}}r_{k}^{p}-C_{l},\hskip
2.84544pt\text{where}\hskip
2.84544ptr_{k}^{p}=\arg\min_{r_{k}^{p}}L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}}).$

As the minimizers of
$L_{c}(\mathbf{r},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\phi}})$ and
$L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\hat{\boldsymbol{\mu}})$
are equal at point $\hat{\boldsymbol{\mu}}$ due to (19) and (20), we observe
that both above derivatives are equal.

  4. 4.

$\sum_{l\in\mathcal{L}}\min_{\mathbf{r}_{l}}L_{l}(\mathbf{r}_{l},\hat{\mu}_{l},\hat{\boldsymbol{\phi}}_{l},\boldsymbol{\mu}^{j-1})$
is a piecewise linear function of $\hat{\boldsymbol{\mu}}$, and thus, it is a
continuous function of $\hat{\boldsymbol{\mu}}$.

0\. Initialization Choose a feasible vector $\mathbf{r}^{0}$, $m=0$;

repeat

1\. Find the upper-bound (21) using $\mathbf{r}^{m}$;

2\. Apply Algorithm 2 to find $\mathbf{r}$;

3\. $m=m+1$ and $\mathbf{r}^{m}=\mathbf{r}$;

until  _variables in $\mathbf{r}^{m}$ converge_;

Algorithm 3 Multi-path routing algorithm for non-separable cost functions

Building on the above arguments, Algorithm 2 is a block successive lower-bound
maximization method, which satisfies all four convergence conditions given in
[43, Assumption 2]. Algorithm 2 converges to the global optimal solution of
the concave problem (9) [43, Theorem 2], which has an identical objective
function to (7) at the optimal point as a result of strong duality [44, p.
226–p. 227]. Once Algorithm 2 converges, $\boldsymbol{\mu}^{j}$ and
$r_{k}^{p}=\arg\min_{r_{k}^{p}}L_{l}(\mathbf{r}_{l},\mu_{l}^{j},\boldsymbol{\phi}_{l}^{j},\boldsymbol{\mu}^{j-1})$
by Algorithm 2 satisfy (15a)–(15d). The KKT conditions for (7) are (15b)–(15d)
in addition to $\frac{\partial
L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})}{\partial
r_{k}^{p}}=0$. Due to (16) and (17), when Algorithm 2 converges, minimizers of
$L_{l}(\mathbf{r}_{l},\mu_{l}^{j},\boldsymbol{\phi}_{l}^{j},\boldsymbol{\mu}^{j-1})$
and $L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})$ are
identical, and thus, (15a) ensures $\frac{\partial
L_{c}(\mathbf{r},\boldsymbol{\mu}^{j},\boldsymbol{\phi}^{j})}{\partial
r_{k}^{p}}=0$. Hence, the primal and dual variables obtained by Algorithm 2
satisfy the KKT conditions for (7). ∎

######  Remark 2.

If the cost function is convex and has a gradient that is Lipschitz
continuous, but the function is not separable in $r_{k}^{p}$, i.e.,
$\psi(\mathbf{r})$ cannot be written as
$\psi(\mathbf{r})=\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\psi_{k}^{p}(r_{k}^{p})$,
we use the quadratic upper-bound given in [46, eq. (12)], which is separable
in variables. For an arbitrary convex cost function with a Lipschitz
continuous gradient like $\psi(\mathbf{r})$, we have the following upper-
bound:

$\displaystyle\psi(\mathbf{\mathbf{r}})\leq\psi(\mathbf{r}^{m})+\nabla\psi(\mathbf{r}^{m})^{T}(\mathbf{r}-\mathbf{r}^{m})+\frac{\gamma}{2}\left\lVert\mathbf{r}-\mathbf{r}^{m}\right\rVert_{2}^{2},$
(21)

where $\gamma$ is the Lipschitz constant, and
$\mathbf{r}^{m}=\\{r_{k}^{p,m}\\}_{p\in\mathcal{P}_{k},k=1:K}$ is the
$m^{\text{th}}$ iterate in the successive upper-bound minimization. We start
from an initial point $\mathbf{r}^{0}$ in the feasible set and find the upper-
bound (21). Then, we apply Algorithm 2 to solve the problem with the upper-
bound (21) to the global optimal solution in a parallel fashion.

When the upper-bound (21) is substituted for the cost function, the first KKT
condition is

$\displaystyle\frac{\partial
L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})}{\partial
r_{k}^{p}}=$
$\displaystyle\alpha_{k,l}^{p,j}\left(\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}\mid_{r_{k}^{P}=r_{k}^{p,m}}+\gamma(r_{k}^{p}-r_{k}^{p,m})\right)+\mu_{l}-\alpha_{k,l}^{p,j}\phi_{k}^{p}=0,$
(22)

instead of (15a). Once the problem with the upper-bound (21) is solved, we use
the obtained solution to update $\mathbf{r}^{m}$ in the upper-bound (21). We
repeat this approach until $\mathbf{r}^{m}$ converges. We summarize this
approach in Algorithm 3.

In iteration $m$, the value of the upper-bound (21) and its gradient are
$\psi(\mathbf{r}^{m})$ and $\nabla\psi(\mathbf{r}^{m})$, respectively, which
are equal to the value and the gradient of the non-separable cost function
$\psi(\mathbf{r})$. Furthermore, the upper-bound in (21) is continuous, and
thus, all four convergence conditions given in [43, Assumption 2] and listed
in the Appendix are satisfied. Due to [43, Theorem 2] and the convexity of the
non-separable cost function, the obtained solution by Algorithm 3, which
implements BSUM, is identical to the solution of the original problem with the
non-separable cost function.

######  Remark 3.

When the cost function is convex and separable, but not strictly convex, we
add a proximal term to the cost function and make it locally strongly convex
as follows:

$\displaystyle\underset{\mathbf{r}}{\min}$
$\displaystyle\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\psi_{k}^{p}(r_{k}^{p})+\frac{\kappa}{2}\left\lVert\mathbf{r}-\mathbf{r}^{m}\right\rVert_{2}^{2}$
(23) s.t. $\displaystyle\eqref{eq:linkcap},r_{k}\geq 0,\hskip
14.22636ptp\in\mathcal{P}_{k},\forall k,$

where $\kappa$ is a small positive constant. We use Algorithm 2 to solve the
above problem when we use the following equation:

$\displaystyle\frac{\partial
L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})}{\partial
r_{k}^{p}}=\alpha_{k,l}^{p,j}\frac{\partial\psi_{k}^{p}(r_{k}^{p})}{\partial
r_{k}^{p}}+\alpha_{k,l}^{p,j}\kappa(r_{k}^{p}-r_{k}^{p,m})$
$\displaystyle+\mu_{l}-\alpha_{k,l}^{p,j}\phi_{k}^{p}=0,$

instead of (15a) to find $r_{k}^{p}$, where $r_{k}^{p,m}$ is the value of
$r_{k}^{p}$ in the $m^{\text{\text{th}}}$ iteration of solving (23). We
successively solve (23) with Algorithm 2 and update $\mathbf{r}^{m}$ until
$\mathbf{r}^{m}$ converges. Similar to Remark 2, one can show that the cost
function with the proximal term in (23) satisfies the four convergence
conditions in [43, Assumption 2] and the global minimum is obtained after
successive minimizations, since each local minimum is also global for a convex
function.

##  IV Simultaneous Resource Reservations in the Backhaul and RAN

In this section, we study the joint link capacity and AP transmission resource
reservation based on the user demand and downlink statistics. Prior to the
observation of user demands, based on the formulated model in (5), the network
operator finds the optimal amount of reserved resources in the backhaul and
APs such that neither the link capacity nor AP capacity is exceeded.

###  IV-A Resource Reservation in the Backhaul

Let us drop the equality constraint (3) from (5) and substitute
$\sum_{p\in\mathcal{P}_{k}}r_{k}^{p}$ for $r_{k}$. Then, we have

$\displaystyle\underset{\mathbf{r,t}}{\min}$
$\displaystyle\sum_{k=1}^{K}\Big{[}-\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p},d_{k})]$
(24) $\displaystyle\hskip
28.45274pt+\theta_{k}\sum_{p\in\mathcal{P}_{k}}\int_{0}^{r_{k}^{p}}z_{k}^{p}(v_{k}^{p},t_{k}^{p})\>(r_{k}^{p}-v_{k}^{p})dv_{k}^{p}\Big{]}$
s.t.
$\displaystyle\eqref{eq:linkcap},\eqref{eq:nodecap},r_{k}^{p},t_{k}^{p}\geq
0,\hskip 14.22636ptp\in\mathcal{P}_{k},\forall k.$

We solve the above problem using the proposed BCD algorithm. With the fixed
$\mathbf{t}$, we minimize (24) with respect to $\mathbf{r}$ and update it.
With updated $\mathbf{r}$, we minimize (24) with respect to $\mathbf{t}$ and
update it. We underline the iterates of the BCD algorithm. In the
$i+1^{\text{th}}$ iteration of the BCD algorithm, fixing
$\underline{\mathbf{t}}^{i}$, we minimize with respect to $\mathbf{r}$. Then,
the minimization problem in (24) reduces to the following convex one:

$\displaystyle\underset{\mathbf{r}}{\min}$
$\displaystyle\sum_{k=1}^{K}\Big{[}-\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p},d_{k})]$
$\displaystyle\hskip
28.45274pt+\theta_{k}\sum_{p\in\mathcal{P}_{k}}\int_{0}^{r_{k}^{p}}z_{k}^{p}(v_{k}^{p},\underline{t}_{k}^{p,i})\>(r_{k}^{p}-v_{k}^{p})dv_{k}^{p}\Big{]}$
s.t. $\displaystyle\eqref{eq:linkcap},r_{k}^{p}\geq 0,\hskip
14.22636ptp\in\mathcal{P}_{k},\forall k.$

It is observed that although the expected outage is separable in $r_{k}^{p}$,
variables are coupled in the first term of the objective function. Therefore,
we substitute the global quadratic upper-bound given in (21) for the expected
supportable traffic demand. First, let us calculate the Lipschitz constant for
the gradient of
$-\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p},d_{k})]$. The second
derivative of $-\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p},d_{k})]$
is

$\displaystyle-\frac{\partial^{2}\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p},d_{k})]}{\partial
r_{k}^{p}\partial
r_{k}^{p^{\prime}}}=f_{k}(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p}).$

The Hessian matrix for
$-\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p},d_{k})]$ is
$|\mathcal{P}_{k}|\times|\mathcal{P}_{k}|$ dimensional, where all entries are
$f_{k}(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p})$. The eigenvalues of the Hessian
matrix are all zeros except one of them, which is
$|\mathcal{P}_{k}|f_{k}(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p})$. Therefore, the
Lipschitz constant is $|\mathcal{P}_{k}|$. We now place the Lipschitz constant
in (21) and find the upper-bound which is separable in $r_{k}^{p}$ as follows:

$\displaystyle-\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p},d_{k})]\leq-\mathbb{E}[\min(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p,m},d_{k})]$
$\displaystyle+(F_{k}(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p,m})-1)(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p}-\sum_{p\in\mathcal{P}_{k}}r_{k}^{p,m})$
$\displaystyle+\frac{|\mathcal{P}_{k}|}{2}\sum_{p\in\mathcal{P}_{k}}(r_{k}^{p,m}-r_{k}^{p})^{2},$
(25)

where $r_{k}^{p,m}$ is the $m^{\text{\text{th}}}$ iterate. We substitute
upper-bound (25) for the expected supportable demand and the optimization
problem in each iteration becomes:

$\displaystyle\underset{\mathbf{r}}{\min}$
$\displaystyle\sum_{k=1}^{K}\Bigg{[}(F_{k}(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p,m})-1)(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p}-\sum_{p\in\mathcal{P}_{k}}r_{k}^{p,m})\;$
(26)
$\displaystyle+\frac{|\mathcal{P}_{k}|}{2}\sum_{p\in\mathcal{P}_{k}}(r_{k}^{p}-r_{k}^{p,m})^{2}$
$\displaystyle+\theta_{k}\sum_{p\in\mathcal{P}_{k}}\int_{0}^{r_{k}^{p}}z_{k}^{p}(v_{k}^{p},\underline{t}_{k}^{p,i})\>(r_{k}^{p}-v_{k}^{p})dv_{k}^{p}\Bigg{]}$
s.t. $\displaystyle\eqref{eq:linkcap},r_{k}^{p}\geq 0,\hskip
14.22636ptp\in\mathcal{P}_{k},\forall k.$

We leverage Algorithm 3 to solve (IV-A) in a parallel fashion. In each
iteration of Algorithm 3, Algorithm 2 is called to solve the problem in (26).
Moreover, Algorithm 1 is called within Algorithm 2 and it needs to solve
$\frac{\partial
L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})}{\partial
r_{k}^{p}}=0$. We rewrite (15a) for the above optimization problem in the
$j^{\text{th}}$ iteration of Algorithm 2 as follows:

$\displaystyle\frac{\partial
L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})}{\partial
r_{k}^{p}}=\alpha_{k,l}^{p,j}(F_{k}(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p,m})-1)+\mu_{l}$
$\displaystyle+\alpha_{k,l}^{p,j}|\mathcal{P}_{k}|(r_{k}^{p}-r_{k}^{p,m})+\theta_{k}\>\alpha_{k,l}^{p,j}Z_{k}^{p}(r_{k}^{p},\underline{t}_{k}^{p,i})-\alpha_{k,l}^{p,j}\phi_{k}^{p}=0.$

We observe that for each $\mu_{l}$, we are able to obtain $r_{k}^{p}$
numerically using $r_{k}^{p,m}$, independent of the other variables. The
solution obtained by Algorithm 3 is unique due to the strong convexity of (26)
and is global minima as explained in Remark 2. After Algorithm 3 converges, we
set $\underline{\mathbf{r}}^{i+1}=\mathbf{r}^{m}$.

######  Remark 4.

Suppose that the number of paths that are available to user $k$ and share
downlink $w\in\mathcal{W}_{k}$ is $\varphi_{k}^{w}$. When multiple paths for
serving a user share one downlink, we substitute the quadratic upper-bound
(21) for the outage (6) as follows:

$\displaystyle\eqref{eq:mul}\leq\sum_{w\in\mathcal{W}_{k}}\int_{0}^{\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}r_{k}^{p,m}}z_{k}^{w}(v_{k}^{p},\underline{t}_{k}^{w,i})(\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}r_{k}^{p,m}-v_{k}^{p})dv_{k}^{p}$
$\displaystyle+\sum_{w\in\mathcal{W}_{k}}Z_{k}^{w}(\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}r_{k}^{p,m},\underline{t}_{k}^{w,i})(\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}(r_{k}^{p}-r_{k}^{p,m}))$
$\displaystyle+\sum_{w\in\mathcal{W}_{k}}\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}\frac{\varphi_{k}^{w}}{2}(r_{k}^{p}-r_{k}^{p,m})^{2},$ (27)

where $\varphi_{k}^{w}$ is the Lipschitz constant. In this case, the objective
function of (26) is obtained from adding (25) and the RHS of (27). We rewrite
(15a) in the $j^{\text{th}}$ iteration of Algorithm 2 for this case as
follows:

$\displaystyle\frac{\partial
L_{l}(\mathbf{r}_{l},\mu_{l},\boldsymbol{\phi}_{l},\boldsymbol{\mu}^{j-1})}{\partial
r_{k}^{p}}=\alpha_{k,l}^{p,j}(F_{k}(\sum_{p\in\mathcal{P}_{k}}r_{k}^{p,m})-1)$
$\displaystyle+\alpha_{k,l}^{p,j}|\mathcal{P}_{k}|(r_{k}^{p}-r_{k}^{p,m})+\theta_{k}\>\alpha_{k,l}^{p,j}Z_{k}^{w}(\sum_{p:\\{p\in\mathcal{P}_{k},w\in
p\\}}r_{k}^{p,m},\underline{t}_{k}^{w,i})$
$\displaystyle+\theta_{k}\varphi_{k}^{w}\alpha_{k,l}^{p,j}(r_{k}^{p}-r_{k}^{p,m})+\mu_{l}-\alpha_{k,l}^{p,j}\phi_{k}^{p}=0.$

###  IV-B Resource Reservation in RAN

When we minimize (24) with respect to $\mathbf{t}$ in the BCD algorithm, we
use $\underline{\mathbf{r}}^{i+1}$ obtained by Algorithm 3. We propose a dual
approach to minimize with respect to $\mathbf{t}$. The objective function of
(24) is separable in $t_{k}^{p}$. We are able to parallelize the algorithm
across APs since each AP has a separate transmission resource capacity
constraint. However, the problem in (24) is not necessarily convex in
$t_{k}^{p}$ for an arbitrary $z_{k}^{p}(v_{k}^{p},t_{k}^{p})$. To tackle the
potential non-convexity of the problem, we use the BSUM method and convexify
the problem locally. We iteratively solve a sequence of convex approximations.
Suppose that for each outage term in the objective function of (24), we add a
proximal term $\frac{\zeta_{k}^{p,j}}{2}\left\lVert
t_{k}^{p}-t_{k}^{p,j}\right\rVert_{2}^{2}$, $\zeta_{k}^{p,j}>0$, to make it
locally strongly convex. In the proximal term, $t_{k}^{p,j}$ is the value of
$t_{k}^{p}$ in the $j^{\text{\text{th}}}$ iteration of successively minimizing
(24) with respect to $\mathbf{t}$. The objective function with the proximal
terms is an upper-bound of the original objective function. We find the
Lagrangian for (24) with respect to $\mathbf{t}$ with proximal terms in the
objective function as follows:

$\displaystyle
L_{t}(\mathbf{t},\boldsymbol{\lambda},\boldsymbol{\beta})=\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\Bigg{(}\theta_{k}\int_{0}^{\underline{r}_{k}^{p,i+1}}z_{k}^{p}(v_{k}^{p},t_{k}^{p})\>(\underline{r}_{k}^{p,i+1}-v_{k}^{p})dv_{k}^{p}$
$\displaystyle+\frac{\theta_{k}\zeta_{k}^{p,j}}{2}\left\lVert
t_{k}^{p}-t_{k}^{p,j}\right\rVert_{2}^{2}\Bigg{)}+\sum_{b\in\mathcal{B}}\lambda_{b}(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}t_{k}^{p}-C_{b})$
$\displaystyle-\sum_{k=1}^{K}\sum_{p\in\mathcal{P}_{k}}\beta_{k}^{p}t_{k}^{p},$

where $\underline{\mathbf{r}}^{i+1}$ block is kept fixed. We can decompose the
above Lagrangian across APs as follows:

$\displaystyle L_{t,b}(\mathbf{t}_{b},\lambda_{b},\boldsymbol{\beta}_{b})=$
$\displaystyle\theta_{k}\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}\Bigg{(}\int_{0}^{\underline{r}_{k}^{p,i+1}}z_{k}^{p}(v_{k}^{p},t_{k}^{p})\>(\underline{r}_{k}^{p,i+1}-v_{k}^{p})dv_{k}^{p}$
$\displaystyle+\frac{\zeta_{k}^{p,j}}{2}\left\lVert
t_{k}^{p}-t_{k}^{p,j}\right\rVert_{2}^{2}\Bigg{)}$
$\displaystyle+\lambda_{b}\>(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}t_{k}^{p}-C_{b})-\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}\beta_{k}^{p}t_{k}^{p},$
(28)

where
$\mathbf{t}_{b}=\\{t_{k}^{p}\\}_{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p},k=1:K}$
and
$\boldsymbol{\beta}_{b}=\\{\beta_{k}^{p}\\}_{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p},k=1:K}\geq\mathbf{0}$.
To develop an algorithm to solve each subproblem with respect to
$\mathbf{t}_{b}$, we use KKT conditions. We write the first-order optimality
conditions with respect to $\mathbf{t}$ as follows:

$\displaystyle\frac{\partial
L_{t,b}(\mathbf{t}_{b},\lambda_{b},\boldsymbol{\beta}_{b})}{\partial
t_{k}^{p}}=\theta_{k}\int_{0}^{\underline{r}_{k}^{p,i+1}}\frac{\partial
z_{k}^{p}(v_{k}^{p},t_{k}^{p})}{\partial
t_{k}^{p}}\>(\underline{r}_{k}^{p,i+1}-v_{k}^{p})dv_{k}^{p}$
$\displaystyle+\lambda_{b}+\theta_{k}\zeta_{k}^{p,j}\>(t_{k}^{p}-t_{k}^{p,j})-\beta_{k}^{p}=0,$
(29a)
$\displaystyle\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}t_{k}^{p}\leq
C_{b},$ (29b)
$\displaystyle\lambda_{b}(\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}t_{k}^{p}-C_{b})=0,\>\>\lambda_{b}\geq
0,$ (29c) $\displaystyle\beta_{k}^{p}t_{k}^{p}=0,\>\>t_{k}^{p}\geq
0,\>\>\beta_{k}^{p}\geq 0.$ (29d)

From (29a), we observe that a given dual variable $\lambda_{b}$, which
corresponds to AP $b$, identifies the reserved resource $t_{k}^{p}$ for all
downlinks created by that AP. The proposed dual algorithm works as follows:
implement a bisection search on $\lambda_{b}$ in the non-negative orthant and
find each $t_{k}^{p}:t_{k}^{p}\geq 0$, which is associated with the AP $b$,
from (29a) when $\beta_{k}^{p}=0$. Continue the bisection search until one
$\lambda_{b}$ is obtained such that for the obtained $\lambda_{b}$, we have
$\sum_{k=1}^{K}\sum_{p:\\{p\in\mathcal{P}_{k},b\in\mathcal{U}_{k}^{p}\\}}t_{k}^{p}=C_{b}$.
If there is no such $\lambda_{b}$, we set $\lambda_{b}=0$ and solve (29a) and
(29d) without (29b)–(29c). Once the optimized variables are obtained, we
update $t_{k}^{p,j}$ and $j=j+1$. We repeat the same process until
$\mathbf{t}^{j}=\\{t_{k}^{p,j}\\}_{p\in\mathcal{P}_{k},k=1:K}$ converges. As
it is explained in Remark 3, after a sequence of upper-bound minimizations and
updating the proximal terms in the objective function, a KKT (local
stationary) solution to the original problem is obtained. If expected outage
terms for downlinks are non-increasing in $\mathbf{t}$, one can show that the
successive upper-bound minimization converges to the global minima with
respect to $\mathbf{t}$. After $\mathbf{t}^{j}$ converges, we set
$\underline{\mathbf{t}}^{i+1}=\mathbf{t}^{j}$.

###  IV-C The Proposed BCD Algorithm

To solve the problem in (24) to a KKT point, we optimize with respect to two
blocks of variables, $\mathbf{r}$ and $\mathbf{t}$, alternatively with the
Gauss-Seidel update style. Therefore, if we choose $\mathbf{r}$ to update
first, with $\underline{\mathbf{r}}^{i+1}$, we optimize with respect to
$\mathbf{t}$, and then, we update $\underline{\mathbf{t}}^{i+1}$. We keep
optimizing with respect to $\mathbf{r}$ and $\mathbf{t}$ alternatively until
both blocks converge. The summary of the overall BCD approach is given in
Algorithm 4.

0\. Initialization Feasible initializations for $\underline{\mathbf{r}}^{0}$
and $\underline{\mathbf{t}}^{0}$, $i=0$;

repeat

1\. Apply Algorithm 3 to solve (24) and find $\underline{\mathbf{r}}^{i+1}$;

2\. Solve (24) with respect to $\mathbf{t}$ and find
$\underline{\mathbf{t}}^{i+1}$;

3\. $i=i+1$;

until
_$\left\lVert\underline{\mathbf{r}}^{i}-\underline{\mathbf{r}}^{i-1}\right\rVert_{2}^{2}+\left\lVert\underline{\mathbf{t}}^{i}-\underline{\mathbf{t}}^{i-1}\right\rVert_{2}^{2}$
is small enough_;

Algorithm 4 The proposed BCD algorithm to solve (24)

######  Proposition 3.

Algorithm 4 converges to a KKT solution to (24).

###### Proof.

First, the objective function of (24) is continuously differentiable. Second,
feasible sets of two blocks of variables are separate in (24). Hence, updating
one block of variables does not change the other block. Third, in each
iteration of Algorithm 4, a KKT solution is obtained. Therefore, according to
[45, Proposition 3.7.1], the proposed (BCD) Algorithm 4 converges to a KKT
solution. ∎

##  V Numerical Tests

In this section, we demonstrate the performance of our proposed approach
against two heuristic algorithms.

###  V-A Simulation Setup

The considered network for evaluations is shown in Fig. 2, which includes both
the backhaul and radio access parts. A data center is connected to routers of
the network through three gateway routers, GW $1$, GW $2$, and GW $3$. The
network includes $57$ APs and $11$ network routers. APs are distributed on the
X-Y plane and they are connected to each other and routers via wired links.
The backhaul network has $162$ links. Wired link capacities are identical in
both directions. Backhaul link capacities are determined as

  * •

Links between the data center and routers: $4$ Gnats/s;

  * •

Links between routers: 2 Gnats/s;

  * •

Links between routers and APs: $2$ Gnats/s;

  * •

2-hop to the routers: $400$ Mnats/s;

  * •

3-hop to the routers: $320$ Mnats/s;

  * •

4-hop to the routers: $160$ Mnats/s.

The considered paths originate from the data center and are extended toward
users. We consider $200$ users are distributed randomly in the same plane of
APs; however, they are not shown in Fig. 2. User AP associations are
determined by the highest long-term received power. We consider three wireless
connections, which have the highest received power, to serve each user. There
are three paths for carrying data from a data center to APs. The distribution
of the demand is log-normal:

$\displaystyle d_{k}\sim\frac{1}{d_{k}\sigma_{k}\sqrt{2\pi}}\exp(-\frac{(\ln
d_{k}-\eta_{k})^{2}}{2\sigma_{k}^{2}}).$ (30)

In addition, it is assumed that $\eta_{k}$ is realized randomly from a normal
distribution for each user. The power allocations in APs are fixed. The
dispensed resource in an AP is bandwidth. The channel between each user and an
AP is a Rayleigh fading channel. The CDF of the wireless channel capacity,
which is parameterized by the allocated bandwidth $t_{k}^{p}$, is given as
follows [35]:

$\displaystyle
Z_{k}^{p}(v_{k}^{p},t_{k}^{p})=1-\exp(\frac{1-2^{v_{k}^{p}/t_{k}^{p}}}{\overline{\text{SNR}_{k}^{p}}}),$

where $\overline{\text{SNR}_{k}^{p}}$ is the average SNR. The PDF of the
wireless channel capacity is

$\displaystyle
z_{k}^{p}(v_{k}^{p},t_{k}^{p})=\frac{\ln(2)2^{v_{k}^{p}/t_{k}^{p}}\exp(\frac{1-2^{v_{k}^{p}/t_{k}^{p}}}{\overline{\text{SNR}_{k}^{p}}})}{\overline{\text{SNR}_{k}^{p}}t_{k}^{p}}.$
(31)

Benchmark heuristic algorithms are the single-path and the average-based
approaches. In the single-path approach, each user is served through one path
from a data center to a user. Moreover, the average-based algorithm only
considers the mean of the user demand and the average achievable rate of a
downlink. To compare algorithms, with an identical network, we measure the
objective function of (5), the sum of user expected supportable rates, the
aggregate expected outage of downlinks and the amount of traffic that each
algorithm can reserve for users. One datastream is associated with each user.
In total, we have $600$ paths in the backhaul. We use C to implement
algorithms.

Figure 2: A wireless data network consists of $57$ APs and $11$ routers.

###  V-B Learning Probability Density Functions

The optimization problem in (5) takes into account PDFs of user demands and
achievable rates of downlinks. When PDFs are not given, one can use a data-
driven approach to learn PDFs used in (5) based on collected observations.
Upon the collection of user demands and achievable rates of downlinks, one can
estimate the PDFs using a recursive non-parametric estimator. In order to
estimate PDFs in an online streaming fashion, one can use efficient recursive
kernel estimators, such as the Wolverton and Wagner estimator [47]. Suppose
that independent random variables $X_{1},X_{2},\dots,X_{n}$ are observations
that are collected from an identical PDF $\chi$ with respect to Lebesgue’s
measure. The estimated PDF is

$\displaystyle\hat{\chi}_{n,\mathbf{h}_{n}}=\frac{1}{n}\sum_{k=1}^{n}\frac{1}{h_{k}}K(\frac{X_{k}-x}{h_{k}}),$

where $\mathbf{h}_{n}=(h_{1},h_{2},\dots,h_{n})$, $h_{1}>\dots>h_{n}$ and
$K(\cdot)$ is a kernel function. The advantage of the above estimator is that
it can be written in a recursive form as follows:

$\displaystyle\hat{\chi}_{n+1,\mathbf{h}_{n+1}}=\frac{n}{n+1}\hat{\chi}_{n,\mathbf{h}_{n}}+\frac{1}{(n+1)h_{n+1}}K(\frac{X_{n+1}-x}{h_{n+1}}),$

which makes it suitable for real-time applications. The bandwidth selection in
[48] can be used for the above estimator. The bandwidth $h_{k}$ is selected in
[48] as $h_{k}=k^{-\gamma},k\in\\{1,\dots,n\\}$, where
$\gamma=\frac{1}{2\beta+1}$ and $\beta>0$.

###  V-C Simulation Results

Before demonstrating the performance of Algorithm 4, we depict the convergence
of Algorithm 2 in Fig. 3. The convergence of Algorithm 2 for different means
of the user demand is depicted in Fig. 3. It is observed that Algorithm 2 has
a fast convergence for the large network of Fig. 2 with $600$ paths. Numerical
results show that the number of required iterations for Algorithm 3 to
converge for the simulation setting described above is at most $60$. The CPU
time for Algorithm 3 is measured and is given in Table I.

Figure 3: The convergence of Algorithm 2. Table I: CPU time for $60$ iterations of Algorithm 3. Mean of $\eta_{k}$ |  $1$ Mnats/s |  $2$ Mnats/s |  $3$ Mnats/s |  $4$ Mnats/s |
---|---|---|---|---|---
CPU time |  $0.051$ s |  $0.082$ s |  $0.104$ s |  $0.121$ s |
Mean of $\eta_{k}$ |  |  $5$ Mnats/s |  $6$ Mnats/s |  |
CPU time |  |  $0.144$ s |  $0.158$ s |  |

(a)

(b)

(c)

Figure 4: (a) Reserved rates by Algorithm 4 (multi-path) and the single-path
approach when wireless channels are deterministic. (b) The expected
supportable rates for users by Algorithm 4 and the single-path approach when
wireless channels are deterministic. (c) The objective function of problem (5)
with the single-path approach and Algorithm 4 when wireless channels are
stochastic.

(a)

(b)

Figure 5: Stochastic wireless channels: performance of the single-path
approach and Algorithm 4 in terms of (a) the aggregate expected supportable
traffic; and (b) aggregate reserved rates.

(a)

(b)

Figure 6:  Stochastic wireless channels: performance of the single-path
approach and Algorithm 4 in terms of (a) expected outage of downlinks; and (b)
reserved bandwidth.

First, let us assume that transmission rates on downlinks are deterministic
functions of bandwidth in APs. Therefore, no outage (rate loss) is considered.
For each downlink, the transmission rate and the allocated bandwidth are
connected to each other as $r_{k}^{p}=\delta_{k}^{p}t_{k}^{p}$, where
$\delta_{k}^{p}$ is the spectral efficiency of the downlink of path $p$ to
serve user $k$. Furthermore, suppose that $\sigma_{k}=3.6$ and the capacity of
each backhaul link listed previously is divided by $4$. When the bandwidth
budget of each AP increases from $15$ MHz to $40$ MHz, the aggregate reserved
rates for users by Algorithm 4 (multi-path) and the single-path approach are
shown in Fig. 4(a). The aggregate expected supportable rates of users with
both approaches are depicted in Fig. 4(b). It is observed that Algorithm 4
outperforms the single-path approach. Both approaches utilize all available
bandwidth in APs.

Consider the distribution of each wireless channel (downlink) achievable rate
follows (31) and backhaul link capacities are as listed previously. Suppose
that the available bandwidth in each AP increases by a step size of $10$ MHz,
where $\theta_{k}=1/2$ and $\sigma_{k}=0.6$. The objective function of the
problem in (5) by Algorithm 4 and the single-path approach are compared in
Fig. 4(c). Our proposed Algorithm 4 outperforms the single-path approach. It
is observed that with the increase of mean for $\eta_{k}$ and the AP bandwidth
budget, the objective function increases.

The expected supportable demands of users, depicted in Fig. 5(a), increases
when the mean of $\eta_{k}$ and the AP bandwidth budget increase. It is
observed from Fig. 5(a) that the aggregate expected supportable traffic for
users obtained by Algorithm 4 is greater than that by the single-path
approach. In Fig. 5(b), we observe that the aggregate reserved rates for users
increases with the increase of mean for $\eta_{k}$. Furthermore, it increases
when the bandwidth budgets of APs increase. From Fig. 6(a), we observe that
the aggregate expected outage increases as the mean of $\eta_{k}$ increases
and decreases when the AP bandwidth budget increases. We observe from Fig.
6(b) that the bandwidth reservation by Algorithm 4 is almost equal to that by
the single-path approach. Numerical results show that $5$ iterations are
sufficient for the convergence of Algorithm 4.

Next, we evaluate the performance of Algorithm 4 against the average-based
approach when both the demand and downlink achievable rates are stochastic.
The average-based algorithm is oblivious to the user demand and the downlink
achievable rate distributions. It only considers the average of each user
demand and the average achievable rate of a downlink. The average-based
approach uses the same set of paths used by Algorithm 4. The bandwidth budget
in each AP is 40 MHz. Furthermore, $\sigma_{k}=0.6$ and $\theta_{k}=1/3$. The
demand and downlink achievable rate distributions are as given in (30) and
(31), respectively. Both approaches are set to make reservations for users
assuming the mean of $\eta_{k}$ is $2$ Mnast/s. We generate $100$ scenarios in
which user demands and downlink capacities are random. For each scenario, we
measure how much the user demands are satisfied using the reserved resources
in the network by both approaches. After collecting results for $100$
scenarios, we plot the empirical CDF for the supply demand ratio in Fig. 7. It
is observed that when the mean of demand exceeds what it was supposed to be,
the resource reservation made by Algorithm 4 is more robust and supports
random demands better. The total reserved link capacities in the backhaul by
Algorithm 4 is $1.0979\times 10^{4}$ Mnats/s and is $7.74\times 10^{3}$
Mnats/s by the average-based approach. Furthermore, the total reserved
bandwidth in RAN by Algorithm 4 is $2.043\times 10^{3}$ MHz and is
$1.968\times 10^{3}$ MHz by the average-based approach.

Figure 7:  The probability of being able to support the user demands up to a
certain percentage.

##  VI Concluding Remarks and Future Directions

In this paper, we studied link capacity and transmission resource reservation
in wireless data networks prior to the observation of user demands. Using the
statistics of user demands and achievable rates of downlinks, we formulated an
optimization problem to maximize the sum of user expected supportable traffic
while minimizing the expected outage of downlinks. We demonstrated that this
problem is non-convex in general. To solve the problem approximately, an
efficient BCD approach is proposed which benefits from distributed and
parallel computation when each block of variables is chosen to be updated. We
demonstrated that despite the non-convexity of the problem, our proposed
approach converges to a KKT solution to the problem. We verified the
efficiency and the efficacy of our proposed approach against two heuristic
algorithms developed for joint resource reservation in the backhaul and RAN.

In future work, we consider multi-tenant networks and reservation-based
network slicing. In addition to users, tenants have different requirements
[49], and maximum isolation between sliced resources should be enforced [50].
The demand distribution of users may change over time and the network
resources should be sliced for tenants accordingly. However, the slice
reconfiguration for each tenant involves cost and overhead. Based on the cost
of reconfiguration and newly arrived statistics, we formulate the problem from
a sparse optimization perspective and propose an efficient approach based on
iteratively solving a sequence of group Least Absolute Shrinkage and Selection
Operator (LASSO) problems [49].

## References

  * [1] N. Reyhanian, H. Farmanbar, and Z.-Q. Luo,  “Resource reservation in backhaul and radio access network with uncertain user demands,”  in Proc. IEEE Signal Process. Adv. Wireless Commun. (SPAWC), May 2020, pp. 1–5.
  * [2] S. Albasheir and M. Kadoch,  “Enhanced control for adaptive resource reservation of guaranteed services in LTE networks,”  IEEE Internet Things J., vol. 3, no. 2, pp. 179–189, Apr. 2015\.
  * [3] K. Kaur, A. Dua, A. Jindal, N. Kumar, M. Singh, and A. Vinel,  “A novel resource reservation scheme for mobile PHEVs in V2G environment using game theoretical approach,”  IEEE Trans. Veh. Technol., vol. 64, no. 12, pp. 5653–5666, Dec. 2015.
  * [4] E. Van Den Berg, T. Zhang, J. Chennikara, P. Agrawal, and T. Kodama,  “Time series-based localized predictive resource reservation for handoff in multimedia wireless networks,”  in Proc. IEEE Int. Conf. Commun., 2001, vol. 2, pp. 346–350.
  * [5] R. E. Gomory and T. C. Hu,  “An application of generalized linear programming to network flows,”  J. Soc. Ind. Appl. Math., vol. 10, no. 2, pp. 260–283, 1962.
  * [6] R. Dai, L. Li, S. Wang, and X. Zhang,  “Planning traffic-oblivious survivable WDM networks using differentiated reliable partial SRLG-disjoint protection,”  in Proc. Sym. Photon. Optoelectronics, Jun. 2010, pp. 1–7.
  * [7] P. Kumar, Y. Yuan, C. Yu, N. Foster, R. Kleinberg, and R. Soulé,  “Kulfi: Robust traffic engineering using semi-oblivious routing,”  arXiv preprint arXiv:1603.01203, 2016.
  * [8] C. Cicconetti, V. Gardellin, L. Lenzini, E. Mingozzi, and A. Erta,  “End-to-end bandwidth reservation in IEEE 802.16 mesh networks,”  in Proc. IEEE Int. Conf. Mobile Adhoc and Sensor Syst., Oct. 2007, pp. 1–6.
  * [9] D. Applegate and E. Cohen,  “Making routing robust to changing traffic demands: algorithms and evaluation,”  IEEE/ACM Trans. Netw., vol. 14, no. 6, pp. 1193–1206, Dec. 2006\.
  * [10] N. Moehle, X. Shen, Z.-Q. Luo, and S. Boyd,  “A distributed method for optimal capacity reservation,”  J. Optim. Theory Appl., vol. 182, no. 3, pp. 1130–1149, May 2019\.
  * [11] D. Ma, B. Sheng, S. Jin, X. Ma, and P. Gao,  “Short-term traffic flow forecasting by selecting appropriate predictions based on pattern matching,”  IEEE Access, vol. 6, pp. 75629–75638, Nov. 2018.
  * [12] M. Yan, G. Feng, J. Zhou, Y. Sun, and Y.-C. Liang,  “Intelligent resource scheduling for 5G radio access network slicing,”  IEEE Trans. Veh. Technol., vol. 68, no. 8, pp. 7691–7703, Jun. 2019\.
  * [13] Q. He, A. Moayyedi, G. Dán, G. P. Koudouridis, and P. Tengkvist,  “A meta-learning scheme for adaptive short-term network traffic prediction,”  IEEE J. Sel. Areas Commun., Jun. 2020.
  * [14] L. U. Khan, I. Yaqoob, N.-H. Tran, Z. Han, and C.-S. Hong,  “Network slicing: Recent advances, taxonomy, requirements, and open research challenges,”  IEEE Access, vol. 8, pp. 36009–36028, Feb. 2020.
  * [15] J. Prados-Garzon, A. Laghrissi, M. Bagaa, T. Taleb, and J. M. Lopez-Soler,  “A complete LTE mathematical framework for the network slice planning of the EPC,”  IEEE Trans. Mobile Comput., vol. 19, no. 1, pp. 1–14, Jan. 2020\.
  * [16] Y. Li, J. Liu, B. Cao, and C. Wang,  “Joint optimization of radio and virtual machine resources with uncertain user demands in mobile cloud computing,”  IEEE Trans. Multimedia, vol. 20, no. 9, pp. 2427–2438, Sep. 2018\.
  * [17] T. Hößler, P. Schulz, E. A. Jorswieck, M. Simsek, and G. P. Fettweis,  “Stable matching for wireless URLLC in multi-cellular, multi-user systems,”  IEEE Trans. Commun., vol. 68, no. 8, pp. 5228–5241, Aug. 2020.
  * [18] D. P. Bertsekas,  Linear network optimization: algorithms and codes,  MIT press, 1991.
  * [19] J. Tsitsiklis and D. Bertsekas,  “Distributed asynchronous optimal routing in data networks,”  IEEE Trans. Autom. Control, vol. 31, no. 4, pp. 325–332, Apr. 1986\.
  * [20] D. P. Bertsekas,  Network optimization: continuous and discrete models,  Athena Scientific Belmont, MA, 1998.
  * [21] H. Zhang and V. W. S. Wong,  “A two-timescale approach for network slicing in C-RAN,”  IEEE Trans. Veh. Technol., vol. 69, no. 6, pp. 6656–6669, Jun. 2020\.
  * [22] L. Li and A. J. Goldsmith,  “Capacity and optimal resource allocation for fading broadcast channels. II. outage capacity,”  IEEE Trans. Inf. Theory, vol. 47, no. 3, pp. 1103–1127, Mar. 2001\.
  * [23] X. Liao, J. Shi, Z. Li, L. Zhang, and B. Xia,  “A model-driven deep reinforcement learning heuristic algorithm for resource allocation in ultra-dense cellular networks,”  IEEE Trans. Veh. Technol., vol. 69, no. 1, pp. 983–997, Jan. 2019\.
  * [24] V. Sciancalepore, X. Costa-Perez, and A. Banchs,  “RL-NSB: Reinforcement learning-based 5G network slice broker,”  IEEE/ACM Trans. Netw., vol. 27, no. 4, pp. 1543–1557, Aug. 2019\.
  * [25] Y. Liang and V. V. Veeravalli,  “Gaussian orthogonal relay channels: Optimal resource allocation and capacity,”  IEEE Trans. Inf. Theory, vol. 51, no. 9, pp. 3284–3289, Sep. 2005\.
  * [26] S.-J. Kim and G. B. Giannakis,  “Optimal resource allocation for MIMO ad–hoc cognitive radio networks,”  IEEE Trans. Inf. Theory, vol. 57, no. 5, pp. 3117–3131, May 2011\.
  * [27] W.-C. Liao, M. Hong, Y.-F. Liu, and Z.-Q. Luo,  “Base station activation and linear transceiver design for optimal resource management in heterogeneous networks,”  IEEE Trans. Signal Process., vol. 62, no. 15, pp. 3939–3952, Jul. 2014.
  * [28] L. Xiao, M. Johansson, and S. P. Boyd,  “Simultaneous routing and resource allocation via dual decomposition,”  IEEE Trans. Commun., vol. 52, no. 7, pp. 1136–1144, Jul. 2004\.
  * [29] A. A. El-Sherif and A. Mohamed,  “Joint routing and resource allocation for delay minimization in cognitive radio based mesh networks,”  IEEE Trans. Wireless Commun., vol. 13, no. 1, pp. 186–197, Jan. 2014.
  * [30] K. Wang, K. Yang, and C. S. Magurawalage,  “Joint energy minimization and resource allocation in C-RAN with mobile cloud,”  IEEE Trans. Cloud Comput., vol. 6, no. 3, pp. 760–770, Sep. 2018\.
  * [31] S. Matoussi, I. Fajjari, S. Costanzo, N. Aitsaadi, and R. Langar,  “5G RAN: Functional split orchestration optimization,”  IEEE J. Sel. Areas Commun., vol. 38, no. 7, pp. 1448–1463, Jul. 2020.
  * [32] J. Liu, Y. Pang, H. Ding, L. Cai, H. Zhang, and Y. Fang,  “Optimizing IoT Energy Efficiency on Edge (EEE): a cross-layer design in a cognitive mesh network,”  arXiv preprint arXiv:1901.05494, 2019.
  * [33] H. Kordbacheh, H. Dalili Oskouei, and N. Mokari,  “Robust cross-layer routing and radio resource allocation in massive multiple antenna and OFDMA-based wireless ad-hoc networks,”  IEEE Access, vol. 7, pp. 36527–36539, Mar. 2019.
  * [34] K. Karakayali, J. H. Kang, M. Kodialam, and K. Balachandran,  “Joint resource allocation and routing for OFDMA-based broadband wireless mesh networks,”  in Proc. IEEE Int. Conf. Commun. (ICC), Jun. 2007, pp. 5088–5092.
  * [35] S. Choudhury and J. D. Gibson,  “Information transmission over fading channels,”  in Proc. IEEE Global Commun. Conf., Nov. 2007, pp. 3316–3321.
  * [36] D. Wu and R. Negi,  “Effective capacity: a wireless link model for support of quality of service,”  IEEE Trans. Wireless Commun., vol. 2, no. 4, pp. 630–643, Jul. 2003\.
  * [37] O. Ertug,  “Asymptotic ergodic capacity of multidimensional vector-sensor array MIMO channels,”  IEEE Trans. Wireless Commun., vol. 7, no. 9, pp. 3297–3300, Sep. 2008.
  * [38] W.-C. Liao, M. Hong, H. Farmanbar, and Z.-Q. Luo,  “A distributed semiasynchronous algorithm for network traffic engineering,”  IEEE Trans. Signal Inf. Process. Netw., vol. 4, no. 3, pp. 436–450, Sep. 2018.
  * [39] H. K. Nguyen, Y. Zhang, Z. Chang, and Z. Han,  “Parallel and distributed resource allocation with minimum traffic disruption for network virtualization,”  IEEE Trans. Commun., vol. 65, no. 3, pp. 1162–1175, Mar. 2017\.
  * [40] Z. Wu, Z. Fei, Y. Yu, and Z. Han,  “Toward optimal remote radio head activation, user association, and power allocation in C-RANs using benders decomposition and ADMM,”  IEEE Trans. Commun., vol. 67, no. 7, pp. 5008–5023, Jul. 2019.
  * [41] D. P. Palomar and M. Chiang,  “A tutorial on decomposition methods for network utility maximization,”  IEEE J. Sel. Areas Commun., vol. 24, no. 8, pp. 1439–1451, Aug. 2006.
  * [42] P. Kemmer, A. K. Strauss, and T. Winter,  “Dynamic simultaneous fare proration for large-scale network revenue management,”  J. Oper. Res. Soc., vol. 63, no. 10, pp. 1336–1350, 2012.
  * [43] M. Razaviyayn, M. Hong, and Z.-Q. Luo,  “A unified convergence analysis of block successive minimization methods for non-smooth optimization,”  SIAM J. Optim., vol. 23, no. 2, pp. 1126–1153, 2013.
  * [44] S. Boyd and L. Vandenberghe,  Convex optimization,  Cambridge university press, 2004.
  * [45] D. P. Bertsekas,  Nonlinear Programming,  Athena Scientific, 3rd ed., 2016.
  * [46] M. Hong, M. Razaviyayn, Z.-Q. Luo, and J.-S. Pang,  “A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing,”  IEEE Signal Process. Mag., vol. 33, no. 1, pp. 57–77, Dec. 2015\.
  * [47] C. Wolverton and T. V. Wagner,  “Asymptotically optimal discriminant functions for pattern classification,”  IEEE Trans. Inf. Theory, vol. 15, no. 2, pp. 258–265, Mar. 1969\.
  * [48] F. Comte and N. Marie,  “Bandwidth selection for the Wolverton–Wagner estimator,”  J. Stat. Planning Inference, vol. 207, pp. 198–214, 2020.
  * [49] N. Reyhanian, H. Farmanbar, and Z.-Q. Luo,  “Data-driven adaptive network resource slicing for multi-tenant networks,”  in Proc. IEEE Int. Conf. Acoust. Speech Signal Process (ICASSP), Jun. 2021.
  * [50] N. Reyhanian and B. Maham,  “Statistical slice selection in multi-tenant networks with maximum isolation of reserved resources,”  in Proc. 54th Asilomar Conf. Signals, Syst. Comput., Pacific Grove, CA, Nov. 2020.

[Block Successive Upper-Bound Minimization] Notations in this Appendix are
identical to [43] and are not related to those defined in the paper. According
to the BSUM algorithm [43, Theorem 2], when an upper-bound satisfies four
conditions, the solution acquired by the BSUM converges to a local minima to
the problem. Here, we give a brief description of the BSUM approach. Suppose
that $u_{i}(\mathbf{x},\mathbf{x}^{t-1})$ is an upper-bound for an arbitrary
objective function $f(\mathbf{x})$ at the point $\mathbf{x}^{t-1}$. In
iteration $t$, one selected block (say, block i) is optimized by solving the
following subproblem:

$\displaystyle\underset{\mathbf{x}_{i}}{\min}$ $\displaystyle
u_{i}(\mathbf{x}_{i},\mathbf{x}^{t-1})$ (32) s.t.
$\displaystyle\mathbf{x}_{i}\in\mathcal{X}_{i},$

where $\mathcal{X}_{i}$ is the feasible set of block $\mathbf{x}_{i}$.
Conditions on the upper-bound are listed in [43, Assumption 2] as follows:

  1. 1.

$u_{i}(\mathbf{y}_{i},\mathbf{y})=f(\mathbf{y}),\leavevmode\nobreak\
\forall\mathbf{y}\in\mathcal{X},\forall i,$

  2. 2.

$u_{i}(\mathbf{x}_{i},\mathbf{y})\geq
f(\mathbf{y}_{1},\dots,\mathbf{y}_{i-1},\mathbf{x}_{i},\mathbf{y}_{i+1},\dots,\mathbf{y}_{n}),\leavevmode\nobreak\
\forall\mathbf{x}_{i}\in\mathcal{X}_{i},\forall\mathbf{y}\in\mathcal{X},\forall
i,$

  3. 3.

$u_{i}^{{}^{\prime}}(\mathbf{x}_{i},\mathbf{y};\mathbf{d}_{i})|_{\mathbf{x}_{i}=\mathbf{y}_{i}}=f^{\prime}(\mathbf{y};\mathbf{d}),\>\forall\mathbf{d}=(0,\dots,\mathbf{d}_{i},\dots,0)\>\>\text{s.t.}\>\>\mathbf{y}_{i}+\mathbf{d}_{i}\in\mathcal{X}_{i},\forall
i,$

  4. 4.

$u_{i}(\mathbf{x}_{i},\mathbf{y})$ is continuous in
$(\mathbf{x}_{i},\mathbf{y}),\forall i$.

When problem (32) is solved sequentially for different $i$ and there exists a
unique solution for each subproblem, $\mathbf{x}^{t}$ converges to a KKT point
of $f(\mathbf{x})$.