# Near-Optimal Packet Scheduling in Multihop Networks with End-to-End Deadline
Constraints

Christos Tsanikidis<EMAIL_ADDRESS>Columbia UniversityNew
YorkNYUSA and  Javad Ghaderi  Columbia UniversityNew YorkNYUSA
<EMAIL_ADDRESS>
###### Abstract.

Scheduling packets with end-to-end deadline constraints in multihop networks
is an important problem that has been notoriously difficult to tackle.
Recently, there has been progress on this problem in the worst-case traffic
setting, with the objective of maximizing the number of packets delivered
within their deadlines. Specifically, the proposed algorithms were shown to
achieve $\Omega(1/\log(L))$ fraction of the optimal objective value if the
minimum link capacity in the network is $\mathrm{C}_{\min}=\Omega(\log(L))$,
where $L$ is the maximum length of a packet’s route in the network (which is
bounded by the packet’s maximum deadline). However, such guarantees can be
quite pessimistic due to the strict worst-case traffic assumption and may not
accurately reflect real-world settings. In this work, we aim to address this
limitation by exploring whether it is possible to design algorithms that
achieve a constant fraction of the optimal value while relaxing the worst-case
traffic assumption. We provide a positive answer by demonstrating that in
stochastic traffic settings, such as i.i.d. packet arrivals, near-optimal,
$(1-\epsilon)$-approximation algorithms can be designed if
$\mathrm{C}_{\min}=\Omega\big{(}\frac{\log(L/\epsilon)}{\epsilon^{2}}\big{)}$.
To the best of our knowledge, this is the first result that shows this problem
can be solved near-optimally under nontrivial assumptions on traffic and link
capacity. We further present extended simulations using real network traces
with non-stationary traffic, which demonstrate that our algorithms outperform
worst-case-based algorithms in practical settings.

Scheduling algorithms; Approximation algorithms; Multihop traffic; Deadline
scheduling

††conference: ; ;

##  1\. Introduction

In recent years, the scheduling of real-time traffic in communication networks
has become increasingly important, due to growing number of emerging real-time
applications such as video streaming and video conferencing, vehicular
networks, cyber-physical networks, and Internet-of-Things (Popovski et al.,
2022). The connectivity within these networks is often complex, with packets
generated at specific sources requiring the traversal of multiple links in
order to reach their destinations successfully. Additionally, in many real-
time applications, the timely delivery of packets is a primary concern, as
packets that fail to meet specific deadlines, on the total time from
generation of a packet at its source until delivery to its destination, are
typically discarded by the application. Meeting the deadline constraints
requires a departure from traditional scheduling algorithms (e.g. based on
MaxWeight or Backpressure (Tassiulas and Ephremides, 1992)) that have been
designed for maximizing throughput and cannot provide deadline guarantees on
packet delivery.

Despite the importance of the problem due to its broad applicability, there is
very limited work on techniques with attractive theoretical guarantees in
multihop networks. This is due to the complexity of the problem, including the
need to make online decisions, the exponential growth in the number of
scheduling decisions that involve the path a packet takes as well as the
specific time slots the packet occupies for transmission on the links, and the
stringent deadline constraints. In particular, prior literature on scheduling
packets with deadlines in multihop networks has focused on either worst-case
traffic (Tsanikidis and Ghaderi, 2022; Gu et al., 2021; Deng et al., 2019; Mao
et al., 2014) with pessimistic approximation ratios, or stochastic traffic
with guarantees in the case that the link capacities are relaxed to be only
satisfied on average (Singh and Kumar, 2018), or in an asymptotic regime where
the link capacities and packet arrival rates scale to infinity (Singh and
Kumar, 2018). As a result, an important open question remains: is it possible
to design algorithms that have attractive performance guarantees (e.g.
constant approximation ratio that does not depend on parameters of the
network), for finite link capacities and given stochastic traffic? In this
work, we answer this question positively by providing algorithms that are
near-optimal while only requiring minimum link capacities that are logarithmic
in the maximum length of a packet’s route $L$ (with $L\leq d_{max}$, where
$d_{max}$ is the maximum deadline of any packet). Further, our algorithms have
polynomial computational complexity, are amenable to distributed
implementations, and can be easily adapted to different traffic distribution
assumptions. Additionally, our work features distinct techniques compared to
prior work in the field (Gu et al., 2021; Deng et al., 2019; Singh and Kumar,
2018; Tsanikidis and Ghaderi, 2022).

The near-optimality of our algorithms is with respect to the commonly studied
average approximation ratio, which is the fraction of the optimal objective
value our algorithm obtains on average (this modifies the worst-case
approximation ratio in (Tsanikidis and Ghaderi, 2022; Deng et al., 2019; Gu et
al., 2021) to a stochastic setting)111In the literature, similar results are
presented in terms of the competitiveness, which is the reciprocal of the
approximation ratio, that is to say, $a$-approximation ratio corresponds to a
$1/a$-competitiveness.. In the weighted-packet case, unlike prior techniques,
e.g. (Tsanikidis and Ghaderi, 2022; Gu et al., 2021), our results do not
depend on the weights of the packets. Additionally, our results hold for
finite time horizons, as opposed to prior works in the stochastic setting that
require infinite time horizons, e.g. (Singh and Kumar, 2018; Singh et al.,
2013).

In our view, this work fills an important gap in the literature regarding the
feasibility of solving the real-time scheduling problem near-optimally in
practical multihop networks which are characterized by finite link capacities,
finite time horizons, and stochastic traffic.

###  1.1. Related Work

The related work can be divided into three categories: scheduling traffic with
deadlines in single-hop networks, scheduling traffic with deadlines in multi-
hop networks, and works on online stochastic programming.

Single-hop networks. In single-hop networks, traffic between any source-
destination node needs to traverse only one link. If the network has a single
link, there is extensive literature that guarantees a constant approximation
ratio (for example $\frac{e-1}{e}$), e.g., (Chin et al., 2006; Jeż, 2013).
Further, when all packets can be successfully delivered within their deadlines
(called the underloaded regime), simple algorithms such as Earliest-Deadline-
First become optimal (Baruah et al., 1992; Liu and Layland, 1973). In the
absence of interference among links, these techniques can be applied to wired
networks, as each link can be handled independently due to the one-hop
traffic. However, in networks that suffer from interference between links, the
required solutions often become more involved, as the decisions across links
become coupled. The problem has been studied in a variety of works, e.g., (Hou
et al., 2009; Hou, 2013; Kang et al., 2013; Tsanikidis and Ghaderi, 2020,
2021; Singh et al., 2013) under different traffic considerations, benchmarks,
and types of interferences. Many of the prior works consider frame-based
traffic, in which, packets usually are assumed to arrive at the beginning of
the frame and expire at the end of the frame (e.g. (Hou et al., 2009; Hou,
2013)). Recently, advancements on the problem in the more general traffic case
have been made (Kang et al., 2014; Tsanikidis and Ghaderi, 2020),
specifically, it has been shown that obtaining constant fractions of the
optimal objective value is still possible, with approximation ratios ranging
from $0.5$ to $\frac{e-1}{e}$, and complexity that depends on the interference
graph (Tsanikidis and Ghaderi, 2020, 2022)222Although the performance metric
in many of these works is in terms of the real-time capacity region, this
metric can be related to the approximation ratio we study here (Tsanikidis and
Ghaderi, 2021)..

Multi-hop networks. The problem becomes significantly more challenging when
multihop traffic is considered. The past work can be divided into four groups:
(i) heuristics without theoretical guarantees (e.g., (Li and Eryilmaz, 2012;
Liu et al., 2019)), (ii) algorithms for restrictive topologies such as single-
destination tree (Bhattacharya et al., 1997; Mao et al., 2014), (iii)
approaches that provide approximation ratios that diminish as parameters of
the network scale (Li and Eryilmaz, 2012; Li et al., 2009; Liu et al., 2019;
Wang et al., 2011; Mao et al., 2014; Deng et al., 2019; Gu et al., 2021;
Tsanikidis and Ghaderi, 2022), and (iv) techniques with guarantees for
relaxations of the problem (Andrews and Zhang, 1999; Andrews et al., 2000;
Singh and Kumar, 2018) or when certain parameters of the network are scaled to
infinity (Singh and Kumar, 2018). Below, we highlight works with theoretical
guarantees in the latter two groups, which are more relevant to our work.

The recent works in (Deng et al., 2019; Gu et al., 2021; Tsanikidis and
Ghaderi, 2022) provide the best existing guarantees for the problem in the
worst-case traffic setting. The work (Deng et al., 2019) studied the problem
without any packet weights, and provided the first algorithm that can achieve
$\Omega(1/\log L)$ approximation, when the minimum link capacity grows to
infinity ($C_{\min}\rightarrow\infty$). This result matches the lower-bound
provided in (Mao et al., 2014) and hence is asymptotically optimal for the
worst-case traffic. The work (Gu et al., 2021) introduced several variants of
an algorithm called GLS, which in the best case, with
$\mathrm{C}_{\max}=\mathrm{C}_{\min}=\Omega(\log L)$, yield a $\Omega(1/\log
L)$-approximation as well. This line of research was improved further in
(Tsanikidis and Ghaderi, 2022) to more general and improved techniques for
weighted packets. In particular, for $\mathrm{C}_{\min}=\Omega(\log(\rho L))$
an algorithm with $\Omega(1/\log(\rho L))$ approximation ratio was provided,
where $\rho=\frac{w_{\max}}{w_{\min}}$ is the ratio of maximum to minimum
packet weight, $w_{\max}$ and $w_{\min}$.

In addition to the work in the worst-case traffic setting, there is work
(Singh and Kumar, 2018) in the case of stochastic traffic, which considers a
relaxed version of the problem under average link capacity constraints.
However, this relaxation considerably simplifies the problem. In fact, some of
the primary challenges of the studied problem in our work, revolve around
handling the strict capacity constraints as opposed to average constraints.
The work (Singh and Kumar, 2018) shows that the performance loss due to such
relaxation becomes asymptotically small in the regime that the link capacities
and arrival rates are all scaled to infinity. However, the performance in real
networks with finite capacities and finite arrival rates is not clear.

There is also work that relaxes the notion of strict deadlines considered in
our paper, e.g., (Andrews et al., 2000; Andrews and Zhang, 1999). The work in
(Andrews et al., 2000) focuses on underloaded networks, with fixed routes,
where the total packet arrival rate to each link is less than the link’s
capacity (assumed to be 1), and packets do not have strict deadlines. In this
case, all packets can be delivered in a bounded time. The paper proposes
algorithms with guaranteed bounded delay as a function of the arrival rates
and the network size. The work in (Andrews and Zhang, 1999) considers
deadlines but they are allowed to be violated by some factor. The system is
again assumed to be underloaded. It provides algorithms where the violation
factor is bounded by a function of arrival rates and the number of links in
the network. In contrast, in our work, packets exceeding their deadlines have
no utility and therefore are discarded. Further, we do not assume an
underloaded system, e.g., the system might be overloaded, in which case, even
the optimal algorithm has to drop some packets. Finally, we consider different
weights or rewards for packets as opposed to unweighted packets in (Andrews
and Zhang, 1999; Andrews et al., 2000). Note that the underloaded system
assumption considerably simplifies the problem. For example, in the case of
single-hop networks with unit capacity, as mentioned earlier, Earliest-
Deadline-First is optimal when the system is underloaded. However, when
packets have weights and the system is overloaded, the approximation ratio
achieved by any online algorithm is strictly less than one (Jeż, 2013).

Although the above works have made significant advances on the problem, they
raise the following lingering question: Can we do better than an approximation
ratio that becomes increasingly worse as network parameters, such as $L$,
become larger, in the case of finite link capacities and for the given
traffic? In this work, we address this question by providing algorithms that,
in the presence of stochastic traffic, can provide near-optimal performance.

Online Stochastic Programming. Beyond the above related literature, there has
been a line of research (Devanur et al., 2019; Agrawal and Devanur, 2014;
Kesselheim et al., 2014; Banerjee and Freund, 2020) for online and sequential
assignment of limited resources to requests. Requests are drawn from an i.i.d.
distribution of request types, which characterize how these requests can be
assigned to resources. Then, each request of a given request type can be
assigned to the same fixed set of resources regardless of the time of arrival.
While there is similarity between this model and our problem, for example, by
considering “requests” as the arriving packets, and “resources” as each pair
of (link, time slot) combinations, unfortunately, however, these techniques
cannot be applied to our problem. This is because these works assume that
requests of a given request type can be assigned to the same resources,
whereas, in our problem, packets of a given packet type cannot be assigned to
the same (link, time slot) pair. For example, a packet arriving at time
$t^{\prime}$ after $t$ cannot be assigned to the (link, time slot) pair
$(l,t)$. Further, these techniques often require randomizing across all
decision variables, which in our case, is prohibitive, as there are
exponentially many options for scheduling a packet in the network.

Finally, we point out that there is a rich literature on the broader problem
of time-sensitive scheduling and/or routing, e.g., (Leighton et al., 1994;
Wang et al., 2014; Sun et al., 2021; Srinivasan and Teo, 1997; Fountoulakis et
al., 2023; Fan et al., 2019). In (Leighton et al., 1994), the time to schedule
a set of packets (makespan) in a network was studied. Deadline-constrained
scheduling has also been studied jointly with other objectives such as energy
consumption (Wang et al., 2014; Fan et al., 2019) or the age of information
(Sun et al., 2021; Fountoulakis et al., 2023).

###  1.2. Contributions

The main contributions of this paper can be summarized as follows.

Near-Optimal Static Algorithm. We design the first algorithm, to the best of
our knowledge, that guarantees near-optimal performance for scheduling packets
arriving as a stochastic process with arbitrary hard deadlines in multihop
networks with finite and strict link capacities. Assuming the knowledge of the
packet arrival rates, our algorithm relies on solving a single linear program
with polynomial number of variables and constraints. Subsequently, packets are
treated independently from each other and are scheduled over different links
according to forwarding probabilities calculated through the solution of the
linear program. Our algorithm provides $(1-3\epsilon)$-approximation when the
minimum link capacity $\mathrm{C}_{\min}$ satisfies $\mathrm{C}_{\min}\geq
2\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log(L/\epsilon)$, when packet
arrivals are generated by a set of Bernoulli (or Binomial) processes, with any
deadline and weight. Our result does not require relaxing the capacity
constraints or taking a limit on arrival rate or capacity. Furthermore, our
guarantees hold for any finite time horizon of length
$\mathrm{T}\geq\frac{2}{\epsilon}d_{max}^{2}$.

Near-Optimal Dynamic Algorithm. When the knowledge of the packet arrival rates
is not available, we provide a dynamic algorithm by dividing time into
$\lceil\log(1/\epsilon)\rceil$ phases of geometrically increasing lengths. At
the beginning of each phase, the packet arrivals in prior phases are used to
estimate the packet arrival rate of different packet types (flows), which are
subsequently used in the linear program by the static algorithm for the
incoming phase. Our dynamic algorithm preserves the
$(1-O(\epsilon))$-approximation, when the time horizon is larger than a
certain threshold.

Extensions to Non-Stationary and Dependent Traffic. Our techniques are
versatile and can be adapted to other distributional assumptions. For
instance, in the case that packet arrivals exhibit dependence across time
slots or within a time slot, our method achieves near-optimal performance when
$\mathrm{C}_{\min}\geq
2D\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log(L/\epsilon)$, where $D$
represents the degree of dependence among packets (in the i.i.d. scenario,
$D=1$). This allows extending our results to more general distributions on the
number of arrivals per time slot. Additionally, our techniques can be extended
to periodic traffic distributions and non-stationary distributions. The
combination of these extensions results in a wide range of stochastic arrival
processes that can be effectively addressed through our method.

Empirical Evaluation Using Real Datasets. In order to evaluate the performance
of our algorithms in practical settings, in addition to extensive synthetic
simulations, we provide simulations using real traffic traces, and over real
networks. The results indicate that, despite the presence of highly non-
stationary traffic in the network traces, our algorithms demonstrate
significant performance improvement over the worst-case-based algorithms.

###  1.3. Notations

We use $[n]$ to denote the set $\\{1,2,\cdots,n\\}$. Further we denote
$[n]_{0}:=[n]\cup\\{0\\}$. We use $\mathbb{N}:=\\{1,2,3,\cdots\\}$, and
$\mathbb{R}^{+}:=\\{x\in\mathbb{R}|x\geq 0\\}$. We define $a\land
b:=\min\\{a,b\\}$, and $(a)^{+}:=\max(a,0)$.

##  2\. Model and Definitions

Consider a communication network consisting of a set of nodes $\mathcal{V}$
and a set of communication links $\mathcal{L}$ between the nodes, defining a
directed graph $\mathcal{G}=(\mathcal{V},\mathcal{L})$. We assume that time is
slotted, i.e., $t=1,2,3,\cdots$. Each link $\ell=(u,v)\in\mathcal{L}$ has a
capacity $\mathrm{C}_{\ell}$ which is the maximum number of packets per time
slot that can be transmitted over the link from node $u\in\mathcal{V}$ to node
$v\in\mathcal{V}$. To simplify future discussions, we further define
$\overline{\mathcal{G}}:=(\mathcal{V},\overline{\mathcal{L}})$ to denote
$\mathcal{G}$ with the addition of self-loops to the set of edges, i.e.,
$\overline{\mathcal{L}}=\mathcal{L}\cup\\{(u,u):u\in\mathcal{V}\\}$.
Scheduling a packet over a “self-loop link” then has the interpretation that
the packet remains at the same node, and we define $\mathrm{C}_{\ell}=\infty$
for any self-loop link $\ell$. We further use $\mathtt{Out}(v)$ to denote the
set of outgoing links of a node $v$ in $\overline{\mathcal{G}}$, i.e.,
$\mathtt{Out}(v)=\\{(v,u):(v,u)\in\overline{\mathcal{L}}\\}$, and similarly
use $\mathtt{Inc}(v)$ to denote its incoming links, i.e.,
$\mathtt{Inc}(v)=\\{(u,v):(u,v)\in\overline{\mathcal{L}}\\}$. Note that by
definition, self-loop link $(v,v)$ belongs to both $\mathtt{Out}(v)$ and
$\mathtt{Inc}(v)$.

The network is shared by a set of packet types (flows)
$[\mathrm{J}]:=\\{1,2,\cdots,\mathrm{J}\\}$. A packet type (flow)
$j\in[\mathrm{J}]$ is characterized by its source $s_{j}\in\mathcal{V}$, its
destination $z_{j}\in\mathcal{V}$, an end-to-end deadline
$d_{j}\in\mathbb{N}\cup\\{0\\}$, and a weight $w_{j}\in\mathbb{R}^{+}$.
Packets of different types arrive during a time horizon of length
$\mathrm{T}$. A packet of type $j$ arriving at the beginning of time slot $t$
has to reach its destination before the end of time slot $d_{j}+t$ to yield a
reward $w_{j}$, otherwise it is discarded. Packets of type $j\in[\mathrm{J}]$
arrive according to a stochastic process $\\{a_{j}^{t}\\}_{t\in[\mathrm{T}]}$
with possibly time-dependent arrival rate $\lambda_{j}^{t}$, where
$a_{j}^{t}\geq 0$ is the number of packet arrivals of type $j$ at time $t$,
and $\lambda_{j}^{t}=\mathbb{E}[{a_{j}^{t}}]$. We denote the maximum deadline,
the maximum number of arrivals, and the minimum arrival rate of any packet
type by $d_{max}:=\max_{j\in[\mathrm{J}]}d_{j}$,
$a_{max}:=\max_{t\in[\mathrm{T}],j\in[\mathrm{J}]}a_{j}^{t}$ and
$\lambda_{\min}:=\min_{t\in[\mathrm{T}],j\in[\mathrm{J}]}\lambda_{j}^{t}$,
respectively.

A packet of type $j$ arriving at time $t$, which is scheduled over the
network, should be transmitted over a sequence of links, at specified time
slots, that deliver it from its source $s_{j}$ to its destination $z_{j}$
before the end of time slot $t+d_{j}$. This sequence of links and time slots
is characterized through the notion of relative route-schedule, defined in
Definition 1 below.

###### Definition 0 (Relative Route-Schedule).

A (relative) route-schedule $k$ for a packet of type $j$ is a directed walk
with sequence of edges $[k_{0}\quad k_{1}\quad\cdots\quad k_{d_{j}}]$ in graph
$\overline{\mathcal{G}}$, where $k_{\tau}\in\overline{\mathcal{L}}$ is the
link over which the packet is scheduled at the $\tau$-th time slot following
its arrival, with $k_{0}\in\mathtt{Out}(s_{j})$ and
$k_{d_{j}}\in\mathtt{Inc}(z_{j})$. We use
$\mathcal{K}_{j}\subseteq\overline{\mathcal{L}}^{d_{j}+1}$ to denote the set
of all valid route-schedules for packet-type $j\in[\mathrm{J}]$.

For notational convenience, we define $k_{\tau}=\varnothing$ for $\tau>d_{j}$,
where $k\in\mathcal{K}_{j}$. Our objective is to maximize the weighted sum of
the packets that are successfully delivered from their sources to their
destinations within their deadlines. Given a random instance of the packet
arrival sequence, this optimization can be formulated as an integer program
over the time horizon $\mathrm{T}$, defined in (1a)-(1d) below. We refer to
this optimization as $\mathrm{RI}(\mathrm{T})$ which stands for Random Integer
problem over the time horizon of length $\mathrm{T}$. In
$\mathrm{RI}(\mathrm{T})$, the optimization is over the scheduling decisions
$\mathbf{y}=\\{y_{jk}^{nt}\\}$, where each $y_{jk}^{nt}$ is a binary variable
indicating whether the $n$-th arriving packet of type $j$ at time $t$ is
scheduled using a relative route-schedule $k\in\mathcal{K}_{j}$ or not.

(1a) $\displaystyle\max_{\mathbf{y}}\quad$
$\displaystyle\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\sum_{n=1}^{a_{j}^{t}}\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt}\quad\quad\quad(:=\mathrm{RI}(\mathrm{T}))$
(1b) s.t. $\displaystyle\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt}\leq
1,\quad\forall t\in[\mathrm{T}],\forall j\in[\mathrm{J}],\forall
n\in[a_{j}^{t}],$ (1c)
$\displaystyle\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}\sum_{j=1}^{\mathrm{J}}\sum_{n=1}^{a_{j}^{\tau}}\sum_{k:k_{t-\tau}=\ell}y_{jk}^{n\tau}\leq\mathrm{C}_{\ell},\quad\forall\ell\in\mathcal{L},\forall
t\in[\overline{\mathrm{T}}],$ (1d) $\displaystyle
y_{jk}^{nt}\in\\{0,1\\},\quad\forall j\in[\mathrm{J}],\forall
t\in[\mathrm{T}],\forall k\in\mathcal{K}_{j}.$

The objective function (1a) is the sum of the weights of arriving packets
which are successfully scheduled in the network. Constraints (1b) and (1d)
require that each arriving packet is assigned to “at most” one route-schedule
(or not scheduled at all). Constraints (1c) require each link’s capacity to be
enforced at every time slot $t\in[\overline{\mathrm{T}}]$ that packets can
exist in the system, where $\overline{\mathrm{T}}:=\mathrm{T}+d_{max}$.
Specifically, in the left-hand-side of (1c), we count the number of packets
that are scheduled to be transmitted on link $\ell$ at time $t$. By the
definition of deadline, these packets could have only arrived in the last
$d_{max}+1$ time slots, i.e., at times $\tau\in\\{t,t-1,\cdots,t-d_{max}\\}$,
and are counted in the left-hand-side of (1c) if they are scheduled using some
route-schedule $k$ which, $t-\tau$ slots after the packet arrival at time
$\tau$, transmits it on $\ell$ at time $t$ (i.e. the route schedules with
$k_{t-\tau}=\ell$). Refer to Figure 1 for an illustration of the notations
through an example.

(a) Example network graph $\mathcal{G}$

(b) Graph $\overline{\mathcal{G}}$ with self-loops

Figure 1. A network $\mathcal{G}$ with nodes $\\{u,a,b,v\\}$ and $4$ links
with unit capacity, and the corresponding graph $\overline{\mathcal{G}}$ with
self-loops. Suppose there are 2 packet types $[\mathrm{J}]=\\{1,2\\}$ with
sources $s_{1}=u,s_{2}=a$, destinations $z_{1}=z_{2}=v$, and deadlines
$d_{1}=2$ and $d_{2}=1$. Two possible route-schedules for packet-types $1$ and
$2$ respectively are: $k=[(u,u)\quad(u,a)\quad({a},{v})]$ and
$k^{\prime}=[({a},v)\quad(v,v)]$. Suppose there are arrivals of type $1$ and
type $2$ at time slots $t_{0}$ and $t_{0}+2$ respectively (i.e.
$a_{1}^{t_{0}}=a_{2}^{t_{0}+2}=1$), then, scheduling both packets on
$k,k^{\prime}$ (i.e. $y_{1k}^{1t}=y_{2k^{\prime}}^{1(t+2)}=1$) is not possible
due to insufficient capacity on link $(a,v)$ at time $t+2$ (violating
constraint (1c)).

Our goal is to provide algorithms that guarantee good performance in the
average sense, as formalized through the definition below.

###### Definition 0.

Suppose the optimal objective value of $\mathrm{RI}(\mathrm{T})$ is
$W_{\mathrm{RI}(\mathrm{T})}^{\star}$. An algorithm $\mathrm{ALG}$ provides
$\gamma$-approximation to $\mathrm{RI}(\mathrm{T})$, if the objective value
attained using $\mathrm{ALG}$, $W_{\mathrm{ALG}}$, satisfies:

$\mathbb{E}\big{[}W_{\mathrm{ALG}}\big{]}\geq\gamma\mathbb{E}\big{[}W_{\mathrm{RI}(\mathrm{T})}^{\star}\big{]},$

where the expectation in
$\mathbb{E}\big{[}W_{\mathrm{RI}(\mathrm{T})}^{\star}\big{]}$is with respect
to the randomness in the random instance $\mathrm{RI}(\mathrm{T})$ (i.e., with
respect to the arrival sequence $\\{a_{j}^{t}\\}$), and
$\mathbb{E}\left[W_{\mathrm{ALG}}\right]$ is with respect to the randomness in
the random instance, and, if applicable, the random decisions of
$\mathrm{ALG}$.

##  3\. Algorithms and Main Results

In this section, we introduce the two main algorithms under the assumption of
fixed packet arrival rates, i.e., $\lambda_{j}^{t}\equiv\lambda_{j}$.
Extensions to this assumption are discussed in Section 5.1. We first introduce
Algorithm 1 (Flow-Based Probabilistic Forwarding) for the case of known packet
arrival rates, and subsequently Algorithm 2 (Dynamic Learning with
Probabilistic Forwarding) for the case where the packet arrival rates are
unknown. For each algorithm, we state its corresponding approximation ratio.

###  3.1. Algorithm 1: Flow-Based Probabilistic Forwarding

We first introduce Algorithm 1, a scheduling algorithm that probabilistically
forwards packets at each time slot based on their age (time since their
arrival), type, and their current node in the system. Specifically, Algorithm
1 forwards a type-$j$ packet, which is currently at node $v$, over an outgoing
link $\ell\in\mathtt{Out}(v)$ with a probability determined by forwarding
variables $f_{j\ell}^{\tau}$, which depend on the packet type $j$ and the
number of time slots $\tau\in\\{0,1,\cdots,d_{j}\\}$ that the packet has been
in the system so far (i.e., its age). To determine the forwarding
probabilities, Algorithm 1 solves a Linear Program (LP), as a preprocessing
step. In the following, we first explain the LP and subsequently discuss the
forwarding process based on its solution.

The forwarding variables $\mathbf{f}=\\{f_{j\ell}^{\tau}\\}$ are computed as
the solution to the following LP, which we refer to as $\mathrm{F}_{S}$:

(2a) $\displaystyle\max_{\mathbf{f}}\quad$
$\displaystyle\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0}\quad(:=\mathrm{F}_{S})$
(2b) s.t. $\displaystyle\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0}\leq
1,\quad\forall j\in[\mathrm{J}],$ (2c)
$\displaystyle\sum_{\ell\in\mathtt{Inc}(v)}f_{j\ell}^{\tau-1}=\sum_{\ell\in\mathtt{Out}(v)}f_{j\ell}^{\tau},\quad\forall
v\in\mathcal{V},\tau\in[d_{j}],\forall j\in[\mathrm{J}],$ (2d) $\displaystyle
f_{j\ell}^{0}=0,\ \ \forall\ell\not\in\mathtt{Out}(s_{j}),\quad
f_{j\ell}^{d_{j}}=0,\ \ \forall\ell\not\in\mathtt{Inc}(z_{j}),$ (2e)
$\displaystyle\sum_{j=1}^{\mathrm{J}}\sum_{\tau=0}^{d_{j}}\lambda_{j}f_{j\ell}^{\tau}\leq\frac{\mathrm{C}_{\ell}}{1+\epsilon},\quad\forall\ell\in\mathcal{L},$
(2f) $\displaystyle f_{j\ell}^{\tau}\geq 0,\quad\forall
j\in[\mathrm{J}],\forall\ell\in\overline{\mathcal{L}},\forall\tau\in[d_{max}]_{0}.\quad$

The objective of $\mathrm{F}_{S}$ (2a) is to maximize the expected weighted
number of arriving packets (i.e. packets with age $\tau=0$) that are
forwarded, which corresponds to the packets attempted to be scheduled.
Constraint (2b) states that the total forwarding probability of arriving
packets of type $j$ out of its source node is at most $1$. If the constraint
is not tight, there is a chance that the packet will not be forwarded, and
therefore be dropped. The constraints (2c) require that packets that have been
in the system for $\tau-1$ time slots, which have been forwarded into some
node $v$, must be forwarded out of the node $v$, or remain in $v$ (through a
self-loop), in the following slot, with their age increasing to $\tau$. This
ensures that packets are continuously forwarded, hop by hop, until their
expiration333Due to the definition of route-schedules $k$ using the self-
looped graph $\overline{\mathcal{G}}$, packets of type $j$ that are delivered
earlier than their deadline to their destination $z_{j}$, are “scheduled” over
$(z_{j},z_{j})$ until their expiry.. Constraints (2d) ensure that all arriving
packets (i.e., with age $\tau=0)$ can only be forwarded out of their sources
(first constraint), and must only be delivered into their destinations by
their deadline (second constraint). Finally, (2e) enforces the capacity
constraint for each link. It requires that a link’s average total traffic, due
to the forwarding of packets of different ages and packet types, cannot be
more than the link’s capacity, divided by $(1+\epsilon)$ for some
$\epsilon>0$. The notations and constraints of $\mathrm{F}_{S}$ are
illustrated through an example in Figure 2. Note that $\mathrm{F}_{S}$ is
independent of the time horizon $T$, unlike $\mathrm{RI}(\mathrm{T})$.

(a) Forwarding variables $f_{j\ell}^{\tau}$, $\tau\in\\{0,1,2=d_{j}\\}$

(b) Links’ capacity consumption

Figure 2. For the network in Figure 1, consider a packet type $j$ with arrival
rate $\lambda_{j}=0.5$, source $s_{j}=u$, destination $z_{j}=v$ and deadline
$d_{j}=2$. Figure 2(a) shows a feasible forwarding variables assignment
$\mathbf{f}$ for $\mathrm{F}_{S}$ (2a)-(2f) (with, e.g., $f_{j(u,a)}^{0}=0.7$
and $f_{j(v,v)}^{2}=0.3$). The capacity constraints (2e) should sum the total
flow on the link over different ages $\tau$, which is illustrated in Figure
2(b). For example, on link $(a,v)$, the total flow is
$\big{(}f_{j(a,v)}^{1}+f_{j(a,v)}^{2}\big{)}\overline{a}_{j}=(0.3+0.2)0.5$.

0.1 Input: Packet types $\\{(s_{j},z_{j},d_{j},w_{j},\lambda_{j})\\}$, and an
$\epsilon>0$.

0.2Find optimal solution $\mathbf{f}^{\star}=\\{f_{j\ell}^{t\star}\\}$ to the
LP $F_{S}$ (2a)-(2f).

0.3for  _each time $t=1,2,\cdots,\mathrm{T}$_ do

0.4 for  _each unexpired packet of type $j$ at each node $v\in\mathcal{V}$_ do

0.5 if  _age of packet is $\tau=0$_ then

0.6 Select link $\ell\in\mathtt{Out}(v)$ w.p. $f_{j\ell}^{0\star}$. Else, drop
w.p. $\big{(}1-\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0\star}\big{)}$.

0.7 end if

0.8 else if  _age of packet is $1\leq\tau\leq d_{j}$_ then

0.9 Select a link $\ell\in\mathtt{Out}(v)$ w.p.
$\frac{f_{j\ell}^{\tau\star}}{\sum_{\ell\in\mathtt{Out}(v)}f_{j\ell}^{\tau\star}}$

0.10 end if

0.11 Forward packet over selected link $\ell$ if there is available capacity
on $\ell$ at $t$, else drop.

0.12 end for

0.13

0.14 end for

Algorithm 1 Flow-Based Probabilistic Forwarding ($\mathrm{FBPF})$

After solving the LP $F_{S}$, we obtain a solution
$\mathbf{f}^{\star}=\\{f_{j\ell}^{\tau\star}\\}$, which we use in Algorithm 1
to schedule unexpired packets at every time slot $t$ and at every node $v$
(Algorithm 1-Algorithm 1 in Algorithm 1). Specifically, for new packets (i.e.,
with age $\tau=0$, Algorithm 1), we select their first link
$\ell\in\mathtt{Out}(s_{j})$ with probability $f_{j\ell}^{0\star}$, or drop
the packet with probability
$1-\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0\star}$ (Algorithm 1). For
packets with age $\tau>0$, we always select one of the outgoing links
probabilistically (Algorithm 1). We then forward the packet over the selected
link if there is available capacity (Algorithm 1). Overall each packet is
forwarded hop by hop using variables $\\{f_{j\ell}^{\tau\star}\\}$, until it
reaches its destination. Although the algorithm only rejects packets when they
arrive, once a packet is admitted, in subsequent time slots, it might also be
dropped due to the capacity constraints (Algorithm 1).

Theorem 1 states the performance guarantee of Algorithm 1 in the case that the
arrivals $a_{j}^{t}$ for each type $j$ are i.i.d. Bernoulli or Binomial.

###### Theorem 1.

Given an $\epsilon\in(0,1/3)$, $\mathrm{FBPF}$ provides
$(1-3\epsilon)$-approximation to $\mathrm{RI}(\mathrm{T})$ when
$\mathrm{C}_{\min}\geq
2\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log(L/\epsilon)$ and
$\mathrm{T}\geq\frac{2d_{max}^{2}}{\epsilon}$.

###### Remark 1.

Based on Theorem 1, we require $\mathrm{C}_{\min}\geq\Omega(\frac{\log
L/\epsilon}{\epsilon^{2}})$ for $(1-\epsilon)$-approximation to
$\mathrm{RI}(\mathrm{T})$ on average. This is in vast contrast to prior
literature (Tsanikidis and Ghaderi, 2022; Deng et al., 2019) that for
$\mathrm{C}_{\min}=\Omega(\log({\rho L}))$, resulted in $\Omega(1/\log({\rho
L}))$ approximation ratios for the worst case, where $\rho=w_{max}/w_{\min}$
was the ratio of maximum weight to minimum weight of packets. These findings
indicate that, in situations where the traffic is stochastic, the results
presented in prior work can be overly pessimistic and significant improvements
can be made. Further, compared to (Singh and Kumar, 2018) with stochastic
traffic, our results are not asymptotic, i.e., we do not require scaling the
arrival rates, capacities, and time horizon to infinity in order to obtain the
stated performance guarantees.

###### Remark 2.

Algorithm 1 admits a distributed implementation following the computation of
$\mathrm{F}_{S}$. Every intermediate node only requires storing the forwarding
variables for each packet type based on their age.

###### Remark 3.

Our results extend to unreliable networks, where the capacity of link $\ell$
at each time is random with average $\mathrm{C}_{\ell}p_{\ell}$ for some
$p_{\ell}\in(0,1]$. Algorithm 1 can be directly applied by replacing
$\mathrm{C}_{\ell}$ with $\mathrm{C}_{\ell}p_{\ell}$ in $F_{S}$ (2a)-(2f). In
the case that the link’s capacity is i.i.d.
$Binomial(\mathrm{C}_{\ell},p_{\ell})$, we obtain a
$(1-\epsilon)$-approximation ratio when
$\min_{\ell}\\{\mathrm{C}_{\ell}p_{\ell}\\}=\Omega(\frac{1}{\epsilon^{2}}\log\frac{L}{\epsilon})$.
Conceptually, the effective number of packets that can be scheduled on link
$\ell$ now becomes $\mathrm{C}_{\ell}p_{\ell}$, and therefore, we require this
product to satisfy a similar bound as in Theorem 1.

###### Remark 4.

Theorem 1 can be generalized to non-stationary arrival rates, dependent packet
arrivals, and more general distributions than Bernoulli and Binomial. This is
described in Section 5.1. In particular, as a result of these extensions, in
cases where a very short horizon is of interest, a non-stationary solution
will still guarantee a near-optimal performance, in contrast to Algorithm 1.

###  3.2. Algorithm 2: Dynamic Learning with Probabilistic Forwarding

In this section, we provide a dynamic variant of Algorithm 1 which does not
require the knowledge of the arrival rates
$\\{\lambda_{j}:j\in[\mathrm{J}]\\}$. Instead, it estimates the arrival rates
dynamically over time and uses the estimates as input to Algorithm 1. The
dynamic variant is described in Algorithm 2.

Algorithm 2 works by dividing the time horizon $\mathrm{T}$ into $(\Phi+1)$
phases, $\phi=0,1,\cdots,\Phi$, where $\Phi=[\log_{2}(1/\epsilon)]$ for some
suboptimality parameter $\epsilon\in(0,1)$. The phases have geometrically
increasing duration, with phase $\phi$ for $\phi\geq 1$ having a duration
$d_{max}+\mathrm{T}_{\phi}$, where
$\mathrm{T}_{\phi}=\mathrm{T}_{0}2^{\phi-1}$, and phase $\phi=0$ having
duration $\mathrm{T}_{0}=\epsilon(\mathrm{T}-d_{max}\Phi)$. An illustration of
the time division into phases is provided in Figure 3. During the first phase
$\phi=0$, due to the lack of any statistics, the algorithm can either remain
idle, or schedule packets in an arbitrary way (for example, according to
algorithms in (Tsanikidis and Ghaderi, 2022; Gu et al., 2021)). At the end of
the first phase, the algorithm estimates $\hat{\lambda}_{j}$ (see Algorithm 2
of Algorithm 2) for the first time, using the number of packet arrivals of
different packet types during the phase. The optimal forwarding probabilities
$\mathbf{\hat{f}}=\\{\hat{f}_{j\ell}^{\tau\star}\\}$ are then found by solving
$\mathrm{F}_{S}$ (as defined in (2a)-(2f)) using the current estimates
$\hat{\lambda}_{j}$ instead of the actual arrival rates $\lambda_{j}$ (after
invoking Algorithm 2). The packets during the following phase are then
scheduled according to the forwarding probabilities
$\mathbf{\hat{f}}=\\{\hat{f}_{j\ell}^{\tau\star}\\}$ as in Lines 1-1 of
Algorithm 1 (Algorithm 2). More specifically, in Algorithm 2, we calculate a
shrinkage parameter which is used to further scale down the right-hand-side of
(2e) (Algorithm 2), based on the confidence of our estimates. We allow a
tuning parameter $\zeta\geq 0$ to reduce the effect of the shrinkage, with
$\zeta=0$ yielding no additional shrinkage, as in Algorithm 1 to obtain
Theorem 1. Similarly, in each new phase, we repeat the above by updating
estimates $\hat{\lambda}_{j}$ and therefore the forwarding probabilities.

Figure 3. The division of time horizon into phases of geometrically increasing
duration in Algorithm 2.

0.1

0.2Input: Suboptimality parameter $\epsilon\in(0,1)$. Tuning parameter
$\zeta\geq 0$.

0.3$\Phi\leftarrow\left[\log_{2}(1/\epsilon)\right]$ $\blacktriangleright$
number of phases

0.4$\mathrm{T}_{0}\leftarrow\epsilon(\mathrm{T}-d_{max}\Phi)$

0.5Run any scheduling algorithm for $\mathrm{T}_{0}$ slots while collecting
statistics.

0.6for  _each time phase $\phi=1,2,\cdots,\Phi$_ do

0.7

0.8 Update estimates $\hat{\lambda}_{j}$ of $\lambda_{j}$ using historical
samples of $a_{j}^{t}$ for all $j\in[\mathrm{J}]$

0.9 $\mathrm{T}_{\phi}\leftarrow\mathrm{T}_{0}2^{\phi-1}$

0.10
$\eta\leftarrow(1-\zeta\sqrt{\log(2\mathrm{J}/\epsilon)/(\lambda_{\min}\mathrm{T}_{\phi})})$

0.11 Run $\mathrm{FBPF}$, with estimates $\hat{\lambda}_{j}$ and replacing
$\mathrm{C}_{\ell}$ with $\eta\mathrm{C}_{\ell}$ in (2e), for
$(d_{max}+\mathrm{T}_{\phi})$ time slots.

0.12 end for

Algorithm 2 Dynamic Learning with Probabilistic Forwarding ($\mathrm{DLPF}$)

We remark that in Algorithm 2, estimating the arrival rates can be performed
based on a simple empirical average of the number of packet arrivals of
different types using the past observations. The performance of Algorithm 2 is
described in Theorem 2 in the case that the arrivals $\\{a_{j}^{t}\\}$ are
i.i.d. Bernoulli or Binomial random variables. For sufficiently large horizon
$\mathrm{T}$, the algorithm is near-optimal, and only requires solving $F_{S}$
a total of $[\log_{2}(1/\epsilon)]$ times.

###### Theorem 2.

Algorithm 2 for parameters $\epsilon>0$ and $\zeta=1$, yields a
$(1-10\epsilon)$-approximation to $\mathrm{RI}(\mathrm{T})$ when
$\mathrm{C}_{\min}\geq
2\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log{\frac{L}{\epsilon}}$ and

$\mathrm{T}=\Omega\left(\frac{a_{\max}^{2}}{2\lambda_{\min}^{2}\epsilon^{2}}\log(\frac{2\mathrm{J}}{\epsilon})+2d_{max}^{2}\log(1/\epsilon)/\epsilon\right).$

###### Remark 5.

Algorithm 2 provides a baseline approach under stationary model assumptions.
For non-stationary packet arrival rates, a more natural approach is to update
the estimates of $\hat{\lambda}_{j}$ periodically, for example, by using a
fixed $\mathrm{T}_{\phi}$, or combining past estimates $\hat{\lambda}_{j}$
such that more recent samples are prioritized. Further the estimates in
Algorithm 2 could leverage prior distributional assumptions, e.g., through a
Bayesian approach, or use more complex machine learning methods.

###### Remark 6.

An alternative simpler two-phase algorithm, which splits the horizon into two
phases $\mathrm{T}_{0}$ and $\mathrm{T}_{1}$ with appropriate lengths, can
achieve $1-O(\epsilon)$ approximation ratio, but, for small $\epsilon$, the
required time horizon will be roughly a factor $1/\epsilon$ lager than the
stated $\mathrm{T}$ in Theorem 2.

###### Remark 7.

Although a small $\lambda_{\min}$ suggests the need of a large horizon in
Theorem 2, this is a consequence of a multiplicative $(1-O(\epsilon))$
approximation ratio considered here for uniformity across our results. In
practice this typically is not required as validated in our simulations in
Section 6. In particular, packets with very small $\lambda_{\min}$ effectively
have a small impact on the performance of the algorithm, and can be handled as
special cases if they have big weights, or are otherwise entirely ignored.
This can be formalized by modifying the approximation ratio definition, such
that it is satisfied up to additive constants. Finally, we remark that
although in our analysis we assume the knowledge of $\lambda_{\min}$ for
simplicity, we can extend the analysis and remove this assumption by replacing
$\lambda_{\min}$ with its estimate $\hat{\lambda}_{\min}$. In practice,
$\zeta$ can either be set to $0$ or be set such that the knowledge of
$\lambda_{\min}$ is not required.

##  4\. Proofs of Main Theorems

In this section, we present the proof of Theorem 1 and the outline for the
proof of Theorem 2.

###  4.1. Proof of Theorem 1 for Algorithm 1

The proof of Theorem 1 has 4 steps. In Step 1, we construct an intermediate
LP, referred to as $\mathrm{EI}(\mathrm{T})$, whose optimal objective value is
higher than that of $\mathrm{RI}(\mathrm{T})$. In Step 2, we argue that this
LP admits a simple near-optimal stationary solution, thus simplifying
$\mathrm{EI}(\mathrm{T})$ significantly into a new LP called $\mathrm{E}_{S}$.
In Step 3, we convert the LP $\mathrm{E}_{S}$ into a flow-based LP
$\mathrm{F}_{S}$, which is the LP defined in (2a)-(2f) for Algorithm 1.
Finally, in Step 4, we argue that by using the forwarding variables given by
$\mathrm{F}_{S}$ appropriately, we obtain Theorem 1. We describe each step
below in detail.

Step 1.  We first construct a deterministic LP below, which we refer to as the
Expected Instance, $\mathrm{EI}(\mathrm{T})$:

(3a) $\displaystyle\max_{\mathbf{x}}\quad$
$\displaystyle\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}^{t}\sum_{k\in\mathcal{K}_{j}}x_{jk}^{t}\quad\quad\quad(:=\mathrm{EI}(\mathrm{T}))$
(3b) s.t. $\displaystyle\sum_{k\in\mathcal{K}_{j}}x_{jk}^{t}\leq
1,\quad\forall t\in[\mathrm{T}],\forall j\in[\mathrm{J}],$ (3c)
$\displaystyle\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}\sum_{j=1}^{\mathrm{J}}\lambda_{j}^{\tau}\sum_{k\in
K_{j}:k_{t-\tau}=\ell}x_{jk}^{\tau}\leq\mathrm{C}_{\ell},\quad\forall\ell\in\mathcal{L},\forall
t\in[\overline{\mathrm{T}}],$ (3d) $\displaystyle x_{jk}^{t}\geq
0,\quad\forall j\in[\mathrm{J}],\forall t\in[\mathrm{T}],\forall
k\in\mathcal{K}_{j}.$

In $\mathrm{EI}(\mathrm{T})$, the variables $\mathbf{x}=\\{x_{jk}^{t}\\}$ can
be interpreted as probabilistic (relaxed) versions of $\\{y_{jk}^{nt}\\}$ (in
$\mathrm{RI}(\mathrm{T})$), i.e., $x_{jk}^{t}$ can be viewed as the
probability of setting $y_{jk}^{nt}=1$ for each $n=1,\ldots,a_{j}^{t}$. Then,
$\mathrm{EI}(\mathrm{T})$ can be viewed as an LP for maximizing the expected
reward if every packet type $j$ has
$\mathbb{E}\big{[}a_{j}^{t}\big{]}=\lambda_{j}^{t}$ arrivals at each time slot
$t$, and the capacity constraints are only satisfied in expectation.

Lemma 1 establishes a relationship between $\mathrm{RI}(\mathrm{T})$ and
$\mathrm{EI}(\mathrm{T})$, for general distributions on $\\{a_{j}^{t}\\}$.

###### Lemma 0.

Let $W_{\mathrm{EI}(\mathrm{T})}$ denote the optimal value of the expected
instance $\mathrm{EI}(\mathrm{T})$, and $W_{\mathrm{RI}(\mathrm{T})}$ denote
the optimal value of $\mathrm{RI}(\mathrm{T})$. Then, for any distribution on
$\\{a_{j}^{t},t\in[\mathrm{T}],j\in[\mathrm{J}]\\}$, with
$\lambda_{j}^{t}<\infty$, we have

$W_{\mathrm{EI}(\mathrm{T})}\geq\mathbb{E}\left[W_{\mathrm{RI}(\mathrm{T})}\right].$

###### Proof of Lemma 1.

The proof is standard, and provided in Section A.1 for completeness. ∎

###### Remark 8.

$W_{\mathrm{EI}(\mathrm{T})}$ only depends on the arrival rates
$\lambda_{j}^{t}$ and not on the exact distribution of $\\{a_{j}^{t}\\}$, and
therefore, it serves as a universal upper bound on
$\mathbb{E}\left[W_{\mathrm{RI}(\mathrm{T})}\right]$ for all distributions
with identical arrival rates.

We note that the expected-instance relaxations have been adopted successfully
in several other problems to obtain a bound on an integer program (Jiang et
al., 2020; Sun et al., 2020; Devanur et al., 2019; Gallego and Van Ryzin,
1994). Here, we aim to construct a scheduling algorithm that obtains a total
reward close to the optimal value of $\mathrm{EI}(\mathrm{T})$, which implies
an approximation to $\mathrm{RI}(\mathrm{T})$ due to Lemma 1. However, we
cannot directly solve $\mathrm{EI}(\mathrm{T})$, as this LP has a number of
constraints which grow with the time horizon $\mathrm{T}$, and moreover, the
number of variables could be exponential, due to the exponentially large
number of route schedules. We address these issues in Step 2 and Step 3 below.

Step 2. To resolve the issue of the number of capacity constraints (3c) in
$\mathrm{EI}(\mathrm{T})$ scaling with $\mathrm{T}$, we show that in the
stationary case of $\lambda_{j}^{t}\equiv\lambda_{j}$, a time-independent
(stationary) solution is near-optimal for $\mathrm{EI}(\mathrm{T})$, thus
allowing us to simplify $\mathrm{EI}(\mathrm{T})$. We state the result in
Lemma 2 below.

###### Lemma 0.

In the stationary case, i.e., $\lambda_{j}^{t}\equiv\lambda_{j}$, there is a
stationary solution $x_{ik}^{t}=x_{ik}^{\star}$ for $\mathrm{EI}(\mathrm{T})$
with objective value $\overline{W}_{\mathrm{EI}(\mathrm{T})}$ that satisfies

(4) $W_{\mathrm{EI}(\mathrm{T})}\geq\overline{W}_{\mathrm{EI}(\mathrm{T})}\geq
W_{\mathrm{EI}(\mathrm{T})}\left(1-2\frac{d_{max}^{2}}{\mathrm{T}}\right).$

###### Proof of Lemma 2.

The first inequality in (4) follows directly since limiting our solutions to
static solutions is an additional constraint on the maximization, and
therefore it can only result in a lower objective value.

We prove the second inequality in (4) in three steps. We first provide an
overview of the steps. In Step A, we derive from $\mathrm{EI}(\mathrm{T})$, a
new optimization problem $\mathrm{\overline{EI}}(\mathrm{T})$, which is
identical in most of the capacity constraints, except for the capacity
constraints for times $t<d_{max}$ or $t>\mathrm{T}$ which are modified such
that they contain more terms (therefore $\mathrm{\overline{EI}}(\mathrm{T})$
is more constrained, and any solution to it, will also be feasible for
$\mathrm{EI}(\mathrm{T})$). More specifically, these constraints are modified
appropriately such that $\mathrm{\overline{EI}}(\mathrm{T})$ is symmetrical
with respect to time shifts of the solution. In Step B, we leverage the
symmetry of $\mathrm{\overline{EI}}(\mathrm{T})$ to argue that it admits a
stationary optimal solution. In Step C, we argue that the optimal value for
$\mathrm{\overline{EI}}(\mathrm{T})$ is close to that of
$\mathrm{EI}(\mathrm{T})$. Since any feasible solution to
$\mathrm{\overline{EI}}(\mathrm{T})$ is also feasible for
$\mathrm{EI}(\mathrm{T})$, and $\mathrm{\overline{EI}}(\mathrm{T})$ admits a
stationary optimal solution with objective value close to the optimal of
$\mathrm{EI}(\mathrm{T})$, then there is a near-optimal static solution to
$\mathrm{EI}(\mathrm{T})$ (as in the statement of Lemma 2).

Step A. Note that $\mathrm{EI}(\mathrm{T})$ is nearly symmetrical w.r.t. time
shifts. With a modification to this optimization program, we arrive to a
symmetrical optimization program $\mathrm{\overline{EI}}(\mathrm{T})$.
Consider the capacity constraints (3c) in $\mathrm{EI}(\mathrm{T})$. If
$t<d_{max}$ or $t>\mathrm{T}$, due to the lower or upper limit of the outer
summation being clipped to $0$ or $\mathrm{T}$, the symmetry of the problem is
broken. We can make the program cyclically symmetrical, by appropriately
extending the outer summation when it is clipped, and using the modulo
operator. Specifically, consider the modified capacity constraints for
$\mathrm{\overline{EI}}(\mathrm{T})$:

(5)
$\sum_{\tau=(t-d_{max})}^{t}\sum_{j=1}^{\mathrm{J}}\lambda_{j}\sum_{k:k_{t-\tau}=\ell}x_{jk}^{\tau\
\mathsf{mod}\ (\mathrm{T}+1)}\leq\mathrm{C}_{\ell},\forall\ell,\forall
t\in\mathrm{T}.$

Note that this modification only impacts the constraints (3c) for $t<d_{max}$
or $t>\mathrm{T}$ in $\mathrm{EI}(\mathrm{T})$. These constraints are replaced
with the stronger constraints (5) for $t\in[d_{max}]$ in
$\mathrm{\overline{EI}}(\mathrm{T})$. This modification is visualized in
Figure 4.

Step B. Here, we prove the existence of a stationary optimal solution in
$\mathrm{\overline{EI}}(\mathrm{T})$. To see this, consider any optimal
solution to $\mathrm{\overline{EI}}(\mathrm{T})$ (which might be time-
dependent), denoted with $\mathbf{x}^{\star}=\\{x_{jk}^{t\star}\\}$. Let us
denote a time-shift of all variables of $\mathbf{x}^{\star}$ with $t^{\prime}$
slots by $\mathbf{x}[t^{\prime}]$, i.e.,
$x_{lk}^{t}[t^{\prime}]:=x_{jk}^{((t+t^{\prime})\textsf{mod}\mathrm{T})\star}$.
Note that $\mathbf{x}[0]=\mathbf{x}^{\star}$. As the program is time-shift
invariant, $\mathbf{x}[t^{\prime}]$ should also be optimal for all
$t^{\prime}$. Note that there are at most $\mathrm{T}$ distinct shifts for
$\mathbf{x}[\cdot]$ which are optimal solutions. Since
$\mathrm{\overline{EI}}(\mathrm{T})$ is an LP, the average of a set of optimal
solutions yields a feasible solution that is also optimal. Hence,
$\overline{\mathbf{x}}^{\star}:=\frac{1}{\mathrm{T}}\sum_{t^{\prime}=1}^{\mathrm{T}}\mathbf{x}[t^{\prime}]$
is an optimal solution. This solution is stationary, as it can easily be
verified that it corresponds to the time-average of the original variables,
i.e.,
$\overline{x}_{jk}^{\star}=\frac{\sum_{t=1}^{\mathrm{T}}x_{jk}^{t}}{\mathrm{T}}.$

Step C. Finally, we finish the proof of Lemma 2. Note that any feasible
solution of $\mathrm{EI}(\mathrm{T})$ can be converted to a feasible solution
for $\mathrm{\overline{EI}}(\mathrm{T})$ by setting all variables $x_{jk}^{t}$
for times $t<d_{max}$ or $t>\mathrm{T}-d_{max}$ to $0$. Indeed, this results
in a feasible solution for $\mathrm{\overline{EI}}(\mathrm{T})$ since the only
constraints that are different between $\mathrm{\overline{EI}}(\mathrm{T})$
and $\mathrm{EI}(\mathrm{T})$, only involve the variables at these boundary
time-slots. Therefore we need to set the variables to $0$ for a range of
$2d_{max}$ time slots. The maximum loss in reward, however, due to setting
these variables to $0$, is, per time slot, at most
$W_{\mathrm{EI}(\mathrm{T})}d_{max}/\mathrm{T}$. Indeed, consider for example
time slot $0$. If the total reward could be higher than
$W_{\mathrm{EI}(\mathrm{T})}d_{max}/\mathrm{T}$, then scheduling every
$d_{max}$ time slots only, we could obtain a reward higher than
$W_{\mathrm{EI}(\mathrm{T})}d_{max}/\mathrm{T}\times\mathrm{T}/d_{max}=W_{\mathrm{EI}(\mathrm{T})}$,
which leads to a contradiction.

Since at every time slot we can lose a maximum reward of
$W_{\mathrm{EI}(\mathrm{T})}d_{max}/\mathrm{T}$, in $2d_{max}$ time slots we
lose a maximum reward of at most
$W_{\mathrm{EI}}(\mathrm{T})2\frac{d_{max}^{2}}{\mathrm{T}}$.

(a) Range of capacity constraints for $\mathrm{EI}(\mathrm{T})$

(b) Modified range for $\mathrm{\overline{EI}}(\mathrm{T})$

Figure 4. Figure 4(a) shows the time horizon for $\mathrm{EI}(\mathrm{T})$
when $d_{max}=2$. The starting point of each arrow indicates the earliest time
slot $\tau$ with arrivals that impact time slot $t$, which is indicated by the
end of an arrow, as seen in Equation 3c. For example time slot $t=2$ can have
arrivals as early as $\tau=0$. The symmetry of the problem across time slots
breaks at $t=0,1,\mathrm{T}+1,\overline{\mathrm{T}}$ (shaded). Adding terms
$x_{jk}^{\tau}$ for these time slots in (3c) allows us to create a set of
stricter constraints (5) that make $\mathrm{EI}(\mathrm{T})$ symmetrical. This
is shown in Figure 4(b)), where, e.g., now at the capacity constraint for
$t=0$, we have terms $x_{jk}^{\tau}$ from
$\tau=\mathrm{T}+1,\overline{\mathrm{T}}$.

∎

Due to Lemma 2, we can search for stationary solutions, which simplifies
$\mathrm{EI}(\mathrm{T})$ significantly. We refer to the simplified LP as the
stationary expected instance program $\mathrm{E}_{S}$, defined in (6a)-(6d)
below:

(6a) $\displaystyle\max_{\mathbf{x}}\quad$
$\displaystyle\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}\sum_{k\in\mathcal{K}_{j}}x_{jk}\quad(:=\mathrm{E}_{S})$
(6b) s.t. $\displaystyle\sum_{k\in\mathcal{K}_{j}}x_{jk}\leq 1,\quad\forall
j\in[\mathrm{J}],$ (6c)
$\displaystyle\sum_{j=1}^{\mathrm{J}}\sum_{\tau=0}^{d_{j}}\lambda_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}x_{jk}\leq\mathrm{C}_{\ell},\quad\forall\ell\in\mathcal{L},$
(6d) $\displaystyle x_{jk}\geq 0,\quad\forall j\in[\mathrm{J}],\forall
k\in\mathcal{K}_{j}.$

Note that compared to $\mathrm{EI}(\mathrm{T})$, the capacity constraints (6c)
are much simpler, as constraints (3c) in $\mathrm{EI}(\mathrm{T})$ for
$t\in[d_{max},\mathrm{T}]$ become identical, and the constraints for
$t<d_{max}$ or $t>\mathrm{T}$ are weaker. Therefore, we only need to keep the
constraints (3c) for only one $t\in[d_{max},\mathrm{T}]$, e.g., for
$t=d_{max}$.

Step 3. The simplified LP $\mathrm{E}_{S}$ still has exponentially many
variables (there are exponentially many route-schedules $k$). To resolve this,
we work with the flow-based reformulation of $\mathrm{E}_{S}$ through the
forwarding variables defined in Section 3 for Algorithm 1. Analogous
reformulations have been utilized to a great extent in network flow problems
(Magnanti and Orlin, 1993) (for more details, refer to Appendix D). The
conversion results in the flow-based LP $\mathrm{F}_{S}$, defined earlier in
(2a)-(2f), but with $\epsilon=0$. To see that, recall that $x_{jk}$ is the
probability of scheduling packet-type $j$ over route-schedule $k$, and
forwarding variable $f_{j\ell}^{\tau}$ can be interpreted as the fraction of
type-$j$ packets scheduled over link $\ell$ at age $\tau$. Then, (6a) and (6b)
are equivalent to (2a) and (2b), respectively, since $\sum_{k\in
K_{j}}x_{jk}\equiv\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0}$, where,
recall that $\mathtt{Out}(s_{j})$ denotes the outgoing links of $s_{j}$ in
$\overline{\mathcal{G}}$. Similarly, capacity constraint (6c) is equivalent to
(2e) with $\epsilon=0$. The remaining constraints (2c), (2d), (2f) are the
flow conservation laws. At the end, similarly to network flow problems, we can
map the solution $\\{f_{j\ell}^{\tau}\\}$ to a solution for $\\{x_{jk}\\}$ due
to the flow decomposition theorem (Goldberg and Rao, 1998). Note that in
$\mathrm{F}_{S}$, we have reduced the capacity in constraint (2e) by a factor
$(1+\epsilon)$. It is easy to verify that solving $\mathrm{F}_{S}$, with
$\epsilon>0$, results in at most a factor $\frac{1}{1+\epsilon}$ reduction in
the optimal objective value, and consequently,
$W_{\mathrm{F}_{S}}\geq\frac{W_{\mathrm{E}_{S}}}{1+\epsilon}$.

Step 4. In this step, we tie our analysis together as follows. If we could
schedule the packets successfully according to forwarding probabilities
$\\{f_{j\ell}^{\tau\star}\\}$ we would obtain total expected reward close to
$\mathrm{EI}(\mathrm{T})$ (Lemma 2 and Step 3). However, scheduling according
to these probabilities might violate the capacity constraints of
$\mathrm{RI}(\mathrm{T})$. We argue that due to the gap between the expected
traffic and the actual capacity on a link, imposed in (2e) by $\epsilon>0$,
the probability that traffic surpasses any link’s capacity is significantly
reduced. In fact, we obtain an “almost feasible” solution to
$\mathrm{RI}(\mathrm{T})$, if the capacity is greater than a certain
threshold. Then, by dropping packets attempted to be scheduled on a link with
full capacity, we obtain a feasible solution to $\mathrm{RI}(\mathrm{T})$.
Further, the drop probability of any packet is small (i.e. bounded by
$\epsilon$) as each link over which the packet might be transmitted is
unlikely to be at full capacity. The above result is stated formally in Lemma
3.

###### Lemma 0.

Assume, for each $j\in[\mathrm{J}]$, $\\{a_{j}^{t},t\geq 0\\}$ are i.i.d.
Bernoulli (or Binomial) random variables. Then using forwarding variables
$\\{f_{j\ell}^{\tau\star}\\}$ for scheduling packets (Lines 1-1 of Algorithm
1), the probability of a packet being dropped is at most $\epsilon$, if
$\mathrm{C}_{\min}\geq
2\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log\frac{L}{\epsilon}$.

###### Proof of Lemma 3.

In the proof, we leverage the following concentration inequality from (Chung
and Lu, 2006).

###### Lemma 0.

Suppose $X_{1},X_{2},\cdots,X_{n}$ are independent, but not necessarily
identically distributed random variables, with $X_{i}\leq B$, for all
$i\in[n]$. Let $X:=\sum_{i=1}^{n}X_{i}$ and
$S:=\sum_{i=1}^{n}\mathbb{E}\left[X_{i}^{2}\right]$. Then

(7) $\Pr{[X\geq\mathbb{E}\left[X\right]+\lambda]}\leq
e^{-\frac{\lambda^{2}}{2S+2B\lambda/3}}.$

Using Lemma 4, we prove that exceeding the capacity of any link happens with
small probability, which will then allow us to bound the probability of a
packet being dropped. First, in accordance with the notation in Lemma 4, we
define the variables:

$X_{j\tau}=\mathds{1}(\text{a type-$j$ packet arrives at time $\tau$ and is
scheduled at $\ell$ at time $t$})$

where $\mathds{1}(\cdot)$ is the indicator function. The total capacity
consumption on $\ell$ at time $t$ can be written as:

(8) $X:=\sum_{j=1}^{\mathrm{J}}\sum_{\tau=t-d_{max}}^{t}X_{j\tau},$

where we assume the case of $d_{max}\leq t\leq\mathrm{T}$ for simplicity, but
the same holds for all $t$ (since we have fewer terms in the other two cases,
the probability of packet drops will be lower). Note that

(9)
$\mathbb{E}\left[(X_{j\tau})^{2}\right]=\mathbb{E}\left[X_{j\tau}\right]=\lambda_{j}f_{j\ell}^{(t-\tau)\star},$

where the first equality is due to $X_{j\tau}\in\\{0,1\\}$ and the second
equality follows by a simple calculation, using the probability of a type-$j$
packet arriving at $\tau$ and being scheduled with total probability
$f_{j\ell}^{(t-\tau)\star}$ over $\ell$, where the optimal forwarding
probabilities $\mathbf{f}^{\star}$ found through $\mathrm{F}_{S}$ are used.
Taking an expectation on the expression in (8) and utilizing (9), and using
the fact that variables $\mathbf{f}$ are a solution to $\mathrm{F}_{S}$ and
thus they satisfy (2e), we have
$S=\mathbb{E}\left[X\right]\leq\frac{\mathrm{C}_{\ell}}{1+\epsilon}:=S_{u}$.
Then, due to the i.i.d. Bernoulli (or Binomial) assumption, all random
variables in the right-hand-side of (8) are independent, and therefore we can
apply Lemma 4, with $B=1$, for $\lambda=\epsilon
S_{u}=\frac{\epsilon\mathrm{C}_{\ell}}{1+\epsilon}$, and using $S\leq S_{u}$.
For the exponent in (7) we have:

$-\frac{\lambda^{2}}{2S+2\lambda/3}=-\frac{\epsilon^{2}S_{u}^{2}}{2S_{u}+2\epsilon
S_{u}/3}\leq-\frac{\epsilon^{2}S_{u}}{2+2\epsilon/3}\leq-\frac{\epsilon^{2}S_{u}}{2(1+\epsilon)}=-\frac{\epsilon^{2}\mathrm{C}_{\ell}}{2(1+\epsilon)^{2}}$

Therefore using inequality (7) we get:

$\Pr{[X\geq\frac{\mathrm{C}_{\ell}}{1+\epsilon}+\frac{\epsilon\mathrm{C}_{\ell}}{1+\epsilon}]}=\Pr{[X\geq\mathrm{C}_{\ell}]}\leq
e^{-\frac{\mathrm{C}_{\ell}}{2}(\frac{\epsilon}{1+\epsilon})^{2}}\overset{(a)}{\leq}e^{\log\frac{\epsilon}{L}}=\epsilon/L$

where in $(a)$ we have used the stated minimum link capacity in Lemma 3. Now,
consider an arbitrary packet which is attempted to be scheduled over a route-
schedule with a maximum of $L$ link-time slot pairs. Then, the probability
that any link in the route-schedule exceeds the capacity constraint can be
bounded through a union bound that gives the statement of the Lemma, i.e., the
probability of dropping the packet is bounded by $\epsilon$.∎

We note that the proof extends to non-stationary arrival rates, and can also
be generalized to dependent packet arrivals (Lemma 3). Since packets are
dropped with probability at most $\epsilon$, the performance loss of the
algorithm is at most a multiplicative factor $\epsilon$ (in addition to factor
$1/(1+\epsilon)$ due to the scaling of the capacities in (2e), as discussed in
Step 3). Overall, the expected obtained reward is close to
$\mathrm{EI}(\mathrm{T})$, which, by Lemma 1, gives a guarantee relative to
the average optimal value of $\mathrm{RI}(\mathrm{T})$, yielding Theorem 1.
This is detailed in the following proof.

###### Proof of Theorem 1.

The expected reward in a time slot is:

$\displaystyle\mathbb{E}\left[W_{\mathrm{ALG}}\right]=\mathbb{E}\left[\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\sum_{n=1}^{a_{j}^{t}}\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt}\right]=\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\mathbb{E}\left[\sum_{n=1}^{a_{j}^{t}}\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt}\right]\overset{(a)}{=}$
$\displaystyle\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}^{t}\mathbb{E}\left[\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt}\right]\overset{(b)}{\geq}\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}^{t}\sum_{\ell\in\mathcal{O}(s_{j})}(1-\epsilon)f_{j\ell}^{0\star}\overset{(c)}{=}\frac{1-\epsilon}{1+\epsilon}\overline{W}_{\mathrm{EI}(\mathrm{T})}\overset{(d)}{\geq}$
$\displaystyle\frac{1-\epsilon}{1+\epsilon}\left(1-2\frac{d_{max}^{2}}{T}\right)W_{\mathrm{EI}(\mathrm{T})}\overset{(e)}{\geq}\frac{1-\epsilon}{1+\epsilon}\left(1-2\frac{d_{max}^{2}}{T}\right)\mathbb{E}\left[W_{\mathrm{RI}(\mathrm{T})}\right]\overset{(f)}{\geq}(1-3\epsilon)\mathbb{E}\left[W_{\mathrm{RI}(\mathrm{T})}\right],$

where $(a)$ follows due to the independence between $\\{a_{j}^{t}\\}$ and
decisions $\\{y_{jk}^{nt}\\}$, $(b)$ follows by Lemma 3 due to the forwarding
probabilities under $\mathrm{FBPF}$, $(c)$ follows due to Step 3, $(d),(e)$
follow due to Lemmas 2,1 respectively, and $(f)$ follows for
$\mathrm{T}\geq\frac{2d_{max}^{2}}{\epsilon}$. Note that the dependence on an
i.i.d Bernoulli or Binomial distribution only appears on $(b)$ and can be
lifted by extending Lemma 3. ∎

###  4.2. Proof Outline of Theorem 2 for Algorithm 2

Here we present the outline for the proof of Theorem 2 for Algorithm 2, which
is divided into 3 steps. The detailed proofs are provided in Section A.2.

Step 1. In the first step, we find the number of samples required to obtain
estimates for $\\{\lambda_{j}\\}$ within a desired multiplicative accuracy
$(1\pm\epsilon_{0})$, for some parameter $\epsilon_{0}$. This is stated in
Lemma 5.

###### Lemma 0.

For any i.i.d. distribution on $\\{a_{j}^{t}\\}$, to estimate
$\\{\lambda_{j}\\}$ such that
$\lambda_{j}(1-\epsilon_{0})\leq\hat{\lambda}_{j}\leq\lambda_{j}(1+\epsilon_{0}),\forall
j$, with probability at least $1-\epsilon$, given
$\epsilon,\epsilon_{0}\in(0,1)$, we require samples of the arrivals from
$\frac{a_{\max}^{2}}{2\lambda_{\min}^{2}\epsilon_{0}^{2}}\log\left(\frac{2\mathrm{J}}{\epsilon}\right)$
time slots.

Step 2. Towards the goal of solving $\mathrm{EI}(\mathrm{T})$ near-optimally
as in Algorithm 1, and due to the results outlined in Section 4.1, we focus on
solving $\mathrm{F}_{S}$ near optimally. As discussed in Algorithm 2 and
Algorithm 2 (Algorithm 2), we substitute $\lambda_{j}$ with estimates
$\hat{\lambda}_{j}$ of $\lambda_{j}$ in $\mathrm{F}_{S}$, to obtain a solution
$\hat{\mathbf{f}}$. Subsequently, by appropriately scaling this solution, we
obtain a near-optimal and feasible solution to $\mathrm{F}_{S}$ with high
probability. This is formalized in Lemma 6.

###### Lemma 0.

Let $\hat{\mathbf{f}}$ be the solution to $\mathrm{F}_{S}$ by using estimates
$\hat{\lambda}_{j}$, satisfying
$\lambda_{j}(1-\epsilon_{0})\leq\hat{\lambda}_{j}\leq\lambda_{j}(1+\epsilon_{0})$,
in $\mathrm{F}_{S}$. Let
$\tilde{\mathbf{f}}=(1-\epsilon_{0})\hat{\mathbf{f}}$. The scaled solution
$\tilde{\mathbf{f}}$, whose objective value we denote with $\tilde{F}_{S}$,
is, with probability at least $1-\epsilon$, a feasible solution to
$\mathrm{F}_{S}$ and

$\tilde{F}_{S}\geq(1-4\epsilon_{0})\mathrm{F}_{S}.$

Step 3. Based on Lemma 6, we can solve $F_{S}$ with progressively better
estimates $\hat{\lambda}_{j}$ in $\log(1/\epsilon)$ phases, to obtain an
increasingly better approximation in each phase for $\mathrm{RI}(\mathrm{T})$,
following the connection between $\mathrm{F}_{S}$ and
$\mathrm{RI}(\mathrm{T})$ discussed in Section 4.1. More specifically, in each
phase we have a different number of samples available and based on Lemma 5 and
Lemma 6, $\epsilon_{0}$ is decreasing across phases (and determined by the
carefully designed phase durations). Aggregating the approximations of each
phase, we obtain an overall near-optimal solution assuming sufficient time is
available and under the capacity constraints in the statement of Algorithm 2.
Refer to Section A.2 for the detailed proofs.

##  5\. Discussion

In this section, we discuss how our technique can be extended to other
distributional assumptions and discuss some complexity and implementation
improvements.

###  5.1. Extension to more general arrival distributions

General stationary distributions. The theoretical results in Theorem 1 and
Theorem 2 extend to general distributions on the arrivals, other than i.i.d.
Binomial and Bernoulli arrivals. The distribution affects the required minimum
link-capacity, and generally, improved results are obtained when packet
arrivals are less dependent on each other. Specifically, if the packet
arrivals of a packet-type in any $(d_{max}+1)$ consecutive time slots can be
partitioned into a set of groups of maximum size $D$ such that each group’s
arrivals are independent of the arrivals of other groups, then our results can
be generalized, requiring capacity $\mathrm{C}_{\min}\geq
2D\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log(L/\epsilon)$ for near-
optimal ($1-3\epsilon$)-approximation (Theorem 2 in Appendix B). Binomial and
Bernoulli arrivals correspond to $D=1$.

We remark that this extension allows a dependence between packets within the
same time slot or across different time slots. For example, if the arrivals of
a packet type are independent across time, and are at most $D$ at any time,
then the above partitioning is possible. Alternatively, packets might depend
across time. For example, when $a_{\max}=1$, with the arrivals of a packet
type at any time, dependent on the arrivals of that type in prior slots, we
obtain near-optimal performance for $\mathrm{C}_{\min}\geq
2(d_{max}+1)\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log(L/\epsilon)$,
that is, $D=(d_{max}+1)$ in this case.

The formal description of these results is provided in Appendix B.

Known non-stationary traffic. Our assumption on identically distributed
arrivals across time slots can be relaxed. In that case, due to the presence
of non-stationarity, we cannot rely on stationary forwarding probabilities
$f_{j\ell}^{\tau}$ (identified earlier in $\mathrm{F}_{S}$) to obtain a near-
optimal solution. As a result, in contrast to the i.i.d. case, the forwarding
variables will need to be time dependent, i.e., $f_{j\ell}^{\tau t}$, in
addition to the dependence on link, type and age of the packet. The modified
LP is referred to as $\mathrm{F}_{NS}$, and is presented in Section C.1, see
(18a)-(18f). For example, the objective in $\mathrm{F}_{S}$ should be modified
to describe the average total reward (rather than the average per time slot
reward) and it becomes:

$\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}^{t}\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0t}.$

Further, the link capacity constraints (2e) should be adjusted such that they
are applied for each time slot independently (as opposed to a single
constraint that captures the average traffic on the link). Note that the
modification for non-stationarity comes with a computational cost, as the
number of variables and constraints now scale also with the horizon
$\mathrm{T}$. The cost of solving $\mathrm{F}_{NS}$ with a large horizon can
be mitigated by splitting the horizon into frames, and solving
$\mathrm{F}_{NS}$ in each independent frame. A trade off between complexity
and performance arises, with larger considered frames generally yielding
better performance but with higher complexity (we provide more details in
Section C.2, e.g., Lemma 2). Following the determination of forwarding
variables $f_{j\ell}^{\tau t}$, we can proceed by probabilistically forwarding
any packet of type $j$ arriving at time $t$ according to variables
$f_{j\ell}^{\tau t}$ (similarly to Algorithm 1). We provide the general
Theorem 1 in Section C.1.

Periodic traffic distribution. In the special case of periodic traffic
distribution, the complexity of the non-stationary approach can be further
reduced such that the number of variables and constraints do not scale with
$\mathrm{T}$, but rather, scale with the length of the period of the periodic
traffic. In this case, our proof techniques outlined in Section 4 (such as
Lemma 2) can be modified to show that a near optimal solution is periodic (as
opposed to stationary). For example, for a period of length
$\widetilde{\mathrm{T}}$, the forwarding variables are $\\{f_{j\ell}^{\tau
t},t\in[\widetilde{\mathrm{T}}]\\}$ and a packet arriving at time $t^{\prime}$
must be forwarded according to variables
$\\{f_{j\ell}^{\tau(t^{\prime}\mathrm{mod}\ \widetilde{\mathrm{T}})}\\}$ . For
more details, refer to Section C.3.

Unknown arrival rates. For the non-stationary traffic when the model is not
known, our techniques can be used in conjunction with prediction methods that
estimate future arrival rates $\lambda_{j}^{t}$. In that case, the predicted
future $\lambda_{j}^{t}$ is used as input to $\mathrm{F}_{NS}$. In particular,
the performance of our methods may be characterized based on the accuracy of
the predictions. Prediction-based methods have been studied recently in (Stein
and Wei, 2023). By leveraging an extension to Lemma 6, we can obtain
performance which depends on the accuracy of the estimated arrival rates
$\widehat{\lambda_{j}^{t}}$. Further, by using frame-based methods, such as
the one proposed in Lemma 2 (Section C.2), we only require predicting the
arrival rates for an incoming frame.

###  5.2. Complexity and implementation

Reducing the number of used route-schedules. In our algorithms, we described a
hop-by-hop randomization, in order to convert the forwarding probabilities to
a particular route-schedule. However, depending on the design needs, other
conversions might be preferred and are available. For example, an iterative
procedure based on (Raghavan and Tompson, 1987) can be used to find a fixed
set of route-schedules, over which the source of the packet can randomize over
for scheduling the packet. Specifically, each packet type will have a maximum
of $d_{max}|\overline{\mathcal{L}}|$ route-schedules over which the packets of
this type are scheduled. This is in contrast to the hop-by-hop approach that
can potentially choose any route-schedule. In certain scenarios, the former
method might be preferred. More details on this alternate process are provided
in Section D.1.

Solving the LP $\mathrm{F}_{S}$. Recall that $\mathrm{F}_{S}$ (2) is an LP
which has $O(d_{max}|\overline{\mathcal{L}}|\mathrm{J})$ number of variables
and constraints and therefore can be solved efficiently in polynomial time.
Following our results, it is interesting to seek even more efficient
implementations of $\mathrm{F}_{S}$ or online solutions of $\mathrm{F}_{S}$,
through, for example distributed primal-dual techniques, or multicommodity-
flow-inspired solutions (Magnanti and Orlin, 1993).

##  6\. Simulation Results

In this section, we report simulation results using both synthetic data and
real network data.

###  6.1. Evaluation using synthetic traffic

Comparisons between algorithms. First, we consider network $\mathcal{G}_{1}$
shown in Figure 5(a). In our first experiment, we examine traffic generated
using $10$ packet types, with source $s_{j}$ and destination $z_{j}$ for each
type $j$ chosen randomly in the network, weights $w_{j}$ chosen randomly in
$(0,1)$, and arrival rates $\lambda_{j}$ uniformly chosen in $(0,25)$. The
arriving packets have a deadline of $10$, and the time horizon is
$\mathrm{T}=5000$.

(a) Simple network $\mathcal{G}_{1}$

(b) IBM network $\mathcal{G}_{IBM}$

Figure 5. Two different networks used in simulations.

We evaluate three variants of our algorithms: $\mathrm{FBPF}$ (Algorithm 1),
and two variants of $\mathrm{DLPF}$ (Algorithm 2), the base method which
divides the horizon into exponentially increasing phases (denoted with
$\mathrm{DLPF}$-EXP), and a method that divides the horizon into phases of
fixed duration of $500$ slots ($\mathrm{DLPF}$-500). With the goal of
verifying the near-optimal behavior of our algorithms as $\mathrm{C}_{\min}$
increases, we additionally provide comparisons with the recent worst-case-
based algorithm GLS-FP (Gu et al., 2021) (Fastest-Path variant of GLS), which
schedules packets greedily on the “fastest path“ (the route-schedule that
delivers the packet as early as possible) from source to destination, when
enough capacity is available, otherwise it drops the packet. This is the
algorithm of choice in the simulations of (Gu et al., 2021), which in the
simulations therein outperforms works such as the one in (Deng et al., 2019).

To evaluate our algorithms with respect to the optimal, we used an upper bound
on the optimal value of $\mathrm{RI}(\mathrm{T})$ based on the optimal value
of $\mathrm{F}_{S}$ ((2a)-(2f)) (multiplied by
$(1-2d_{max}^{2}/\mathrm{T})^{-1}$ due to (4) and our earlier analysis,
outlined in Section 4.1). Hence, we can obtain a lower bound on the
approximation ratio of each algorithm, which would be the metric in our
comparisons here.

In our first experiment, we study the impact of link capacity with fixed
traffic. We considered an equal link capacity
$\mathrm{C}_{\ell}=\mathrm{C}_{\min}$ for all the links in the network and
varied $\mathrm{C}_{\min}$ between $1$ and $25$. The results are shown in
Figure 6(a). We observe that all our algorithms outperform GLS-FP
significantly, and even for very small capacities, obtain more than a fraction
$0.8$ of the optimal value.

In the next experiment, we study the impact of increasing capacity
$\mathrm{C}_{\min}$ while also scaling the traffic intensity proportionally
(to avoid simplification of the problem). We do so by increasing the number of
packet types with the capacity, with $\mathrm{J}=\\{5,10,\cdots,40\\}$ and
$\mathrm{C}_{\min}=\\{1,4,\cdots,25\\}$ (as earlier). The results are depicted
in Figure 6(b). We observe that our algorithms improve the performance with
the scaling of capacity and traffic regardless, whereas GLS-FP does not
improve as much as in Figure 6(a). This is in line with our analysis that
suggests an improvement of the approximation ratio for larger capacities.

The good performance of our algorithms was further validated in an extensive
set of topologies, deadlines and traffic configurations. Refer to Appendix E
for additional simulation results.

Performance over different arrival distributions. As discussed in Section 5.1
(with more details in Appendix B), our methods extend to general packet
arrival distributions, with the characteristics of the distribution only
affecting the required $\mathrm{C}_{\min}$ to obtain near-optimal performance.
Here, we validate the good performance of $\mathrm{FBPF}$, under three
different distributions, for identical arrival rates in each case. In
particular, we consider an i.i.d. process with arrivals at each time slot and
for each packet type from the standard Binomial distribution, a Poisson
distribution, and a scaled Bernoulli distribution which is either $0$ or
$a_{j}$ for packet type $j$ (all with equal $\lambda_{j}$). In particular, we
note that i.i.d. arrivals from the Binomial distribution result in the minimum
dependency degree $D=1$, whereas those from Scaled Bernoulli have maximum
dependency degree, with $D=a_{\max}$. We use IBM’s network $\mathcal{G}_{IBM}$
(Zoo, 2011), shown in Figure 7(a) (similar results are obtained for
$\mathcal{G}_{1}$), which consists of $18$ nodes and $20$ links. Here, we
consider higher arrival rates compared to the previous simulations, randomly
selected in $(0,150)$ for each $j$, as we seek to study the performance for
larger capacities in order to verify that approximation ratios are improved
over all distributions as $\mathrm{C}_{\min}$ increases. The approximation
ratio bounds are shown in Figure 6(c). We notice that the best performance is
obtained under the Binomial distribution, in line with our theoretical results
($D=1$, Theorem 2). The Poisson distribution also exhibits good performance.
The worst performance is obtained (as expected) for the Scaled Bernoulli
$(D=a_{\max})$.

###  6.2. Evaluation using real traffic

(a) Varying capacity

(b) Varying capacity and traffic

(c) Three packet arrival distributions

Figure 6. Comparison of algorithms in terms of obtained approximation ratios.

In this section, we focus on understanding the performance of our algorithms
in real networks with real packet arrival traces. First, we consider IBM’s
network $\mathcal{G}_{IBM}$ (Zoo, 2011), shown in Figure 7(a). To simulate the
packet arrivals, we use a trace of IP traffic collected at the University of
Cauca (Rojas et al., 2019). We focus on the arrivals during the $10$ minutes
of the dataset. As nodes in networks such as $\mathcal{G}_{IBM}$ typically
serve multiple clients (e.g. routers of an ISP, serving the traffic of users),
we assign the traffic of roughly $1000$ IPs (all the IPs in the trace) to the
$18$ nodes of the network. Based on this assignment, we obtain a traffic flow
for each ordered pair of nodes in $\mathcal{G}_{IBM}$. For example in $18$
nodes, there are $18\times 17$ source-destination pairs. Many of these pairs,
following the random assignment of IPs result in little to no traffic. As a
result, we focus on the $30$ source-destination pairs with the highest
traffic. We illustrate three of the final traces of packet arrivals for
different source-destination pairs with different characteristics in Figure 7.
We observe that all three traces have clear non-stationarity and are
characterized by sudden bursts of packets.

(a) Trace $1$

(b) Trace $2$

(c) Trace $3$

Figure 7. Three different traces of packet arrivals for three source-
destination pairs in $\mathcal{G}_{IBM}$.

In our simulations, we assign a random reward uniformly chosen between $(0,1)$
for each source-destination pair and a deadline $10$ time slots to each
packet. To further evaluate the impact of traffic intensity and the number of
flows, we simulate two cases. In one case, we maintain all top $30$ source-
destination pairs, whereas in the second experiment, we maintain a random
subset of $15$ source-destination pairs.

As the traffic is non-stationary, we focus our study on DLPF-30, DLPF-100, two
variants of our dynamic algorithm DLPF, that update their estimates every
$30$, and $100$ time slots respectively (our preliminary experiments included
phase lengths larger than $100$, which were found to attain a lower reward
compared to phase lengths of $30$ and $100$). We compare the algorithms in
terms of their average reward per time slot. The final results for the two
experiments are shown in Figure 8(a) and Figure 8(b). As we can see, the
algorithms maintain a significant gain over GLS-FP, and the gain grows as the
traffic intensity increases and the problem becomes more challenging.

Finally, we verified the good performance of our methods over an additional
network. We used the Hibernia Atlantic (Canada) network from (Zoo, 2011),
which is a network of smaller size compared to $\mathcal{G}_{IBM}$. Under
similar simulation configurations but over the new topology (and the mapping
of the IPs on the new nodes), we simulated the three methods, showing the
results in Figure 8(c). Both of our methods are outperforming GLS-FP in this
topology as well.

(a) $\mathcal{G}_{IBM}$ with 15 source-destination pairs

(b) $\mathcal{G}_{IBM}$ with 30 source-destination pairs

(c) Hibernia Atlantic Network with 15 source-destination pairs

Figure 8. Performance comparison (in terms of average reward per time slot)
using real traffic on $\mathcal{G}_{IBM}$ (Figure 8(a), Figure 8(b)) and the
Hibernia Atlantic Network (Figure 8(c)).

##  7\. Conclusions

This paper shows that it is possible to design algorithms that provide near-
optimal $(1-\epsilon)$ approximation ratio for the problem of scheduling
deadline-constrained multihop traffic in the stochastic setting. Our result
requires minimum link capacity
$\mathrm{C}_{\min}=\Omega(\log(L/\epsilon)/\epsilon^{2})$ to guarantee such
performance. This is a significant improvement over the prior results in the
worst-case traffic setting. Our techniques are general and can be applied
under different distributional assumptions. An interesting future work is to
investigate more efficient and distributed ways of solving the LP
$\mathrm{F}_{S}$, e.g., by using distributed primal-dual methods.

## References

  * (1)
  * Agrawal and Devanur (2014) Shipra Agrawal and Nikhil R Devanur. 2014.  Fast algorithms for online stochastic convex programming. In _Proceedings of the twenty-sixth annual ACM-SIAM symposium on Discrete algorithms_. SIAM, ACM, 1405–1424.
  * Andrews et al. (2000) Matthew Andrews, Antonio Fernández, Mor Harchol-Balter, Tom Leighton, and Lisa Zhang. 2000.  General dynamic routing with per-packet delay guarantees of O (distance+ 1/session rate).  _SIAM J. Comput._ 30, 5 (2000), 1594–1623.
  * Andrews and Zhang (1999) Matthew Andrews and Lisa Zhang. 1999.  Packet routing with arbitrary end-to-end delay requirements. In _Proceedings of the thirty-first annual ACM symposium on Theory of Computing_. 557–565.
  * Banerjee and Freund (2020) Siddhartha Banerjee and Daniel Freund. 2020.  Uniform loss algorithms for online stochastic decision-making with applications to bin packing. In _Abstracts of the 2020 SIGMETRICS/Performance Joint International Conference on Measurement and Modeling of Computer Systems_. 1–2.
  * Baruah et al. (1992) Sanjoy Baruah, Gilad Koren, Decao Mao, Bhubaneswar Mishra, Arvind Raghunathan, Louis Rosier, Dennis Shasha, and Fuxing Wang. 1992.  On the competitiveness of on-line real-time task scheduling.  _Real-Time Systems_ 4 (1992), 125–144.
  * Bhattacharya et al. (1997) Partha P Bhattacharya, Leandros Tassiulas, and Anthony Ephremides. 1997.  Optimal scheduling with deadline constraints in tree networks.  _IEEE Trans. Automat. Control_ 42, 12 (1997), 1703–1705.
  * Chin et al. (2006) Francis YL Chin, Marek Chrobak, Stanley PY Fung, Wojciech Jawor, Jiří Sgall, and Tomáš Tichỳ. 2006.  Online competitive algorithms for maximizing weighted throughput of unit jobs.  _Journal of Discrete Algorithms_ 4, 2 (2006), 255–276.
  * Chung and Lu (2006) F.R.K. Chung and L. Lu. 2006.  _Complex Graphs and Networks_.  American Mathematical Society.
  * Deng et al. (2019) Han Deng, Tao Zhao, and I-Hong Hou. 2019.  Online routing and scheduling with capacity redundancy for timely delivery guarantees in multihop networks.  _IEEE/ACM Transactions on Networking_ 27, 3 (2019), 1258–1271.
  * Devanur et al. (2019) Nikhil R. Devanur, Kamal Jain, Balasubramanian Sivan, and Christopher A. Wilkens. 2019.  Near Optimal Online Algorithms and Fast Approximation Algorithms for Resource Allocation Problems.  _J. ACM_ 66, 1 (2019).
  * Fan et al. (2019) Keke Fan, Ying Wang, Junhua Ba, Wenjing Li, and Qi Li. 2019.  An approach for energy efficient deadline-constrained flow scheduling and routing. In _2019 IFIP/IEEE symposium on integrated network and service management (IM)_. IEEE, 469–475.
  * Fountoulakis et al. (2023) Emmanouil Fountoulakis, Themistoklis Charalambous, Anthony Ephremides, and Nikolaos Pappas. 2023.  Scheduling Policies for AoI Minimization with Timely Throughput Constraints.  _IEEE Transactions on Communications_ (2023).
  * Gallego and Van Ryzin (1994) Guillermo Gallego and Garrett Van Ryzin. 1994.  Optimal dynamic pricing of inventories with stochastic demand over finite horizons.  _Management science_ 40, 8 (1994), 999–1020.
  * Goldberg and Rao (1998) Andrew V Goldberg and Satish Rao. 1998.  Beyond the flow decomposition barrier.  _Journal of the ACM (JACM)_ 45, 5 (1998), 783–797.
  * Gu et al. (2021) Yan Gu, Bo Liu, and Xiaojun Shen. 2021.  Asymptotically Optimal Online Scheduling With Arbitrary Hard Deadlines in Multi-Hop Communication Networks.  _IEEE/ACM Transactions on Networking_ (2021).
  * Hoeffding (1963) Wassily Hoeffding. 1963.  Probability inequalities for sums of bounded random variables.  _Journal of the American statistical association_ 58, 301 (1963), 13–30.
  * Hou (2013) I-Hong Hou. 2013.  Scheduling heterogeneous real-time traffic over fading wireless channels.  _IEEE/ACM Transactions on Networking_ 22, 5 (2013), 1631–1644.
  * Hou et al. (2009) I.-H. Hou, V. Borkar, and P. R. Kumar. 2009.  A Theory of QoS for Wireless. In _IEEE INFOCOM 2009_. 486–494.
  * Jeż (2013) Łukasz Jeż. 2013.  A universal randomized packet scheduling algorithm.  _Algorithmica_ 67, 4 (2013), 498–515.
  * Jiang et al. (2020) Jiashuo Jiang, Xiaocheng Li, and Jiawei Zhang. 2020.  Online stochastic optimization with wasserstein based non-stationarity.  _arXiv preprint arXiv:2012.06961_ (2020).
  * Kang et al. (2013) Xiaohan Kang, Weina Wang, Juan José Jaramillo, and Lei Ying. 2013.  On the performance of largest-deficit-first for scheduling real-time traffic in wireless networks. In _Proceedings of the fourteenth ACM international symposium on Mobile ad hoc networking and computing_. 99–108.
  * Kang et al. (2014) Xiaohan Kang, Weina Wang, Juan José Jaramillo, and Lei Ying. 2014.  On the performance of largest-deficit-first for scheduling real-time traffic in wireless networks.  _IEEE/ACM Transactions on Networking_ 24, 1 (2014), 72–84.
  * Kesselheim et al. (2014) Thomas Kesselheim, Andreas Tönnis, Klaus Radke, and Berthold Vöcking. 2014.  Primal beats dual on online packing LPs in the random-order model. In _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_. 303–312.
  * Lee (2021) Jeonghwa Lee. 2021.  Generalized Bernoulli process with long-range dependence and fractional binomial distribution.  _Dependence Modeling_ 9, 1 (2021), 1–12.
  * Leighton et al. (1994) Frank Thomson Leighton, Bruce M Maggs, and Satish B Rao. 1994.  Packet routing and job-shop scheduling in O (congestion+ dilation) steps.  _Combinatorica_ 14, 2 (1994), 167–186.
  * Li et al. (2009) Hongkun Li, Yu Cheng, Chi Zhou, and Weihua Zhuang. 2009.  Minimizing end-to-end delay: A novel routing metric for multi-radio wireless mesh networks. In _IEEE INFOCOM 2009_. IEEE, 46–54.
  * Li and Eryilmaz (2012) Ruogu Li and Atilla Eryilmaz. 2012.  Scheduling for end-to-end deadline-constrained traffic with reliability requirements in multihop networks.  _IEEE/ACM Transactions on Networking_ 20, 5 (2012), 1649–1662.
  * Liu and Layland (1973) Chung Laung Liu and James W Layland. 1973.  Scheduling algorithms for multiprogramming in a hard-real-time environment.  _Journal of the ACM (JACM)_ 20, 1 (1973), 46–61.
  * Liu et al. (2019) Xin Liu, Weichang Wang, and Lei Ying. 2019.  Spatial–temporal routing for supporting end-to-end hard deadlines in multi-hop networks.  _Performance Evaluation_ 135 (2019), 102007.
  * Magnanti and Orlin (1993) Thomas L Magnanti and James B Orlin. 1993.  Network Flows.  _PHI, Englewood Cliffs NJ_ (1993).
  * Mao et al. (2014) Zhoujia Mao, Can Emre Koksal, and Ness B Shroff. 2014.  Optimal online scheduling with arbitrary hard deadlines in multihop communication networks.  _IEEE/ACM Transactions on Networking_ 24, 1 (2014), 177–189.
  * Popovski et al. (2022) Petar Popovski, Federico Chiariotti, Kaibin Huang, Anders E Kalør, Marios Kountouris, Nikolaos Pappas, and Beatriz Soret. 2022.  A perspective on time toward wireless 6G.  _Proc. IEEE_ 110, 8 (2022), 1116–1146.
  * Raghavan and Tompson (1987) Prabhakar Raghavan and Clark D Tompson. 1987.  Randomized rounding: a technique for provably good algorithms and algorithmic proofs.  _Combinatorica_ 7, 4 (1987), 365–374.
  * Rojas et al. (2019) Juan Sebastián Rojas, Alvaro Rendon, and Juan Carlos Corrales. 2019.  Consumption behavior analysis of over the top services: Incremental learning or traditional methods?  _IEEE Access_ 7 (2019), 136581–136591.
  * Singh et al. (2013) Rahul Singh, I-Hong Hou, and PR Kumar. 2013.  Pathwise performance of debt based policies for wireless networks with hard delay constraints. In _52nd IEEE Conference on Decision and Control_. IEEE, 7838–7843.
  * Singh and Kumar (2018) Rahul Singh and PR Kumar. 2018.  Throughput optimal decentralized scheduling of multihop networks with end-to-end deadline constraints: Unreliable links.  _IEEE Trans. Automat. Control_ 64, 1 (2018), 127–142.
  * Srinivasan and Teo (1997) Aravind Srinivasan and Chung-Piaw Teo. 1997.  A constant-factor approximation algorithm for packet routing, and balancing local vs. global criteria. In _Proceedings of the twenty-ninth annual ACM symposium on Theory of computing_. 636–643.
  * Stein and Wei (2023) Clifford Stein and Hao-Ting Wei. 2023.  Learning-Augmented Online Packet Scheduling with Deadlines.  _arXiv preprint arXiv:2305.07164_ (2023).
  * Sun et al. (2021) Jingzhou Sun, Lehan Wang, Zhiyuan Jiang, Sheng Zhou, and Zhisheng Niu. 2021.  Age-optimal scheduling for heterogeneous traffic with timely throughput constraints.  _IEEE Journal on Selected Areas in Communications_ 39, 5 (2021), 1485–1498.
  * Sun et al. (2020) Rui Sun, Xinshang Wang, and Zijie Zhou. 2020.  Near-optimal primal-dual algorithms for quantity-based network revenue management.  _arXiv preprint arXiv:2011.06327_ (2020).
  * Tassiulas and Ephremides (1992) L Tassiulas and A Ephremides. 1992.  Stability properties of constrained queueing systems and scheduling policies for maximum throughput in multihop radio networks.  _IEEE Trans. Automat. Control_ 37, 12 (1992), 1936–1948.
  * Tsanikidis and Ghaderi (2020) Christos Tsanikidis and Javad Ghaderi. 2020.  On the power of randomization for scheduling real-time traffic in wireless networks. In _IEEE INFOCOM 2020_. 59–68.
  * Tsanikidis and Ghaderi (2021) Christos Tsanikidis and Javad Ghaderi. 2021.  Randomized Scheduling of Real-Time Traffic in Wireless Networks Over Fading Channels. In _IEEE INFOCOM 2021_. 1–10.
  * Tsanikidis and Ghaderi (2022) Christos Tsanikidis and Javad Ghaderi. 2022.  Online scheduling and routing with end-to-end deadline constraints in multihop wireless networks. In _Proceedings of the Twenty-Third International Symposium on Theory, Algorithmic Foundations, and Protocol Design for Mobile Networks and Mobile Computing_. 11–20.
  * Wang et al. (2014) Lin Wang, Fa Zhang, Kai Zheng, Athanasios V Vasilakos, Shaolei Ren, and Zhiyong Liu. 2014.  Energy-efficient flow scheduling and routing with hard deadlines in data center networks. In _2014 IEEE 34th International Conference on Distributed Computing Systems_. IEEE, 248–257.
  * Wang et al. (2011) Qing Wang, Pingyi Fan, Dapeng Oliver Wu, and Khaled Ben Letaief. 2011.  End-to-end delay constrained routing and scheduling for wireless sensor networks. In _2011 IEEE International Conference on Communications (ICC)_. IEEE, 1–5.
  * Zoo (2011) Internet Topology Zoo. 2011.  IBM & Hibernia Networks.  http://www.topology-zoo.org/dataset.html.  Accessed: 2023-02-1.

##  Appendix A Remaining Proofs

###  A.1. Proofs for Algorithm 1

###### Proof of Lemma 1.

Consider possibly non-stationary, general packet arrival distribution, and
without loss of generality assume $\lambda_{j}^{t}>0$. Further, consider the
following solution for the expected instance $\mathrm{EI}(\mathrm{T})$
((3a)-(3d)):
$x_{jk}^{t}=\mathbb{E}\left[\sum_{n=1}^{a_{j}^{t}}y_{jk}^{nt\star}\right]/{\lambda_{j}^{t}}$
, where, $y_{jk}^{nt\star}$ is the solution to the random instance
$\mathrm{RI}(\mathrm{T})$, and the expectation is taken with respect to the
randomness in the definition of $\mathrm{RI}(\mathrm{T})$. This solution is
feasible for $\mathrm{EI}(\mathrm{T})$. Indeed, first, note that constraints
(3b) are satisfied:

$\sum_{k\in\mathcal{K}_{j}}x_{jk}^{t}=\sum_{k\in\mathcal{K}_{j}}\frac{\mathbb{E}\left[\sum_{n=1}^{a_{j}^{t}}y_{jk}^{nt\star}\right]}{\lambda_{j}^{t}}=\frac{\mathbb{E}\left[\sum_{n=1}^{a_{j}^{t}}\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt\star}\right]}{\lambda_{j}^{t}}\overset{(\alpha)}{\leq}\frac{\mathbb{E}\left[a_{j}^{t}\right]}{\lambda_{j}^{t}}=1,$

where in $(\alpha)$, we used inequality
$\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt\star}\leq 1$ from Equation 1b, as
$y_{jk}^{nt\star}$ must be feasible for each random program
$\mathrm{RI}(\mathrm{T})$. Further, for the capacity constraint (3c), we have

(10)
$\displaystyle\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}\sum_{j=1}^{\mathrm{J}}\lambda_{j}^{t}\sum_{k:k_{t-\tau}=\ell}x_{jk}^{\tau}=\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}\sum_{j=1}^{\mathrm{J}}\sum_{k:k_{t-\tau}=\ell}\lambda_{j}^{t}x_{jk}^{\tau}$
(11)
$\displaystyle=\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}\sum_{j=1}^{\mathrm{J}}\sum_{k:k_{t-\tau}=\ell}\lambda_{j}^{t}\frac{\mathbb{E}\left[\sum_{n=1}^{a_{j}^{\tau}}y_{jk}^{n\tau\star}\right]}{\lambda_{j}^{t}}$
(12)
$\displaystyle=\mathbb{E}\left[\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}\sum_{j=1}^{\mathrm{J}}\sum_{k:k_{t-\tau}=\ell}\sum_{n=1}^{a_{j}^{\tau}}y_{jk}^{n\tau\star}\right]\overset{(a)}{\leq}\mathbb{E}\left[\mathrm{C}_{\ell}\right]=\mathrm{C}_{\ell},$

where in $(a)$ we used (1c) due to the feasibility of $y_{jk}^{nt\star}$.
Finally, the non-negativity constraints (3d) are trivially satisfied. The
objective value of the proposed solution $x_{jk}^{t}$ is then:

$\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}^{t}\sum_{k\in\mathcal{K}_{j}}x_{jk}^{t}=\mathbb{E}\left[\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\sum_{n=1}^{a_{j}^{t}}\sum_{k\in\mathcal{K}_{j}}y_{jk}^{nt\star}\right]=\mathbb{E}\big{[}W_{\mathrm{RI}(\mathrm{T})}^{\star}\big{]}.$

Therefore, we proposed a feasible solution to the expected instance which has
objective value $\mathbb{E}\big{[}W_{\mathrm{RI}(\mathrm{T})}^{\star}\big{]}$.
This implies that for the optimal objective value
$W_{\mathrm{EI}}(\mathrm{T})$ we have:
$W_{\mathrm{EI}(\mathrm{T})}\geq\mathbb{E}\big{[}W_{\mathrm{RI}(\mathrm{T})}^{\star}\big{]}$.
∎

###  A.2. Proofs for Algorithm 2

###### Proof of Lemma 5.

For a packet type $j$, let $A_{j}$ denote the event of estimating within
sufficient accuracy its probability:
$\\{\lambda_{j}(1-\epsilon_{0})\leq\hat{\lambda}_{j}\leq\lambda_{j}(1+\epsilon_{0})\\}$.
We will show that $P(\bigcap_{j=1}^{\mathrm{J}}A_{j})\geq 1-\epsilon$, or
equivalently $P(\bigcup_{j=1}^{\mathrm{J}}A_{j}^{c})\leq\epsilon$. Using the
Hoeffding bound (Hoeffding, 1963) we have for the error probability for a
fixed $j$:

$\Pr(A_{j}^{c})=\Pr(\lambda_{j}(1-\epsilon_{0})\leq\hat{\lambda}_{j}\leq\lambda_{j}(1+\epsilon_{0}))^{c}\leq
2\exp\left(-2\lambda_{j}^{2}\epsilon^{2}_{0}T/a_{\max}^{2}\right)\overset{(a)}{\leq}\frac{\epsilon}{\mathrm{J}},$

where $(a)$ holds for

$\mathrm{T}\geq\frac{a_{\max}^{2}}{\lambda_{\min}^{2}\epsilon^{2}_{0}2}\log\left(\frac{2\mathrm{J}}{\epsilon}\right).$

The result is obtained by union bound over all $j$. ∎

###### Proof of Lemma 6.

The proof is conditioned on the assumption of accurate estimates
$\\{\lambda_{j}(1-\epsilon_{0})\leq\hat{\lambda}_{j}\leq\lambda_{j}(1+\epsilon_{0})\\}$
which holds with probability at least $(1-\epsilon)$ based on Lemma 5. In this
case, first, consider the solution $\mathbf{x}^{\star}$ to $\mathrm{E}_{S}$
(identified from $\mathbf{f}^{\star}$, with $\epsilon=0$). Let us refer to the
modification of optimization $\mathrm{E}_{S}$ with probabilities
$\hat{\lambda}_{j}$ as $\hat{E}_{S}$. Then
$\frac{\mathbf{x}^{\star}}{1+\epsilon_{0}}$ is feasible for $\hat{E}_{S}$.
Indeed for the capacity constraints we have:

$\sum_{\tau=0}^{d_{max}}\sum_{j=1}^{\mathrm{J}}\hat{\lambda}_{j}w_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\frac{x_{jk}^{\star}}{1+\epsilon_{0}}\overset{(a)}{\leq}\sum_{\tau=0}^{d_{max}}\sum_{j=1}^{\mathrm{J}}\lambda_{j}w_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}x_{jk}^{\star}\overset{(b)}{\leq}\mathrm{C}_{\ell}$

where $(a)$ follows from Lemma 5, and $(b)$ follows from the feasibility of
$\mathbf{x}^{\star}$ for $\mathrm{E}_{S}$. The remaining constraints follow
trivially as they do not depend on $\lambda_{j}$. Since
$\frac{\mathbf{x}^{\star}}{1+\epsilon_{0}}$ is a feasible solution to
$\hat{E}_{S}$, we have, due to the optimality of $\mathbf{\hat{x}}$ for
$\hat{E}_{S}$:

(13)
$\sum_{j=1}^{\mathrm{J}}w_{j}\hat{\lambda}_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\hat{x}_{jk}\geq\sum_{j=1}^{\mathrm{J}}w_{j}\hat{\lambda}_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\frac{x_{jk}^{\star}}{1+\epsilon_{0}}$

We proceed by also showing that
$\mathbf{\tilde{x}}=(1-\epsilon_{0})\mathbf{\hat{x}}$ is feasible for
$\mathrm{E}_{S}$. Indeed:

$\sum_{\tau=0}^{d_{max}}\sum_{j=1}^{\mathrm{J}}\lambda_{j}w_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}(1-\epsilon_{0})\hat{x}_{jk}\leq\sum_{\tau=0}^{d_{max}}\sum_{j=1}^{\mathrm{J}}\hat{\lambda}_{j}w_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\hat{x}_{jk}\leq
C_{\ell}.$

In the following, we obtain a relation between the objective value for
$\tilde{\mathbf{x}}$ and $\hat{\mathbf{x}}$, and then between
$\hat{\mathbf{x}}$ and $\mathbf{x}^{\star}$, to conclude the Lemma.

(14)
$\frac{1+\epsilon_{0}}{1-\epsilon_{0}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\tilde{x}_{jk}=(1+\epsilon_{0})\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\hat{x}_{jk}\geq\sum_{j=1}^{\mathrm{J}}w_{j}\hat{\lambda}_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\hat{x}_{jk}$

Further, we have

(15)
$\sum_{j=1}^{\mathrm{J}}w_{j}\hat{\lambda}_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\frac{x_{jk}^{\star}}{1+\epsilon_{0}}\geq\frac{1-\epsilon_{0}}{1+\epsilon_{0}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}x_{jk}^{\star}$

Using Inequality (13) and Inequality (15) and plugging into inequality (14) we
get:

$\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}\tilde{x}_{jk}\geq\frac{(1-\epsilon_{0})^{2}}{(1+\epsilon_{0})^{2}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}\sum_{k\in\mathcal{K}_{j}:k_{\tau}=l}x_{jk}^{\star}$

Further note that $\frac{(1-\epsilon_{0})^{2}}{(1+\epsilon_{0})^{2}}\geq
1-4\epsilon_{0}$, yielding the final result. ∎

###### Proof of Theorem 2.

Recall that the time horizon of length $\mathrm{T}$ is divided into $\Phi+1$
phases, with $\epsilon 2^{\Phi}=1$. For simplicity we assume that $\epsilon$
is such that $\Phi=\log(1/\epsilon)$ is an integer. For each phase we ignore
any transmissions of the first $d_{max}$ slots (i.e., as if they were idle)
for the purpose of the analysis. The actual performance of the algorithm can
only be better than that. Then, any phase $\phi$ (with $\phi\geq 1$) is split
into the idle period $I_{\phi}$ of duration $d_{max}$ and an active period of
duration $T_{i}$. The total active time is
$\mathrm{T}^{\prime}=\sum_{\phi=1}^{\Phi}\mathrm{T}_{\phi}$ and the total time
is $\mathrm{T}=\mathrm{T}^{\prime}+d_{max}\Phi$. Recall that
$\mathrm{T}_{0}=\mathrm{T}_{1}=\epsilon\mathrm{T}^{\prime}$ and for the
remaining values $\mathrm{T}_{\phi}=\epsilon\mathrm{T}^{\prime}2^{\phi}$.

Let $\mu=\frac{a_{\max}^{2}}{2\lambda_{\min}^{2}}$. We saw in Lemma 5, that
using $\mu\frac{\log(2\mathrm{J}/\epsilon)}{\epsilon_{0}^{2}}$ samples, we
obtain a $1-4\epsilon_{0}$ approximation to $\mathrm{E}_{S}$ w.p. at least
$(1-\epsilon)$. Since we are scheduling with the corresponding scaled flow
probabilities, similarly to the proof of Theorem 1 in Section 4.1, we have a
loss of factor $3\epsilon$. Combining the various approximation ratio losses,
we obtain the following fraction of the optimal value for each period of
accuracy $\epsilon_{0}$:
$(1-4\epsilon_{0})(1-3\epsilon)(1-\epsilon)\geq(1-4\epsilon_{0})(1-4\epsilon)$

More specifically, in our case, in phase $\phi$ of total duration
$\mathrm{T}_{\phi}$ ($\phi\geq 1$), the algorithm can use all prior phases
$0,1,\cdots,\phi-1$ of total duration
$\sum_{i=0}^{\phi-1}\mathrm{T}_{i}=\mathrm{T}_{0}+\mathrm{T}_{0}(2^{\phi}-1)=\mathrm{T}_{0}2^{\phi}=\mathrm{T}_{\phi}$.
Therefore, during phase $\phi$, we use:
$\epsilon_{0}:=\epsilon_{\phi}=\sqrt{\mu\frac{\log(2\mathrm{J}/\epsilon)}{\mathrm{T}_{\phi}}}$,
and we obtain a fraction $(1-4\epsilon)(1-4\epsilon_{\phi})$ of the optimal.
Overall, we can identify the approximation ratio of Algorithm 2 by analyzing
the quantity.

(16)
$(1-4\epsilon)\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}}(1-4\epsilon_{\phi}),$

Indeed this can be seen as we are obtaining $(1-\epsilon_{\phi})(1-4\epsilon)$
approximation in phase $\phi$, for a fraction $\mathrm{T}_{\phi}/\mathrm{T}$
of the total time. We therefore have:

$\displaystyle\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}}(1-4\epsilon_{\phi})=\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}}-\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}}4\epsilon_{\phi}\geq\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}}-\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}^{\prime}}4\epsilon_{\phi}$
$\displaystyle\overset{(a)}{\geq}(1-\epsilon)\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}^{\prime}}-\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}^{\prime}}4\epsilon_{\phi}\overset{(b)}{=}(1-\epsilon)^{2}-\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}^{\prime}}4\epsilon_{\phi}\geq(1-2\epsilon)-\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}^{\prime}}4\epsilon_{\phi}.$

To see $(a)$, note that by assumption: $\mathrm{T}\geq
d_{max}(1+\frac{1}{\epsilon})\Phi$. We then have:

$\mathrm{T}^{\prime}=\mathrm{T}-d_{max}\Phi\geq\mathrm{T}-\frac{\mathrm{T}}{1+1/\epsilon}=\frac{\mathrm{T}}{1+\epsilon}\geq\mathrm{T}(1-\epsilon).$

from which $(a)$ follows by the reciprocal inequality after dividing with
$(1-\epsilon)$. Further, $(b)$ holds due to
$\sum_{\phi=1}^{\Phi}{\mathrm{T}}_{\phi}/{\mathrm{T}^{\prime}}=\frac{\mathrm{T}^{\prime}-\mathrm{T}_{0}}{\mathrm{T}^{\prime}}=\frac{\mathrm{T}^{\prime}-\epsilon\mathrm{T}^{\prime}}{\mathrm{T}^{\prime}}=1-\epsilon$.
We will proceed by showing
$\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}}\epsilon_{\phi}\leq\epsilon$
which will give us a bound:

(17)
$\sum_{\phi=1}^{\Phi}\frac{\mathrm{T}_{\phi}}{\mathrm{T}}(1-4\epsilon_{\phi})\geq
1-6\epsilon.$

Indeed, we have

$\sum_{\phi=1}^{\Phi}\mathrm{T}_{\phi}\epsilon/\mathrm{T}^{\prime}\overset{(a)}{=}\sqrt{\mu\log(2\mathrm{J}/\epsilon)\epsilon/\mathrm{T}^{\prime}}\sum_{\phi=1}^{\Phi}\sqrt{2}^{(\phi-1)}\overset{(b)}{\leq}\frac{\sqrt{\mu\log(2\mathrm{J}/\epsilon)\epsilon/\mathrm{T}^{\prime}}}{\sqrt{\epsilon}(\sqrt{2}-1)}=\frac{\sqrt{\mu\log(2\mathrm{J}/\epsilon)/\mathrm{T}^{\prime}}}{\sqrt{2}-1}\overset{(c)}{\leq}\epsilon.$

$(a)$ is from definition of $\epsilon_{\phi}$ and $\mathrm{T}_{\phi}$. In
$(b)$, we used the fact that
$\sum_{\phi=1}^{\Phi}\sqrt{2}^{\phi-1}=\frac{\sqrt{2}^{\Phi}-1}{\sqrt{2}-1}=(1/\sqrt{\epsilon}-1)/(\sqrt{2}-1)\leq\frac{1}{\sqrt{\epsilon}(\sqrt{2}-1)}$
since $2^{\Phi/2}\sqrt{\epsilon}=1$. Finally, $(c)$ holds under the assumption
that
$\mathrm{T}>\frac{1}{(\sqrt{2}-1)^{2}}\mu\log(2\mathrm{J}/\epsilon)/\epsilon^{2}+d_{max}\log(1/\epsilon)$
which implies
$\mathrm{T}^{\prime}>\frac{1}{(\sqrt{2}-1)^{2}}\mu\log(2\mathrm{J}/\epsilon)/\epsilon^{2}$.
Therefore, from (17) and due (16), we obtain a fraction
$(1-6\epsilon)(1-4\epsilon)\geq(1-10\epsilon)$ of $\mathrm{RI}(\mathrm{T})$. ∎

##  Appendix B General Distributions with stationary arrival rates

In this section we consider general distributions on the arrival processes
$\\{a_{j}^{t}\\}$ with stationary arrival rates, i.e.,
$\lambda_{j}^{t}\equiv\lambda_{j}$. First, we provide a preliminary definition
which allows us to extend Theorem 1.

###### Definition 0.

Consider a nonnegative integer-valued random variable $Z$ with distribution
$\mathcal{D}$, and suppose $Z$ is bounded by a constant $M$. Consider a
decomposition of $Z$ as the sum of $N$ independent random variables
$\\{Z_{i}:Z_{i}\leq M_{i}\\}$:

$Z=\sum_{i=1}^{N}Z_{i}.$

We define the dependency-degree of the decomposition as $\max_{i}M_{i}$. Then,
define the dependency-degree of distribution $\mathcal{D}$, as the minimum
dependency-degree of all decompositions for $\mathcal{D}$.

In particular, the Binomial and Bernoulli distributions, both have dependency-
degree $1$ (i.e., $M_{i}=1$). A scaled Bernoulli distribution which is either
$0$ or $M$, has the maximum possible dependency-degree of $M$. Further, it is
easy to see that the dependency-degree of the sum of two independent random
variables is equal to the maximum dependency degree of the two random
variables.

We now state Theorem 2, which generalizes Theorem 1.

###### Theorem 2.

Consider stationary arrival processes $\\{a_{j}^{t}\\}$, with the property
that each packet type’s total traffic within a fixed time window,
$\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}a_{j}^{t}$, has at most
dependency degree $D$ (for all $j$ and $t$). Then, given $\epsilon\in(0,1/3)$,
$\mathrm{FBPF}$ provides $(1-3\epsilon)$-approximation to
$\mathrm{RI}(\mathrm{T})$ when $\mathrm{C}_{\min}\geq
2D\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log(L/\epsilon)$ and
$\mathrm{T}\geq\frac{2d_{max}^{2}}{\epsilon}$.

The minimum link-capacity threshold for a given approximation ratio, now
depends on the dependency degree of the distribution of the total arrivals for
the packet type in a window of $(d_{max}+1)$ consecutive time slots. In other
words, a lower dependency-degree results in better performance guarantees,
given $\mathrm{C}_{\min}$.

###### Remark 9.

Due to Theorem 2, we can obtain guarantees for distributions with dependencies
across time slots. For example, if $a_{j}^{t}\in\\{0,1\\}$, and we allow the
arrivals of a packet type at a given time, to depend on the arrivals of that
type at different time slots, then there can be a maximum of $d_{max}+1$
packets depending on each other, which results in dependency degree of
$D=(d_{max}+1)$. An example of such a distribution is given in (Lee, 2021).

For the proof of Theorem 2, the analysis remains mostly identical to the case
of Bernoulli or Binomial arrivals. However we need to generalize Lemma 3 as
follows.

###### Lemma 0.

In the case that $\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}a_{j}^{t}$ has
at most dependency degree $D$, for all $j$ and $t$, using the forwarding
probabilities $\\{f_{j\ell}^{\tau\star}\\}$ for scheduling packets (Lines 1-1
of Algorithm 1), the probability of a packet being dropped is at most
$\epsilon$, if $\mathrm{C}_{\min}\geq
2D\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log\frac{L}{\epsilon}$.

###### Proof.

The proof in Lemma 3 leverages the independence of packet types. More
specifically the proof for Lemma 3 relied on Lemma 4 to show that any link’s
capacity is unlikely to be exceeded. This was then used to obtain through a
union bound the final result. We follow similar reasoning here, however when
using Lemma 4 we need to define a set of independent variables, which now
cannot be the arrival or not of a packet on the link.

Suppose all the packets that can arrive on a link $\ell$ are partitioned into
$N$ groups. We define the total consumption due to group $i$ as $Z_{g_{i}}$.
We then have $Z_{g_{i}}\leq D$ by assumption (although the assumption is
imposed on the arrivals for each packet type, this is propagated to the
arrivals at a link, due to the properties of the dependency-degree). We
proceed similarly to Lemma 3 by analyzing $\mathbb{E}{Z_{g_{i}}^{2}}$ and
$\mathbb{E}{Z_{g_{i}}}$. However, now, in contrast to Lemma 3, it is not true
that $\mathbb{E}{Z_{g_{i}}^{2}}=\mathbb{E}{Z_{g_{i}}}$. Indeed let us consider
a group $g$ with a maximum number of $I$ packets that can arrive in the
considered link, and let these arrivals be indicated through variables
$Y_{1},Y_{2},\cdots,Y_{I}$. Then

$\mathbb{E}{[Z_{g}^{2}]}=\mathbb{E}{\big{(}\sum_{i=1}^{I}Y_{i}\big{)}^{2}}=\mathbb{E}{\big{(}\sum_{i=1}^{I}Y_{i}+\sum_{i\neq
j}Y_{i}Y_{j}\big{)}}=\mathbb{E}{\big{[}\sum_{i=1}^{I}Y_{i}(1+\sum_{j:j\neq
i}^{I}Y_{j})\big{]}}\leq
D\mathbb{E}{\big{[}\sum_{i=1}^{I}Y_{i}}]=D\mathbb{E}{[Z_{g}]}$

Using the above, we prove (based on the notation of Lemma 4, similarly to
Lemma 3):

$S=\sum_{i=1}^{N}\mathbb{E}{[Z_{g_{i}}^{2}]}\leq D\mathbb{E}{Z}\leq
D\mathrm{C}_{\ell}/(1+\epsilon)=DS_{u},$

where $S_{u}$ and $\lambda$ are defined as in the proof of Lemma 3. Then the
exponent in Lemma 4, for $B=D$, becomes

$-\frac{\lambda^{2}}{2S+2D\lambda/3}\leq-\frac{1}{D}\frac{\epsilon^{2}S_{u}^{2}}{2S_{u}+2\epsilon
S_{u}/3}\leq-\frac{1}{D}\frac{\epsilon^{2}\mathrm{C}_{\ell}}{2(1+\epsilon)^{2}}$

from which we can see that $\mathrm{C}_{\ell}$ needs to scaled this time by a
factor of $D$ in order to replicate the results of the proof in Lemma 3. ∎

##  Appendix C Non-stationary distributions

###  C.1. Extension to non-stationary arrival rates

In this section we present details for the non-stationary extension of
Algorithm 1. First, we need to extend $\mathrm{F}_{S}$, to obtain
$\mathrm{F}_{NS}$. In essence, $\mathrm{F}_{NS}$ is obtained directly from
$\mathrm{EI}(\mathrm{T})$, by omitting Step 2 in Section 4.1. As a result,
randomizing according to the probabilities from the solution of
$\mathrm{F}_{NS}$ results in improved approximation ratios over using the
probabilities obtained from $\mathrm{F}_{S}$, as we do not need to leverage
Lemma 2 anymore (which adds a multiplicative factor of $(1-2d_{max}^{2}/T)$ in
the approximation ratio and thus worsens the quality of the guarantee for
small horizons). Alternatively, $\mathrm{F}_{S}$ could be interpreted as an
approximation of $\mathrm{F}_{NS}$.

The analysis follows similarly to the stationary case, by omitting the
intermediate step of Lemma 2. An analogous statement as in Lemma 3 can then be
derived for $\mathrm{F}_{NS}$. The LP $\mathrm{F}_{NS}$ is given below:

(18a) $\displaystyle\max_{\mathbf{f}}\quad$
$\displaystyle\sum_{t=1}^{\mathrm{T}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}^{t}\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0t}\quad(:=\mathrm{F}_{NS})$
(18b) s.t. $\displaystyle\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0t}\leq
1,\quad\forall j\in[\mathrm{J}],\forall t\in\mathrm{T},$ (18c)
$\displaystyle\sum_{\ell\in\mathtt{Inc}(v)}f_{j\ell}^{(\tau-1)t}=\sum_{\ell\in\mathtt{Out}(v)}f_{j\ell}^{\tau
t},\quad\forall v\in\mathcal{V},\tau\in[d_{max}],\forall j\in[\mathrm{J}],$
(18d) $\displaystyle f_{j\ell}^{0}=0,\ \
\forall\ell\not\in\mathtt{Out}(s_{j}),\quad f_{j\ell}^{d_{j}}=0,\ \
\forall\ell\not\in\mathtt{Inc}(d_{j}),$ (18e)
$\displaystyle\sum_{\tau=0}^{d_{max}}\sum_{j=1}^{\mathrm{J}}\lambda_{j}^{t-\tau}f_{j\ell}^{\tau(t-\tau)}\leq\frac{\mathrm{C}_{\ell}}{1+\epsilon},\quad\forall\ell\in\mathcal{L},\forall
t\in[\overline{\mathrm{T}}],$ (18f) $\displaystyle f_{j\ell}^{\tau t}\geq
0,\quad\forall
j\in[\mathrm{J}],\forall\ell\in\overline{\mathcal{L}},\forall\tau\in[d_{max}],t\in[\mathrm{T}].\quad$

In the optimization above, we conventionally consider $f_{j\ell}^{\tau t}=0$
for $t\leq 0$ and $t>\mathrm{T}$ (in order to simplify indexing in Equation
18e). We now state a theorem for the general non-stationary case.

###### Theorem 1.

Consider general arrival processes $\\{a_{j}^{t}\\}$, with the property that
each packet type’s total traffic within a fixed time window,
$\sum_{\tau=(t-d_{max})^{+}}^{t\land\mathrm{T}}a_{j}^{t}$, has at most
dependency degree $D$ (for all $j$ and $t$). Then, given $\epsilon\in(0,1/2)$,
and scheduling using the forwarding variables from $\mathrm{F}_{NS}$, yields a
$(1-2\epsilon)$-approximation to $\mathrm{RI}(\mathrm{T})$ when
$\mathrm{C}_{\min}\geq
2D\left(\frac{1+\epsilon}{\epsilon}\right)^{2}\log(L/\epsilon)$.

Note that Theorem 1 does not impose any constraints on the time horizon. The
use of $\mathrm{F}_{NS}$ for the forwarding variables comes with a cost in
terms of computational complexity. In particular, the number of variables and
constraints required to solve $\mathrm{F}_{NS}$ now scale with the horizon
$\mathrm{T}$ (which was the reason we considered $\mathrm{F}_{S}$ in
stationary environments). This can be mitigated as we discuss next in Section
C.2.

###  C.2. Frame-Based solutions for non-stationary arrival rates

As discussed in Section C.1, in the presence of non-stationary traffic, the
number of variables and constraints for solving $\mathrm{F}_{NS}$ scales with
horizon $\mathrm{T}$. However, a frame-based solution can be adopted to
mitigate the computational cost of solving $\mathrm{F}_{NS}$, while offering
additional benefits which we discuss below.

In a frame-based construction, we aim to partition horizon $[\mathrm{T}]$ into
consecutive time slots, named “frames“, of some maximum duration $H$, with
$H\ll\mathrm{T}$. Then, we solve $\mathrm{F}_{NS}$ for each frame
independently, and in order to forward packets that have arrived within a
frame, we use the forwarding variables obtained from the solution of
$\mathrm{F}_{NS}$ at the beginning of that frame.

Frame-based constructions offer three motivating benefits: (a) We need not
know (or estimate) ahead of time the arrival rates for the entire horizon, but
only for the next frame. This is particularly relevant to prediction-based
methods, since the arrival rates are often predictable in the short-term
future, but difficult to predict over longer horizons. (b) The complexity of
solving the LP is reduced. For example, if we are using an LP solver with
complexity $O(n^{k})$ for some $k>1$, where $n$ the number of variables, then,
as $n\propto\mathrm{T}$, we benefit by solving the problem in frames of size
$H$, which results in $\mathrm{T}/H$ optimization problems, with $n\propto H$
each, and therefore a reduction of the complexity from $O(\mathrm{T}^{k})$ to
$O(\frac{\mathrm{T}}{H}H^{k})=O(TH^{k-1})$. (c) The computational complexity
of finding the forwarding variables, is more evenly split across the time
horizon, rather than requiring the calculation of the forwarding variables for
the entire horizon at the beginning of the execution of the algorithm.

A frame-based method however might lead to a degradation in the approximation
ratio obtained, unless it is meticulously designed, or certain assumptions are
imposed on the traffic. This is because the scheduling decisions within a
frame do not take into consideration the possible arrivals in other frames,
despite the fact that these decisions may impact these future frames. It then
becomes important to characterize the losses due to a frame construction.
Fortunately, these losses can be localized around the boundaries of the frame
(the end or the beginning of the frame), since, due to the bound $d_{max}$ on
the deadlines, the impact of decisions can be shown to be limited to at most
$d_{max}$ slots away from the boundaries of each frame. However, without
additional assumptions, due to the arbitrary non-stationary distributions on
the arrivals, a fixed-frame construction could have significant losses, if
most of the high-weight traffic is concentrated on these boundaries. This can
be resolved by either placing certain mild assumptions on the non-stationarity
of the traffic, or by constructing a more involved frame construction, with
frames of potentially different frame lengths. We illustrate a variable-sized
frame construction.

Lemma 2 characterizes the losses of using a variable-sized frame construction.

###### Lemma 0.

$\mathrm{F}_{NS}$ admits a frame-based solution, with frames of variable
length, of at most $H=(2+4/\epsilon)d_{max}$ time slots each, with

$W_{\mathrm{F}_{NS}}^{\mathrm{FRAME}}\geq(1-\epsilon)W_{\mathrm{F}_{NS}}$

###### Proof.

We assume for simplicity that $\mathrm{T}$ is a multiple of $d_{max}N$ for
some $N$. Then, let us partition $\mathrm{T}$ into groups of $d_{max}$ slots,
referred to as mini-frames. Further, define cycles
$\mathcal{C}_{1},\mathcal{C}_{2},\cdots$, with each cycle consisting of $N$
consecutive mini-frames. Therefore each cycle has duration $d_{max}N$. Assume
we are processing some arbitrary cycle, say $\mathcal{C}_{1}$ without loss of
generality. Solve $\mathrm{F}_{NS}$ on each mini-frame for $\mathcal{C}_{1}$
and $\mathcal{C}_{2}$. Suppose the optimal solutions within each mini-frame
are of reward $W[1],W[2],\cdots,W[N]$ for the first cycle, and
$W^{\prime}[1],\cdots,W^{\prime}[N]$ for the second cycle. Find $i$ with
minimum reward for the first cycle $W_{min}=W[i]$, and similarly $i^{\prime}$
for the second cycle. These minimum reward frames can be assigned to be the
boundaries of the frame and intuitively they should result in limited losses.
Indeed we define the frame to be the set of consecutive mini-frames between
(but not including) mini-frames $i$ and $i^{\prime}$. We argue that the loss
in each cycle due to omitting these mini-frames will be limited. Consider
mini-frame $i$, which is part of cycle $\mathcal{C}_{i}$. Since our frame-
based solution will not schedule packets in this mini-frame, it will suffer a
loss from arrivals in that mini-frame. However outside the mini-frame, the
remaining reward will be at least as large as the non-frame-based solution. As
a consequence of the definition of $i$, we can further argue that these losses
are limited. In particular, the loss cannot be higher than a fraction
$2/(N-1)$ of the reward in the cycle. Indeed, if that was the case, $i$
wouldn’t have been the minimum reward mini-frame in the cycle. Ignoring all
transmissions in this mini-frame can result in a loss of reward of the optimal
policy due to transmissions in this mini-frame, which cannot be higher than a
fraction $2/(N-1)$ of the reward in the entire frame. If $2/(N-1)=\epsilon$,
that is, for $N=1+2/\epsilon$, we have an $\epsilon$ loss in each cycle. The
above scheme, only requires knowing or estimating the distribution of arrivals
over the incoming two cycles. Figure 9 illustrates the construction through an
example. ∎

Figure 9. A frame-based construction, where frames are defined as a sequence
of mini-frames between two mini-frames of minimum reward. The figure indicates
an example with $N=5$ mini-frames in each cycle. In the first cycle, mini-
frame $3$ has minimum reward $W[i]$. In the second cycle, mini-frame $4$ has
minimum reward $W^{\prime}[i^{\prime}]$.

###  C.3. Reducing Complexity for Periodic Arrival Rates

(a) Range of capacity constraints for $\mathrm{EI}(\mathrm{T})$

(b) Modified range for $\mathrm{\overline{EI}}(\mathrm{T})$

Figure 10. We consider the same example as in Figure 4, with $d_{max}=2$, but
now, rather than a stationary traffic, we assume periodic traffic with
$\mathrm{T}=2\widetilde{\mathrm{T}}$, illustrated in Figure 10(a). The broad
symmetry of the problem breaks again at $1,2,\mathrm{T}-1,\mathrm{T}$. After
the modification shown in Figure 10(b), the first $\widetilde{\mathrm{T}}$
slots are indistinguishable from the final $\widetilde{\mathrm{T}}$ slots in
$\mathrm{\overline{EI}}(\mathrm{T})$. For example the arrivals rates impacting
slot $1$ are identical to the arrival rates impacting slot
$1+\widetilde{\mathrm{T}}$. This implies a symmetry in shifts by
$\widetilde{\mathrm{T}}$ slots, which yields equivalent solutions following by
shifting a solution in the time index by $\widetilde{\mathrm{T}}$ (similarly
to the proof of Lemma 2).

Here, we provide details on how a periodic set of forwarding variables can be
used to obtain a near-optimal solution for $\mathrm{RI}(\mathrm{T})$. For
simplicity, it is assumed that $\widetilde{}\mathrm{T}$ divides the horizon
$\mathrm{T}$. Then, to argue that a periodic solution to
$\mathrm{EI}(\mathrm{T})$ is near optimal, we need to provide an argument
similar to that in the proof of Lemma 2 in a stationary setting (i.e.,
generalize for $\widetilde{\mathrm{T}}\geq 1$).

Indeed, it is not hard to verify that the modified program,
$\mathrm{\overline{EI}}(\mathrm{T})$ is now symmetric with respect to shifts
in the index of the variables by $\widetilde{\mathrm{T}}$. Refer to Figure 10
for an illustration of the point. Then, any general optimal solution to
$\mathrm{EI}(\mathrm{T})$, can be made periodic, by noticing that a shift of
all time-indices in the solution of $\mathrm{EI}(\mathrm{T})$ by
$\widetilde{\mathrm{T}}$ results in another optimal solution. Then, averaging
all shifts, yields a solution that only depends on the timeslot within the
period, and hence, is periodic. From that argument as in Lemma 2, we infer
that a periodic near-optimal solution exists for $\mathrm{EI}(\mathrm{T})$,
and hence a simplification similar to $\mathrm{E}_{S}$ can be obtained, which
converted in a flow form, yields the following program $\mathrm{F}_{NSP}$.

(19a) $\displaystyle\max_{\mathbf{f}}\quad$
$\displaystyle\sum_{t=1}^{\widetilde{\mathrm{T}}}\sum_{j=1}^{\mathrm{J}}w_{j}\lambda_{j}^{t}\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0t}\quad(:=\mathrm{F}_{\mathrm{NSP}})$
(19b) s.t. $\displaystyle\sum_{\ell\in\mathtt{Out}(s_{j})}f_{j\ell}^{0t}\leq
1,\quad\forall j\in[\mathrm{J}],\forall t\in[\widetilde{}\mathrm{T}],$ (19c)
$\displaystyle\sum_{\ell\in\mathtt{Inc}(v)}f_{j\ell}^{(\tau-1)t}=\sum_{\ell\in\mathtt{Out}(v)}f_{j\ell}^{\tau
t},\quad\forall v\in\mathcal{V},\tau\in[d_{max}],\forall j\in[\mathrm{J}],$
(19d) $\displaystyle f_{j\ell}^{0}=0,\ \
\forall\ell\not\in\mathtt{Out}(s_{j}),\quad f_{j\ell}^{d_{j}}=0,\ \
\forall\ell\not\in\mathtt{Inc}(d_{j}),$ (19e)
$\displaystyle\sum_{\tau=0}^{d_{max}}\sum_{j=1}^{\mathrm{J}}\lambda_{j}^{(t-\tau)\mathrm{mod}\widetilde{}\mathrm{T}}f_{jl}^{\tau((t-\tau)\mathrm{mod}\widetilde{}\mathrm{T})}\leq\frac{\mathrm{C}_{l}}{1+\epsilon},\quad\forall\ell\in\mathcal{L},t\in[\widetilde{}\mathrm{T}],$
(19f) $\displaystyle f_{j\ell}^{\tau t}\geq 0,\quad\forall
j\in[\mathrm{J}],\forall\ell\in\overline{\mathcal{L}},\forall\tau\in[d_{max}],t\in[\mathrm{T}].\quad$

##  Appendix D Relation between flows and route-schedules

Recall that in this work we randomize over route-schedules for each packet
(e.g. following our analysis in $\mathrm{E}_{S}$). In order to assign
probabilities to all route-schedules compactly, we reformulate
$\mathrm{E}_{S}$, inspired by ideas from the maximum flow and multicommodity
flow problem (Magnanti and Orlin, 1993). First, we need to define a new graph
over which we can refer to flows. To do this, we initially associate route-
schedules to regular routes in a time-expansion of the original network, such
as the one indicated in Figure 11. In the time-expansion each link is copied
$d_{max}+1$ times (with each copy representing the age at which the packet can
be scheduled over the link) and each node is copied $d_{max}+2$ (representing
the location of the packet at different ages).

Then, all route schedules $\mathcal{K}_{j}$ correspond to paths in the time-
expansion from node $(s_{j}:0)$ to $(z_{j}:d_{max}+1)$. However, from standard
results in the network flow literature (Magnanti and Orlin, 1993), we can
combine a composition of weighted paths in any graph, into a flow, and vice
versa, we can decompose a weighted flow into a set of paths (Goldberg and Rao,
1998). In particular, by making sure that the weights are appropriately scaled
(i.e., the total weight of the weighted paths is at most one), we can
interpret the weights as probabilities of scheduling over each path.

Applied to our problem, this allows us to work with either variables directly
on route-schedules $x_{jk}$ ($\mathrm{E}_{S}$), or their flow-inspired
corresponding forwarding variables $f_{j\ell}^{\tau}$. For instance, note that
the example of forwarding-variables defined in Figure 2(a) can also be viewed
as a flow over the time-expansion in Figure 11.

###  D.1. From forwarding variables to route-schedules

As in the network flow literature (Magnanti and Orlin, 1993; Raghavan and
Tompson, 1987), there are multiple ways to assign forwarding variables to
probabilities over route-schedules. In the paper we presented a hop-by-hop
approach, which uses directly the intuitive definition of a forwarding
variable. An alternate mapping however could be through an iterative process
that in each step finds a route-schedule of maximum possible probability,
subtracts the probability assigned to it from the flow, and repeating the
process. This results in exhausting for each packet type at least one edge of
the time-expansion per iteration. Since there is a maximum of
$(d_{max}+1)|\overline{\mathcal{L}}|$ edges in the time-expansion, this
results in a limited number of route-schedules with positive probability per
packet type.

Figure 11. A time-expansion of the network in Figure 2(a) with the packet type
considered in Figure 2(a). Note that there are only $3$ route-schedules for
that packet type, and these correspond to all paths from $u:0$ to $v:3$
(recall that the deadline was $2$). For example path $(u:0,u:1,a:2,v:3)$
corresponds to route schedule $k=[(u,u)\quad(u,a)\quad(a,v)]$. Each route-
schedule can be assigned a probability through the forwarding variables
indicated in both Figure 2(a) as well as the time-expansion. Note however that
in the time-expansion, the forwarding variables consist also of a valid flow
as in the max-flow problem (Magnanti and Orlin, 1993).

##  Appendix E ADDITIONAL SIMULATIONS

###  E.1. Synthetic Datasets

Grid Network $4\times 4$. We simulated the performance of our algorithm in a
grid network with $16$ nodes, in a $4\times 4$ arrangement. In the grid
topology, a node which is not at an edge or corner of the grid, is connected
with $4$ nodes (top,right,bottom,left nodes). Corner nodes only have 2
connections instead. For the simulation parameters, we used $\mathrm{J}=30$,
and randomly assigned sources and destinations, as earlier, and weights
randomly chosen in $\\{10,\cdots,100\\}$. The results are shown in Figure 12
for different $d_{max}$ and $\lambda_{j}$ ranges (i.e. each $\lambda_{j}$ is