\section{Introduction} \label{Sec:Intro}

The rapid growth of information technology has resulted in increasingly complex network systems and poses challenges in obtaining explicit knowledge of system dynamics. For instance, due to security or economic concerns, a number of network systems are built as overlay networks, e.g. caching overlays, routing overlays and security overlays \cite{sitaraman2014overlay}. In these cases, only the overlay part is fully controllable by the network administrator, while the underlay part remains uncontrollable and/or unobservable. The ``black box" components make network control policy design challenging.


In addition to the challenges brought by unknown system dynamics, many of the current network control algorithms (e.g., MaxWeight~\cite{tassiulas1990stability} and Drift-plus-Penalty~\cite{neely2008fairness}) aim at stabilizing the system. However, existing optimization methods for long-term performance metrics (e.g. queue backlog and delay) only work for specific structures \cite{buyukkoc1985cmu, maguluri2016heavy}, while the development for general setting remains insufficient. 


To overcome the above challenges, it is desirable to apply inference and learning schemes. A natural approach is reinforcement learning, which learns and reinforces a good decision policy by interacting with the environment and receiving feedbacks.
Reinforcement learning methods provide a framework for the design of learning policies for general networks. There have been two main lines of work on reinforcement learning methods: model-free reinforcement learning (e.g., Q-learning \cite{watkins1992q}, policy gradient \cite{sutton2000policy}) and model-based reinforcement learning (e.g., UCRL \cite{jaksch2010near}, PSRL \cite{osband2013more}). In this work, we focus on the model-based approach. 



\subsection{Related Work} \label{Sec:RelatedWork}

Existing methods for reducing the average job delay (or equivalently reducing the average queue backlog) of general networks can be classified into three types: equivalent constraint, Lyapunov drift and Markov decision process (MDP) \cite{cui2012survey}. However, most equivalent constraint approaches can only be applied to single-hop networks \cite{cui2012survey}. Therefore, we focus on discussing the latter two classes of algorithms, as they are more universal and can be applied to multi-hop network systems.

We first discuss the Lyapunov drift approach, which generally does not require learning the dynamics for network problems with hidden dynamics. The Lyapunov drift approach has been widely applied in numerous network control problems. For instance, in dynamic server allocation problems, a Lyapunov-drift-based algorithm named the longest connected queue (LCQ) can stabilize the queue backlog without knowing the underlying dynamics (e.g., arrival rates, channel statistics) \cite{tassiulas1993dynamic,ganti2007optimal}. For multiclass routing networks, the Backpressure algorithm was designed under the Lyapunov drift framework and only requires the observation of queue backlogs \cite{tassiulas1990stability,awerbuch1993simple,neely2003dynamic}. Similarly, Lyapunov drift methods that do not require learning the network dynamics have been proposed for switch scheduling \cite{mekkittikul1998practical,keslassy2001analysis,shah2012switched} and inventory control problems \cite{neely2010dynamic, dai2005maximum} etc. Most Lyapunov drift algorithms are shown to be throughput optimal \cite{tassiulas1990stability}, i.e., they can stabilize the system whenever the system is stabilizable. Also, the Lyapunov drift approach usually only requires solving a linear programming problem and does not suffer from the ``curse of dimensionality'', which makes it applicable to large-scale queueing systems. However, without learning the system dynamics, the Lyapunov drift approach generally cannot guarantee minimum queueing delay. For instance, for dynamic server allocation problems, the optimal policy has only been developed for highly symmetric systems (i.e., uniform external arrival rates, uniform connectivity probabilities and uniform success rates) \cite{ganti2007optimal}. Another well-known example is switching system, for which the Maximum Matching policy has been shown to be close to the optimal policy when the system is in the heavy traffic regime~\cite{maguluri2016heavy}. In Section 5.1, 5.2 and 5.4, we conduct numerical experiments and show that our approach significantly outperforms the Lyapunov drift methods.

The MDP approach models the queueing system control problem as an MDP that aims to minimize the long-term average queue backlog. Classical algorithms to solve MDP problems include value iteration and policy iteration \cite{bertsekas2017dynamic}. Note that although the MDP approach can minimize the average job delay, it can only be applied to networks with finite buffer size, which is unrealistic for many practical network models. Therefore, the MDP approach fails to minimize the average job delay of general stochastic networks and is typically applied to MDPs of small scale \cite{yeh2001multiaccess,yeh2003throughput} or special structures \cite{guestrin2003efficient, qu2020scalable}. However, traditional MDP approach requires explicit knowledge of the network dynamics. Since the network problems studied in this work have hidden dynamics, it is necessary to learn the dynamics when applying the MDP approach. A related field is reinforcement learning. The majority of reinforcement learning algorithms apply approximations to reduce the computation complexity, including state representation \cite{de2018integrating,maillard2013selecting}, value function approximation \cite{baird1995residual,jin2020provably} and pruning \cite{livne2020pops}. Reinforcement learning methods have been applied to a wide range of network control problems, like routing \cite{boyan1994packet,peshkin2002reinforcement,chen2018deep,yu2018drom}, spectrum access \cite{naparstek2017deep,wang2018deep,guan2020robust}, switch scheduling \cite{brown2001switch}, state probing \cite{brun2016big} etc. Although these reinforcement learning methods achieve satisfactory performance in simulation, they lack performance guarantees.

We consider model-based reinforcement learning methods, which tend to be more analytically tractable. Conventional model-based reinforcement learning methods like UCRL \cite{jaksch2010near} and PSRL \cite{osband2013more} only work for finite-state-space systems, yet queueing systems are usually modeled to have unbounded buffer sizes. The work in \cite{liang2019optimal} assumes that the system has an admission control scheme to keep queue backlogs finite. The resulting system can be model as a MDP with finite states, where UCRL can be applied. In \cite{theocharous2017posterior}, the authors modify the PSRL algorithm to deal with MDPs with large state space, yet the algorithm requires the MDP to have a finite bias span, which is unrealistic for problems that aim to minimize the average cost with a countably infinite state space.

In summary, existing methods either cannot guarantee optimal delay (the Lyapunov drift approach), or lack performance guarantees (heuristic reinforcement learning) or can only deal with finite state spaces (the MDP approach and theoretical reinforcement learning). In this paper, we aim to develop an algorithm that achieves provable optimality for networks with countably infinite state spaces. To the best of our knowledge, our algorithm is the first to achieve minimum average queue backlog in general networks with unbounded buffers. Our analysis also offers a possible roadmap to solving general MDPs with countably infinite state spaces.




\subsection{Our Contributions} \label{Sec:Contribution}

We apply model-based reinforcement learning to queueing networks with unbounded state spaces and unknown dynamics. Our approach leverages the fact that for a vast class of stable queueing systems, the probability of the queue backlog being large is relatively small. This observation motivates us to focus on learning control policies over a finite subset of states that the system visits with high probability. Our main contributions are summarized as follows. 

\begin{itemize}
    \item We propose a model-based RL algorithm that can deal with unbounded state spaces. In particular, we introduce an auxiliary system with the state space bounded by a threshold $U$. Our approach employs a piecewise policy: for states below the threshold, a model-based RL algorithm is used; for all other states, a simple baseline algorithm is applied.
    
    \item We establish that the episodic average queue backlog under the proposed algorithm can be arbitrarily close to the optimum with a large threshold $U$. In particular, by applying Lyapunov analysis, we characterize the gap to the optimal performance as a function of the threshold $U$. In addition, our proof technique may be of independent interest for analyzing the convergence of other reinforcement learning algorithms for queueing networks.
    
    \item Simulation results on dynamic server allocation, routing, and network switching problems corroborate the validity of our theoretical guarantees. In particular, the proposed algorithm effectively achieves a small average queue backlog, with the gap to optimum diminishing as the threshold $U$ increases.
\end{itemize}

The paper is organized as follows. In Section \ref{Sec:Prob} we formulate the problem and introduce notations. Section \ref{Sec:Approach} gives an outline of our approach, presents required assumptions and the proposed algorithm. We conduct theoretical performance analysis regarding convergence and optimality of the proposed algorithm in Section \ref{Sec:Theory}. We evaluate the algorithm in various settings and the numerical results are given in Section \ref{Sec:Sim}.






\section{System Model} \label{Sec:Prob}

We consider a discrete time network with the general topology of a directed graph $\mathcal{G} = (\mathcal{N}, \mathcal{L})$, where $\mathcal{N}$ is the set of nodes and $\mathcal{L}$ is the set of links. Each node maintains one or more queues for undelivered packets, and each queue has an unbounded buffer. We denote by $D$ the number of queues. Nodes represent the locations where data packets are generated, relayed and/or processed in the network. Links represent the communication links.

During each time slot, external data packets arrive at nodes where they are either processed and then depart the system or are relayed. For relayed packets, the communication links can be stochastic, i.e. data transmissions between nodes can fail. We assume the underlying dynamics are partially or fully unknown. This model captures a large class of queueing networks that involve routing, scheduling and switching.

We suppose that the system has the Markovian property: the probability distributions of the queue backlogs in the next time slot only depend on the current queue backlogs and the control decision. In other words, the system can be modeled as a countable-state-space MDP $\mathcal{M}$ with average cost as follows. 

\begin{itemize}
\item State space $\mathcal{S}$: the set of queue backlog vectors, i.e., $\mathcal{S} \triangleq \mathbb{N}^D$.
\item Action space $\mathcal{A}$: the set of feasible control decisions, which are specified by the problem setting. In this paper we only consider finite action space.
\item State-transition probability $p \left( \bm{Q}' \mid \bm{Q}, a \right)$: the probability of transitioning into state $\bm{Q}'$ from state $\bm{Q}$ with action $a$. For simplicity, we assume that the magnitude that a queue backlog can change during one time slot is upper bounded by a constant $W$. Let $\mathcal{R} \left( \bm{Q} \right)$ denote the set of one-step reachable queue backlog vectors from state $\bm{Q}$.
\item Cost function $c \left( \bm{Q} \right)$: the total backlog of $D$ queues, i.e., $c \left( \bm{Q} \right) \triangleq \sum_{i=1}^{D} Q_i$.
\end{itemize}

As discussed in Section \ref{Sec:RelatedWork}, a large number of queueing networks have been studied extensively and various stabilizing policies that aim to keep the average total queue backlog finite have been developed. Here we take a step further to go beyond stability and our goal is to minimize the average queue backlog.

For readers' convenience, we summarize the notations used in this paper in Table \ref{Tab:Notations}.

\begin{table}[htbp]\caption{Notations}
\begin{center
\begin{tabular}{r p{12cm} }
\toprule
$D$ & The number of queues in the queueing network \\
$\bm{Q}$ & The $D$-dimensional queue backlog vector of the queueing network \\
$Q_i$ & The queue backlog of node $i$ \\
$Q_{\max}$ & The maximum entry of $\bm{Q}$ \\
$\mathcal{R} \left( \bm{Q} \right)$ & The set of one-step reachable queue backlog vectors from $\bm{Q}$ \\
$W$ & The magnitude that a queue backlog can change during one time slot \\
$U$ & The buffer size of the auxiliary system \\
$\mathcal{M}$ & The countable-state-space MDP of the original queueing network \\
$\tilde{\mathcal{M}}$ & The MDP of the auxiliary system \\
$\mathcal{S}$ & The state space of $\mathcal{M}$ \\
$\tilde{\mathcal{S}}$ & The state space of $\tilde{\mathcal{M}}$ \\
$\mathcal{A}$ & The action space of $\mathcal{M}$ \\
$p \left( \cdot \mid \cdot, \cdot \right)$ & The state-transition kernel of $\mathcal{M}$  \\
$\tilde{p} \left( \cdot \mid \cdot, \cdot \right)$ & The state-transition kernel of $\tilde{\mathcal{M}}$  \\
$\pi_0$ & The known stabilizing policy of the queueing network \\
$\pi_{\mathrm{rand}}$ & The random policy that selects an action in $\mathcal{A}$ uniformly \\
$\pi^*$ & The stationary policy that minimizes the expected average total queue backlog of $\mathcal{M}$ \\
$\tilde{\pi}^*$ & The stationary policy that minimizes the expected average total queue backlog of $\tilde{\mathcal{M}}$ \\
$\rho^*$ & The expected average queue backlog of $\mathcal{M}$ when applying $\pi^*$ \\
$\tilde{\rho}^*$ & The expected average queue backlog of $\tilde{\mathcal{M}}$ when applying $\tilde{\pi}^*$ \\
$\tilde{\Phi}^* (\cdot)$ & The Lyapunov function with negative drift regarding $\tilde{\pi}^*$ \\
$\beta$ & The order of $\tilde{\Phi}^* (\cdot)$ \\
$\mathcal{S}^{\mathrm{in}}$ & The set of $\bm{Q}$'s with $\tilde{\Phi}^* \left( \bm{Q} \right) \leqslant (U-W)^\beta$ \\
$\mathcal{S}^{\mathrm{out}}$ & $\mathcal{S} \setminus \mathcal{S}^{\mathrm{in}}$ \\
$p^{\pi + \pi'} (\cdot)$ & The stationary probability distribution in $\mathcal{M}$ with $\pi$ applied to $\mathcal{S}^{\mathrm{in}}$ and $\pi'$ applied to $\mathcal{S}^{\mathrm{out}}$ \\
$\mathbb{E}_{\pi} [\cdot]$ & The expection of a random variable in $\mathcal{M}$ with $\pi$ applied to $\mathcal{S}$ \\
$\mathbb{E}_{\pi + \pi'} [\cdot]$ & The expection of a random variable in $\mathcal{M}$ with $\pi$ applied to $\mathcal{S}^{\mathrm{in}}$ and $\pi'$ applied to $\mathcal{S}^{\mathrm{out}}$ \\
$\tilde{\mathbb{E}}_{\tilde{\pi}} [\cdot]$ & The expection of a random variable in $\tilde{\mathcal{M}}$ with $\tilde{\pi}$ applied to $\tilde{\mathcal{S}}$ \\

\bottomrule
\end{tabular}
\end{center}
\label{Tab:Notations}
\end{table}






\section{Our Approach} \label{Sec:Approach}

Classical model-based reinforcement learning fits a model of state transition kernel to observed data and then solves the dynamic programming problem on the estimated system. The challenge of applying such an approach to countable-state MDPs arises from both of estimating model parameters and solving the estimated MDP, due to the fact that the state space is unbounded. 


Here we introduce an auxiliary system $\tilde{\mathcal{M}}$ with a bounded state space. We only apply reinforcement learning techniques on the constructed bounded state space, while simply apply a \emph{known} stabilizing policy $\pi_0$ (cf. Assumption~\ref{Asp:KnownPolicy}) to the rest of the states. We show that the performance gap between the proposed algorithm and the optimal policy can be made arbitrarily small by designing appropriate $\tilde{\mathcal{M}}$.



\subsection{Overview} \label{Sec:Outline}

We first provide an overview of our approach and defer a detailed description of the algorithm to Section~\ref{Sec:Algorithm}. Our reinforcement learning method operates in an episodic manner. 

We apply a decaying $\epsilon$-greedy method to decide whether an episode should conduct exploration or exploitation. For each episode, we perform exploration (i.e., apply some random policy to collect diverse samples) with probability $\epsilon_i$, and we conduct exploitation (i.e., using current estimates of dynamics to compute an estimated optimal policy) with probability $1 - \epsilon_i$. At the beginning of the learning process, we tend to explore the system to collect samples, and thus $\epsilon_i$ is relatively large when $i$ is small. As the learning process goes on, we gradually obtain enough samples and are close to the optimal policy, and thus $\epsilon_i$ is decreased to exploit the learned policy. The scheme has been applied in network control \cite{alaya2008dynamic,sodagari2011anti,guan2020robust} to achieve a trade-off between exploration and exploitation.

For an exploration episode, we apply a randomized policy $\pi_{\mathrm{rand}}$ that takes action uniformly to obtain samples for the estimation of state-transition probabilities in $\mathcal{M}$. However, since $\mathcal{M}$ has a countably infinite state space, we instead estimate the model for an auxiliary system $\tilde{\mathcal{M}}$ with a bounded state space. The auxiliary system $\tilde{\mathcal{M}}$ has a threshold $U$: the system has the same dynamics as the real one, with the only difference that each queue has a bounded buffer size $U$. For each queue in $\tilde{\mathcal{M}}$, when its queue backlog reaches $U$, new packets to it will be dropped. Mathematically, the state space of $\tilde{\mathcal{M}}$ is given by $\tilde{\mathcal{S}} \triangleq \{ \bm{Q} \in \mathcal{S} : Q_{\max} \leqslant U \}$ where $Q_{\max} \triangleq \max Q_i$. The auxiliary system $\tilde{\mathcal{M}}$ shares the same action space $\mathcal{A}$ and cost function $c \left( \bm{Q} \right)$ as $\mathcal{M}$. With the introduction of $\tilde{\mathcal{M}}$, an exploration episode operates as in Figure \ref{Fig:Approach_1}: we only apply $\pi_{\mathrm{rand}}$ to states in $\tilde{\mathcal{S}}$, while simply applying the \emph{known} stabilizing policy $\pi_0$ to other states.

For an exploitation episode, we first repartition the state space $\mathcal{S}$ due to technical reasons. We denote by $\tilde{\pi}^*$ the optimal policy for $\tilde{\mathcal{M}}$, and define a Lyanopuv function $\tilde{\Phi}^* (\cdot)$ (cf. Assumption~\ref{Asp:DriftOptimalPolicy}). Denoting by $\beta$ the order of $\tilde{\Phi}^* (\cdot)$, we partition $\mathcal{S}$ in the following manner:
$$
\begin{cases}
\mathcal{S}^{\mathrm{in}} \triangleq \left\{ \bm{Q} : \tilde{\Phi}^* \left( \bm{Q} \right) \leqslant (U-W)^\beta \right\} \\
\mathcal{S}^{\mathrm{out}} \triangleq \mathcal{S} \setminus \mathcal{S}^{\mathrm{in}}
\end{cases}
.
$$
As illustrated in Figure \ref{Fig:Approach_2}, during an exploitation episode, we first compute an estimated $\tilde{\pi}^*$ using the estimated dynamics obtained from exploration episodes. We then apply the estimated $\tilde{\pi}^*$ to states in $\mathcal{S}^{\mathrm{in}}$ and $\pi_0$ to states in $\mathcal{S}^{\mathrm{out}}$ throughout the episode.

\begin{figure}
  \centering
  \begin{subfigure}[b]{0.47\textwidth}
      \centering
      \includegraphics[width=\linewidth]{Approach_1}
      \caption{Exploration}
      \label{Fig:Approach_1}
  \end{subfigure}%
  \hfill
  \begin{subfigure}[b]{0.47\textwidth}
      \centering
      \includegraphics[width=\linewidth]{Approach_2}
      \caption{Exploitation}
      \label{Fig:Approach_2}
  \end{subfigure}
  \caption{Schemetic illustration of our approach (when $D=2$).}
  \label{Fig:Approach}
\end{figure}


We will show that the average queue backlog under our algorithm  converges to the optimal average queue backlog $\rho^*$ as we increase $U$. We divide the analysis into two stages: before and after $\tilde{\pi}^*$ is learned. For the first stage, our model-based reinforcement learning approach applies $\epsilon$-greedy exploration. We show that the proposed algorithm gradually obtains $\tilde{\pi}^*$ with high probability. For the second stage, by applying drift analysis on the Markov chain, we show that when $\tilde{\pi}^*$ is applied for states in $\mathcal{S}^{\mathrm{in}}$ and $\pi_0$ for states in $\mathcal{S}^{\mathrm{out}}$, the probability of queue backlog exceeding into $\mathcal{S}^{\mathrm{out}}$ decays exponentially with $U$. In addition, whenever $\bm{Q}$ leaves $\mathcal{S}^{\mathrm{in}}$, the expected accumulated queue backlog before $\bm{Q}$ returns back to $\mathcal{S}^{\mathrm{in}}$ can be upper bounded as a polynomial term in $U$. Together, we show that the gap between our result and the optimal average queue backlog $\rho^*$ is upper bounded by $\mathcal{O} \left( \mathrm{poly}(U) / \exp \left( U \right) \right)$, which diminishes as $U$ increases.


\subsection{Assumptions}


Our algorithm can be applied to a broad class of network problems. To establish rigorous performance guarantee, we need to impose some assumptions on $\mathcal{M}$ and $\tilde{\mathcal{M}}$. We will discuss these assumptions and argue that they are reasonable under many queueing networks.

\medskip

\noindent \textbf{Assumption on $\pi_0$:} As introduced in Section~\ref{Sec:Outline}, we only learn the optimal policy for an auxiliary system with a bounded state space, while for the rest of the states, we apply a \emph{known} stabilizing policy $\pi_0$ to control the queue backlog. We  quantify the stability of $\pi_0$ as follows.
\begin{assumption} \label{Asp:KnownPolicy}
There exists a \emph{known} policy $\pi_0$, a Lyapunov function $\Phi_0: \mathcal{S} \rightarrow \mathbb{R}_{+}$ and constants $a, \alpha, \epsilon_0, B_0 > 0$ such that the following properties hold:
\begin{enumerate}
    \item $\Phi_0 \left( \bm{Q} \right) \leqslant a Q_{\max}^{\alpha}$ for each $\bm{Q} \in \mathcal{S}$ ;
\item $\mathbb{E}_{\pi_0} \left[\Phi_0 \left( \bm{Q}(t+1) \right) - \Phi_0 \left( \bm{Q}(t) \right) \mid \bm{Q}(t) = \bm{Q} \right] \leqslant -\epsilon_0$ for each  $\bm{Q}$ such that $Q_{\max} \geqslant B_0$.
\end{enumerate}

\end{assumption}



The requirement of MDP respecting Lyapunov function under a stable policy is not restrictive. For a stabilizable stochastic network, there exists a policy $\pi_0$ under which the corresponding Markov chain is positive recurrent. Moreover, the property of positive recurrence is known to be equivalent to the existence of the so-called Lyapunov function~\cite{mertens1978necessary}. These observations motivate us to focus on MDPs that satisfy Assumption \ref{Asp:KnownPolicy}. In practice, for stabilizable networks, numerous stable policies have been proposed, including dynamic server allocation problems \cite{baras1985k, buyukkoc1985cmu, tassiulas1993dynamic}, multiclass routing networks \cite{bertsimas1994optimization, humes1997delay, kumar1994performance, kumar1995stability, bertsimas2001performance}, switch scheduling \cite{shah2012switched,keslassy2001analysis,tassiulas1993dynamic} and inventory control \cite{neely2010dynamic, dai2005maximum}, which all satisfy Assumption \ref{Asp:KnownPolicy}.

\medskip

\noindent\textbf{Assumption on $\tilde{\pi}^*$: } From the outline of performance analysis in Section~\ref{Sec:Outline}, we need to show that the probability of queue backlog entering $\mathcal{S}^{\mathrm{out}}$ decays exponentially with respect to $U$. We therefore make a natural assumption on the stability property of $\tilde{\pi}^*$. This assumption is necessary for the proof of Lemma \ref{Lem:pBd}. For clarity, we use $\tilde{\mathbb{E}}_{\tilde{\pi}}[\cdot]$ to denote the expectation with respect to the randomness of the auxilary system $\tilde{\mathcal{M}}$ under the policy $\tilde{\pi}$ (to distinguish from $\mathbb{E}[\cdot]$ for $\mathcal{M}$).
\begin{assumption} \label{Asp:DriftOptimalPolicy}
There exist a Lyapunov function $\tilde{\Phi}^*: \tilde{\mathcal{S}}\rightarrow \mathbb{R}_+$ and constants $\beta \geqslant 1$, $b_1, b_2, \tilde{B}^*, \tilde{\epsilon}^* > 0$, such that for any $U > 0$, the following properties hold:
\begin{enumerate}
\item $Q_{\max}^{\beta} \leqslant \tilde{\Phi}^* \left( \bm{Q} \right) \leqslant b_1 Q_{\max}^{\beta}$ for each $\bm{Q} \in \tilde{\mathcal{S}}$ ;
\item $\max_{\bm{Q}' \in \mathcal{R} \left( \bm{Q} \right)} \abs{\tilde{\Phi}^* \left( \bm{Q}' \right) - \tilde{\Phi}^* \left( \bm{Q} \right)} \leqslant b_2 \left( Q_{\max}^{\beta-1} + Q_{\max}^{'\beta-1} \right)$ for each $\bm{Q} \in \tilde{\mathcal{S}}$ ;
\item $\tilde{\mathbb{E}}_{\tilde{\pi}^*} \left[\tilde{\Phi}^* \left( \bm{Q}(t+1) \right) - \tilde{\Phi}^* \left( \bm{Q}(t) \right) \mid \bm{Q}(t) = \bm{Q} \right] \leqslant -\tilde{\epsilon}^* \cdot Q_{\max}^{\beta-1}$ for each $\bm{Q}\in \tilde{\mathcal{S}}$ such that $Q_{\max} \geqslant \tilde{B}^*$.
\end{enumerate}
\end{assumption}

The assumption is natural in the sense that optimal policies are stable and thus are likely to have good Lyapunov drift properties. For instance, if $\tilde{\Phi}^* \left( \bm{Q} \right) = \sum_{i=1}^D \omega_i Q_i^2$ has a negative drift for each $\bm{Q}$ satisfying $\sum_{i=1}^D Q_i \geqslant B$, one can easily show that $\omega_{min} \cdot Q_{\max}^2 \leqslant \tilde{\Phi}^* \left( \bm{Q} \right) \leqslant(\sum_{i=1}^D \omega_i) \cdot Q_{\max}^2$, $\max_{\bm{Q}' \in \mathcal{R} \left( \bm{Q} \right)} \abs{\tilde{\Phi}^* \left( \bm{Q}' \right) - \tilde{\Phi}^* \left( \bm{Q} \right)} \leqslant W \cdot(\sum_{i=1}^D \omega_i) \cdot (Q_{\max} + Q'_{\max})$, and has a negative drift for each $\bm{Q}$ with $Q_{\max} \geqslant B$. Possible methods to analyze Lyapunov drift properties from queueing stability can be found in \cite{meyn2012markov}. 

\medskip

\noindent\textbf{Assumption on communication properties: } Existing model-based reinforcement learning algorithms with performance guarantees usually require the ``diameter'' (i.e. the upper bound for the shortest first hitting time between any two states) of the MDP to be finite \cite{jaksch2010near}. However, it is unrealistic to assume the MDP diameter to be finite when the state space is unbounded. Instead, we make the following assumption (this assumption is used in the proof of Lemma \ref{Lem:EpiExBacklog}) where $T_{\bm{Q} \to \bm{Q}'}$ denotes the first hitting time from $\bm{Q}$ to $\bm{Q}'$.

\begin{assumption} \label{Asp:PolyTrans}
There exist constants $c, \gamma > 0$, such that for any $U > 0$ and every $\bm{Q}, \bm{Q}' \in \tilde{\mathcal{S}}$,
$$
\min_{\tilde{\pi}} \tilde{\mathbb{E}}_{\tilde{\pi}} \left[ T_{\bm{Q} \to \bm{Q}'} \right] \leqslant c \lVert \bm{Q}' - \bm{Q} \rVert_1^{\gamma} ,
$$
where $\tilde{\pi}$ is a policy that can be applied to $\tilde{S}$.
\end{assumption}

In other words, Assumption~\ref{Asp:PolyTrans} states that there exists a policy $\tilde{\pi}$ such that the first hitting time between two states $\bm{Q}$ and $\bm{Q}'$ in $\tilde{\mathcal{S}}$ is a polynomial function of $\lVert \bm{Q}' - \bm{Q} \rVert_1$. We emphasize that the constant $\gamma$ is independent of the value of $U$. Indeed, with larger $U$'s, $\tilde{\mathbb{E}}_{\tilde{\pi}} \left[ T_{\bm{Q} \to \bm{Q}'} \right]$ might grow larger, since the trajectory from $\bm{Q}$ to $\bm{Q}'$ are more likely to be longer. However, it is possible that the increment is bounded. For instance, certain $M/G/1$ queueing systems \cite{ross1999hitting} and random walk models \cite{khantha1983first} have constant $\gamma$'s even if the state space is unbounded. Directly verifying Assumption~\ref{Asp:PolyTrans} might be computationally difficult. Theorem 11.3.11 from \cite{meyn2012markov} offers a drift analysis approach for justifying Assumption~\ref{Asp:PolyTrans}. 

We remark that the maximum first passage time of the MDP plays an essential role in analyzing the reinforcement learning algorithms in infinite-horizon average-reward MDPs. Existing related works either use MDP diameter \cite{jaksch2010near,zhang2019regret}, or span \cite{bartlett2012regal,ouyang2017learning,fruit2018efficient}, or mixing time \cite{ortner2020regret,abbasi2019politex,wei2020model} to characterize the maximum first passage time of the MDP, where these metrics serve as coefficients in the gap to the optimum. The analysis in our paper also relies on
the assumption of reasonable bounds for the maximum first passage time. 
It would be of great interest to establish performance guarantees under relaxed assumptions.

\medskip

\noindent \textbf{Assumption on error tolerance: } As the learning process proceeds, ideally the estimation for $\tilde{\mathcal{M}}$ becomes increasingly accurate. However, it is unrealistic for us to obtain the exact $\tilde{\mathcal{M}}$. Therefore, we assume that if the estimates for the state-transition kernels are accurately enough (i.e., within a certain error bound), the optimal policy to the estimated $\tilde{\mathcal{M}}$ is the same as $\tilde{\pi}^*$. The assumption is stated as follows.
\begin{assumption} \label{Asp:ErrDist}
There exists a $\Delta p > 0$, such that for any MDP $\tilde{\mathcal{M}}'$ with the same state space, action space and cost function as $\tilde{\mathcal{M}}$, if for each $\bm{Q} \in \tilde{\mathcal{S}}$ and each $a \in \mathcal{A}$, we have
$$
\norm{ \tilde{p} \left( \cdot \mid \bm{Q}, a \right) - \tilde{p}' \left( \cdot \mid \bm{Q}, a \right) }_1 \leqslant \Delta p ,
$$
then the optimal policy for $\tilde{\mathcal{M}}'$ is the same as the optimal policy $\tilde{\pi}^*$ for $\tilde{\mathcal{M}}$.
\end{assumption}

In this paper we focus on network control problems with finite action space. For many networks, when the dynamics (e.g. exogenous arrival rates, service rates, channel capacities) vary slightly, the optimal policies remain the same. For instance, if the arrival rates of the switching networks vary yet remain heavy loads, Maximum Matching still remains a close-to-optimal policy \cite{maguluri2016heavy}. 


\subsection{Algorithm} \label{Sec:Algorithm}

We propose an algorithm called Reinforcement Learning for Queueing Networks (RL-QN). 
RL-QN operates in an episodic manner: at the beginning of episode $k$, we uniformly draw a real number $\xi \in [0, 1]$ to decide whether to explore or exploit during this episode. The length of each episode depends on the observations.
\begin{itemize}
    \item If $\xi \leqslant \epsilon_k \triangleq \ell / \sqrt{k} $ (where $\ell \in (0,1]$ can be tuned to control the exploration frequency), we perform exploration during this episode. For states in $\tilde{\mathcal{S}}$, we apply the randomized policy $\pi_{\mathrm{rand}}$. For states in $\mathcal{S} \setminus \tilde{\mathcal{S}}$, we apply $\pi_0$. 
    
    \item If $\xi > \epsilon_k$, we enter the exploitation stage. We first calculate sample means to estimate the state-transition function $\tilde{p}$ of $\tilde{\mathcal{M}}$. We then apply value iteration on the estimated system $\tilde{M}_k$ to obtain an estimated optimal policy $\tilde{\pi}^*_k$.  During this episode, we apply $\tilde{\pi}^*_k$ for states in $\mathcal{S}^{\mathrm{in}}$ and $\pi_0$ otherwise. 
    
    \item When the number of visits to states in $\mathcal{S}^{\mathrm{in}}$ exceeds $L_k = L \cdot \sqrt{k}$ (where $L$ is a positive constant and can be tuned to adjust the sampling rate), RL-QN enters episode $k+1$ and repeat the process above.
\end{itemize}

The details are presented in Algorithm \ref{Alg:RL-QN}. We use $\min \left\{ U, \bm{Q} \right\}$ to denote a vector with the $i$-th coordinate being $\min\{U,Q_i\}.$

\begin{algorithm}
\caption{The RL-QN algorithm} \label{Alg:RL-QN}
\begin{algorithmic}[1]
  \STATE \textbf{Input: } $\mathcal{A}$, $U > W$, $0 < \ell \leqslant 1$, $L > 0$, $K > 0$


  \FOR{episodes $k \leftarrow 1, 2, \cdots,K$}

    \STATE Set $L_k \leftarrow  L \cdot \sqrt{k}$, $\epsilon_k \leftarrow \ell / \sqrt{k}$ and uniformly draw $\xi \in [0, 1]$.

    \IF{$\xi \leqslant \epsilon_k$} 
      \STATE Set 
              $$
              \pi_k \left( \bm{Q} \right) = 
              \begin{cases}
              \pi_{\mathrm{rand}} \left( \bm{Q} \right) , \text{ for } \bm{Q} \in \tilde{\mathcal{S}} \\
              \pi_0 \left( \bm{Q} \right) , \text{ for } \bm{Q} \in \mathcal{S} \setminus \tilde{\mathcal{S}}
              \end{cases} .
              $$  
    \ELSE
      \STATE For each $\bm{Q}, \bm{Q}' \in \tilde{\mathcal{S}}$ and $a \in \mathcal{A}$, let $\tilde{p}_k \left( \bm{Q}' \mid \bm{Q}, a \right) = \tilde{P} \left( \bm{Q}, a, \bm{Q}' \right) / N \left( \bm{Q}, a \right)$ for $N \left( \bm{Q}, a \right) > 0$ and $\tilde{p}_k \left( \bm{Q}' \mid \bm{Q}, a \right) = 1 / \ \abs{\mathcal{R} \left( \bm{Q} \right)}$ otherwise.
      \STATE Solve the estimated MDP $\tilde{M}_k$ and obtain the estimated optimal policy $\tilde{\pi}^*_k$.
      \STATE Set 
              $$
              \pi_k \left( \bm{Q} \right) = 
              \begin{cases}
              \tilde{\pi}^*_k \left( \bm{Q} \right) , \text{ for } \bm{Q} \in \mathcal{S}^{\mathrm{in}} \\
              \pi_0 \left( \bm{Q} \right) , \text{ for } \bm{Q} \in \mathcal{S}^{\mathrm{out}}
              \end{cases} .
              $$  
    \ENDIF

    \WHILE{visits to states in $\mathcal{S}^{\mathrm{in}}$ is smaller than $L_k$}
      \STATE Take the action $a_t = \pi_k \left( \bm{Q}(t) \right)$ for the real system.
      \STATE Observe the next state $\bm{Q}(t+1)$.
      \IF{$\bm{Q}(t) \in \tilde{\mathcal{S}}$} 
        \STATE Increase $N \left( \bm{Q}(t), a_t \right)$ by $1$.
        \STATE Increase $\tilde{P} \left( \bm{Q}(t), a_t, \min \left\{ U, \bm{Q} \right\} \right)$ by $1$. 
      \ENDIF
      \STATE $t \leftarrow t+1$.
    \ENDWHILE

  \ENDFOR

  \STATE \textbf{Output: } estimated optimal policy $\tilde{\pi}^*_K$

\end{algorithmic}
\end{algorithm}





\section{Performance Analysis} \label{Sec:Theory}

We analyze the performance of our algorithm from both exploration and exploitation perspectives, under Assumptions \ref{Asp:KnownPolicy}-\ref{Asp:ErrDist}. We first prove that RL-QN can learn $\tilde{\pi}^*$ within finite episodes with high probability, which implies that RL-QN explores different states sufficiently to obtain an accurate estimation of $\tilde{\mathcal{M}}$ (cf. Theorem~\ref{Thm:LearningProcess}). We then show that RL-QN exploits the estimated optimal policy and achieves a performance close to the optimal result of $\rho^*$ (cf. Theorem~\ref{Thm:Backlog}).

In this paper, we focus on MDPs such that all states are accessible from each other under the following policies: (a) $\pi_{\mathrm{rand}}+\pi_0:$ applying $\pi_{\mathrm{rand}}$ to $\tilde{\mathcal{S}}$ and $\pi_0$ to $\mathcal{S} \setminus \tilde{\mathcal{S}}$; (b) $\tilde{\pi}^*+\pi_0:$ applying $\tilde{\pi}^*$ to $\mathcal{S}^{\mathrm{in}}$ and $\pi_0$ to $\mathcal{S}^{\mathrm{out}}$; and (c) $\pi^*:$ applying (a truncated version of) $\pi^*$ to $\tilde{\mathcal{S}}$ in $\tilde{\mathcal{M}}$. That is, the corresponding Markov chains under the above policies are irreducible.


\subsection{Convergence to the Optimal Policy} \label{Sec:Conv2Policy}

The following theorem states that, with arbitrarily high probability, RL-QN learns $\tilde{\pi}^*$ within a finite number of episodes. 


\begin{theorem} \label{Thm:LearningProcess} Suppose Assumption~\ref{Asp:ErrDist} holds. For each $\delta\in(0,1)$, there exists $k^* < \infty$ such that \emph{RL-QN} learns $\tilde{\pi}^*$ (i.e. $\tilde{\pi}^*_{k^*} = \tilde{\pi}^*$) within $k^*$ episodes with probability at least $1 - \delta$.
\end{theorem}

\begin{proof} We give the detailed proof in Appendix~\ref{App:PfLearningProcess}. Let us now sketch the main steps in our analysis. By analyzing the mixing property of the underlying Markov chain of applying $\pi_{\mathrm{rand}}$ to $\tilde{\mathcal{S}}$ and $\pi_0$ to $\mathcal{S} \setminus \tilde{\mathcal{S}}$, we first show that there exists a constant $K_0$ such that for episode $k \geqslant K_0$, we are able to obtain $\Theta(\sqrt{k})$ samples for each $(\bm{Q}, a)$ pair.

We then use a result in \cite{weissman2003inequalities} to compute the number of samples required for each $(\bm{Q}, a)$ pair to ensure $\tilde{\pi}^*_k = \tilde{\pi}^*$, under Assumption~\ref{Asp:ErrDist}. 
Finally, we show that by choosing the exploration probability of $\epsilon_k = \ell / \sqrt{k}$, exploration episodes occur infinitely often, which implies that we can reach the number of required exploration episodes within finite number of episodes. 


\end{proof}

Theorem \ref{Thm:LearningProcess} indicates that RL-QN explores (i.e. samples) state-transition functions of each state-action pair $(\bm{Q}, a)$ in $\tilde{\mathcal{M}}$ sufficiently.


\subsection{Gap to Optimum} \label{Sec:Backlog}

Theorem \ref{Thm:LearningProcess} states the sufficient exploration aspect of RL-QN. In reinforcement learning, the trade-off between exploration and exploitation is of significant importance to the algorithm performance. In this section, we show that RL-QN also exploits the learned policy such that the episodic average queue backlog is bounded and can get arbitrarily close to $\rho^*$, the optimal average queue backlog of $\mathcal{M}$, as we increase $U$.

We denote the time step at the end of episode $k$ by $t_k$ and the actual length of episode $k$ by $L_k'$, i.e., $L_k' = t_k - t_{k-1}$ with $t_0 \triangleq 0$. We use $\pi_k^{\mathrm{in}}$ to represent the policy applied to $\mathcal{S}^{\mathrm{in}}$ during episode $k$ and $p^{\tilde{\pi} + \pi_0}(\cdot)$ to denote the stationary distribution of states when applying $\tilde{\pi}$ to states in $\mathcal{S}^{\mathrm{in}}$ and $\pi_0$ to states in $\mathcal{S}^{\mathrm{out}}$.

By Theorem \ref{Thm:LearningProcess}, RL-QN learns $\tilde{\pi}^*$ with high probability. Note that as the exploration probability decays by $1/\sqrt{k}$, the probability of utilizing the learned policy converges to $1$ as the episodes increase. Hence the episodic average queue backlog when $\pi_k^{\mathrm{in}} = \tilde{\pi}^*$ constitutes a large proportion of the overall expected average queue backlog. Therefore, the key step to upper bound the expected average queue backlog is to upper bound the episodic average queue backlog when $\pi_k^{\mathrm{in}} = \tilde{\pi}^*$. We prove that it can be upper bounded with respect to $\tilde{\rho}^*$, the optimal average queue backlog of $\tilde{\mathcal{M}}$, as stated in Lemma~\ref{Lem:EpiExBacklog}. We first define $\mathcal{S}^{\mathrm{in}}_{\mathrm{bd}}$ as the set of states in $\mathcal{S}^{\mathrm{in}}$ that is possible to exit into $\mathcal{S}^{\mathrm{out}}$ at the next time slot, i.e. $\mathcal{S}^{\mathrm{in}}_{\mathrm{bd}} \triangleq \left\{ \bm{Q} \in \mathcal{S}^{\mathrm{in}}: \mathcal{R} \left( \bm{Q} \right) \cap \mathcal{S}^{\mathrm{out}} \neq \emptyset \right\}$.

\begin{lemma} \label{Lem:EpiExBacklog} Under Assumptions~\ref{Asp:KnownPolicy}-\ref{Asp:PolyTrans}, we have
$$
\lim_{k \to \infty} \mathbb{E}_{\tilde{\pi}^* + \pi_0} \left[ \frac{\sum_{t = t_{k-1}+1}^{t_k} \sum_i Q_i(t)}{L_k'} \right] = \tilde{\rho}^* + p^{\tilde{\pi}^* + \pi_0} \left( \mathcal{S}^{\mathrm{in}}_{\mathrm{bd}} \right) \cdot \mathcal{O} \left( U^{1 + \max \{ 2\alpha, \gamma \}} \right) .
$$
\end{lemma}

\begin{proof}
We first define the accumulated regret regarding $\tilde{\rho}^*$ for a given episode $k$ with $\pi_k^{\mathrm{in}} = \tilde{\pi}^*$ as $\sum_{t = t_{k-1}+1}^{t_k} \big(\sum_i Q_i(t) - \tilde{\rho}^*\big)$. We then define $\mathcal{S}^{\mathrm{in}}_{\mathrm{in}}  \triangleq \mathcal{S}^{\mathrm{in}} \setminus \mathcal{S}^{\mathrm{in}}_{\mathrm{bd}}$ and decompose the expected accumulated regret into three parts according to state position: (i) $\bm{Q} \in \mathcal{S}^{\mathrm{in}}_{\mathrm{in}}$, (ii) $\bm{Q} \in \mathcal{S}^{\mathrm{in}}_{\mathrm{bd}}$ and (iii) $\bm{Q} \in \mathcal{S}^{\mathrm{out}}$. 

For the first part (i), we use Bellman equation analysis to show that, every time $\bm{Q}$ enters $\mathcal{S}^{\mathrm{in}}_{\mathrm{in}}$ from $\mathcal{S}^{\mathrm{in}}_{\mathrm{bd}}$ and returns back to $\mathcal{S}^{\mathrm{in}}_{\mathrm{bd}}$, the expected accumulated regret is upper bounded by the span of the solutions to the Bellman equations. We then use Proposition 5.5.1 in \cite{bertsekas2017dynamic} to obtain a polynomial upper bound for the span of 
the solution to the Bellman equation under Assumption~\ref{Asp:PolyTrans}.

The second part (ii) is trivially upper bounded by $DU$. 

For the third part (iii), we  prove that the time it takes to return back $\mathcal{S}^{\mathrm{in}}$ is polynomial in $U$ under  Assumption~\ref{Asp:KnownPolicy} (using techniques as the proof of Theorem 1.1 in Chapter 5 of \cite{bremaud1999lyapunov}). Therefore, by Theorem 6.3.4 in \cite{hou2012homogeneous}, every time $\bm{Q}$ enters $\mathcal{S}^{\mathrm{out}}$ from $\mathcal{S}^{\mathrm{in}}_{\mathrm{bd}}$ and returns back to $\mathcal{S}^{\mathrm{in}}_{\mathrm{bd}}$, the incurred expected accumulated regret is upper bounded by a polynomial function of $U$.

The detailed proof is given in Appendix~\ref{App:PfEpiExBacklog}.
\end{proof}

We then proceed to upper bound $p^{\tilde{\pi}^* + \pi_0} \left( \mathcal{S}^{\mathrm{in}}_{\mathrm{bd}} \right),$ as stated in the following lemma.

\begin{lemma} \label{Lem:pBd} Under Assumption~\ref{Asp:DriftOptimalPolicy}, we have
$$
p^{\tilde{\pi}^* + \pi_0} \left( \mathcal{S}^{\mathrm{in}}_{\mathrm{bd}} \right) = \mathcal{O} \left( \exp \left( - U \right) \right) .
$$
\end{lemma}

\begin{proof}
We show that under Assumption~\ref{Asp:DriftOptimalPolicy}, we can construct a linear Lyapunov function with a negative drift for states with large $\bm{Q}_{\max}$. By applying similar techniques as Lemma 1 in \cite{bertsimas1998geometric}, we establish an upper bound for the tail probability of Lyapunov values, which decays exponentially.

The detailed proof is given in Appendix~\ref{App:pBd}.
\end{proof}

With Theorem~\ref{Thm:LearningProcess}, Lemma \ref{Lem:EpiExBacklog} and Lemma \ref{Lem:pBd}, we establish the following main result of this paper.

\begin{theorem} \label{Thm:Backlog} Suppose Assumptions~\ref{Asp:KnownPolicy}-\ref{Asp:ErrDist} hold. When applying \emph{RL-QN} to $\mathcal{M}$, the asympototic episodic average queue backlog is upper bounded as follows:
\begin{align}
 \lim_{k \to \infty} \mathbb{E} \left[ \frac{\sum_{t = t_{k-1}+1}^{t_k} \sum_i Q_i(t)}{L_k'} \right] \leq \rho^* + \mathcal{O} \left( \frac{U^{1 + \max \{ 2\alpha, \gamma \}}}{\exp \left( U \right) } \right) .  \label{eq:episodic_bound} 
\end{align}
\end{theorem}

\begin{proof}
With Lemma~\ref{Lem:EpiExBacklog} and Lemma~\ref{Lem:pBd}, we show that when $k \to \infty$, for each episode with $\pi_k^{\mathrm{in}} = \tilde{\pi}^*$, the expected average queue backlog is upper bounded by $\tilde{\rho}^* + \mathcal{O} (U^{1 + \max \{ 2\alpha, \gamma \}} / \exp (U))$. We next show that the expected average queue backlog contributed by episodes where $\pi_k^{\mathrm{in}} \neq \tilde{\pi}^*$ becomes negligible as $k \to \infty$. By applying similar techniques as the proof of Lemma~\ref{Lem:EpiExBacklog} and Lemma~\ref{Lem:pBd}, we then obtain that $\tilde{\rho}^* = \rho^* + \mathcal{O} ( U^{1 + \gamma } / \exp \left( U \right) )$.

The detailed proof is given in Appendix~\ref{App:PfBacklog}.
\end{proof}

Theorem~\ref{Thm:Backlog} provides an upper bound on the performance guarantee of RL-QN with respect to the threshold parameter $U$: by increasing $U$, the long-term episodic average queue backlog approaches $\rho^*$ exponentially fast. Recall that the episodic length $L'_k$ increases to $\infty$ as $k \rightarrow \infty$. We conjecture that the same upper bound hold for the overall average queue backlog regarding the time horizon $T,$ i.e., 
\begin{align}
\lim_{T \to \infty} \frac{\mathbb{E} \left[\sum_{t = 1}^{T} \sum_i Q_i(t) \right]}{T} \leq \rho^* + \mathcal{O} \left( \frac{U^{1 + \max \{ 2\alpha, \gamma \}}}{\exp \left( U \right) } \right) . \label{eq:conjecture}
\end{align}
We note that the result of Theorem~\ref{Thm:Backlog} does not imply (\ref{eq:conjecture}) directly. A rigorous proof of the conjecture~(\ref{eq:conjecture}) seems difficult with current techniques. We leave as an interesting future direction to investigate if (\ref{eq:conjecture}) holds.


\subsection{Complexity Analysis} \label{Sec:Complexity}

Here we present the complexity analysis. Our algorithm requires solving an estimated MDP for each exploitation episode. Since episode $k$ has length of at least $L \cdot \sqrt{k}$, and
$$
  \sum_{k=1}^{K} \sqrt{k} \geqslant \frac{2}{3} K^{\frac{3}{2}} ,
$$
we have that given a time horrzon $T$, the number of episodes is upper bounded as
$$
  K_T \leqslant \bigg( \frac{3T}{2L} \bigg)^{\frac{2}{3}} .
$$
It has been shown that to solve a general MDP, the time complexity is at least polynomial in the size of the state space~\cite{littman2013complexity}. The concrete time complexity depends on the solution method and its parameters. In our setting, the state space $\mathcal{S}^{in}$ has a size of $\Theta(U^N)$. Therefore, the time complexity of our algorithm is
$$
\mathcal{O} \Big( \big( T/L \big)^{\frac{2}{3}} \cdot \text{poly} \big( U^N \big) \Big) .
$$

Therefore, for a given system with fixed number of nodes, the time complexity grows at most polynomially in the threshold $U$. When applying our algorithm to problems of larger scales (larger $N$), the complexity suffers from ``curse of dimensionality''. However, our algorithm remains computationally feasible in practice for the following reasons. First, although the ``curse of dimensionality'' persists, our algorithm substantially simplifies the original MDP with unbounded state space to an MDP with finite state space. For instance, the system in Section 5.1 cannot be optimized by traditional methods but is optimized efficiently under our algorithm. Second, we only require solving the estimated MDP sparsely, i.e., no more than $(3T/2L)^{2/3}$ times. Choosing a relatively large $L$ can greatly reduce the time complexity. Third, solving the MDP is independent of other steps of our algorithm, which allows us to employ appropriate methods to solve MDPs according to different application scenarios, computational capacities and performance requirements. For instance, since we aim at optimizing the average cost of the MDP, undiscounted value iteration/policy iteration methods should be applied. However, we apply discounted methods in our numerical experiments and obtain the same solution with significantly less computation time. Also, we may utilize the special structure of the MDP \cite{guestrin2003efficient,koller2013policy} or apply approximation methods \cite{chen1999value,chen2018deep,veatch2010approximate} to alleviate the computational cost. Our numerical results in Section \ref{Sec:SA_2} further validates the conclusion and shows that our algorithm is computationally feasible in practice.





\section{Numerical Experiments} \label{Sec:Sim}

In this section, we evaluate the performance of RL-QN for three classes of queueing systems:  server allocation, routing, switching.


\subsection{Server Allocation (Two Nodes)} \label{Sec:SA_2}

We consider a dynamic server allocation problem: exogenous packets arrive to two nodes according to Bernoulli process with rate $\lambda_1$ and $\lambda_2$ respectively. At each time slot, the central server selects one node to serve. The head of line job in the selected queue $i$ completes the required service and leaves the system with probability $p_i$. The system model and parameters are illustrated in Figure~\ref{Fig:SA_2_Model}.

\begin{figure}[H]
\centering
\begin{minipage}{.5\textwidth}
  \centerline{\includegraphics[width=0.9\linewidth]{SA_2_Model}}
  \caption{System model of the dynamic server allocation model with two nodes.}
  \label{Fig:SA_2_Model}
\end{minipage}%
\begin{minipage}{.5\textwidth}
  \centerline{\includegraphics[width=0.9\linewidth]{SA_2_HT}}
  \caption{Simulation results for the dynamic server allocation model of two nodes.}
  \label{Fig:SA_2_HT}
\end{minipage}
\end{figure}

According to \cite{tassiulas1993dynamic}, whenever the condition
$$
\frac{\lambda_1}{p_1} + \frac{\lambda_2}{p_2} < 1 ,
$$
is satisfied, one stabilizing policy is to always serve the node with the longest queue (LQ). Therefore, we can use LQ policy as $\pi_0$. To evaulate whether this problem satisfies Assumption \ref{Asp:PolyTrans}, we apply $\pi_0$ to the truncated state space $\tilde{\mathcal{S}}$ with $U = 10$ and simulate the system to collect the hitting times, as shown in Figure \ref{Fig:SA_2_HT}. We can see that as the state distance grows, the average first hitting time grows sublinearly, which is obviously upper bounded by linear growth and thus indicates that Assumption \ref{Asp:PolyTrans} is satisfied. Moreover, our numerical experiments show that our algorithm can optimize the average queue backlog for various settings even if we cannot provably verify Assumption \ref{Asp:PolyTrans}. Therefore, the applicability of our algorithm may not be sensitive to Assumption \ref{Asp:PolyTrans}, which helps to show that our algorithm is practical.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.5\linewidth]{SA_2_Result}}
\caption{Average queue backlog under different policies for the server allocation model with two nodes.}
\label{Fig:SA_2_Result}
\end{figure}

We simulate RL-QN for $U=5$, $U=10$ and $U=20$. The results are shown in Figure~\ref{Fig:SA_2_Result}. To obtain a lower bound on the average queue backlog, we assume that the system dynamics, i.e., arrival rates and success rates, are known. We then optimize the MDP and obtain its average queue backlog. This average queue backlog is guaranteed to be a lower bound since, in practice, the system dynamics are unknown and errors that occur during the learning process can downgrade the performance. From Figure~\ref{Fig:SA_2_Result}, we see that the longest queue (LQ) policy stabilizes the queue backlog, yet its average queue backlog is far from the lower bound. All of our RL-QN methods outperform $\pi_0$. When $U = 5$, the average queue backlog converges to $2.38$, while for $U = 10$ and $U = 20$, the average queue backlog becomes $2.24$. This indicates that RL-QN achieves better performance with a larger threshold parameter $U$, as implied by Theorem~\ref{Thm:Backlog}. Moreover, since the cases with $U = 10$ and $U = 20$ achieve similar performance, it is very likely that in practice a small $U$ achieves satisfactory performance, which makes our algorithm more computationally feasible.

\begin{figure}[H]
\centering
\begin{minipage}{.48\textwidth}
  \centerline{\includegraphics[width=0.95\linewidth]{SA_2_TU}}
  \caption{The relationship between $U$ and the relative \\ computation time.}
  \label{Fig:SA_2_TU}
\end{minipage}%
\begin{minipage}{.48\textwidth}
  \centerline{\includegraphics[width=0.95\linewidth]{SA_2_QT}}
  \caption{The relationship between the relative computation time \\ and the average queue backlog.}
  \label{Fig:SA_2_QT}
\end{minipage}
\end{figure}

We then study the time complexity of our algorithm. We apply three different MDP solvers in our algorithm: value iteration, relative value iteration and policy iteration, all with discount factor $0.99$. Figure \ref{Fig:SA_2_TU} shows the relationship between the relative computation time (i.e., the actual computation time of RL-QN devided by the actual computation time of LQ) and the average queue backlog. Even when $U = 20$, the computation time is less than three times of the LQ policy, which shows that increasing $U$ has relatively small impact on the computation time. From Figure \ref{Fig:SA_2_QT}, we see that the average queue backlog drops significantly even with a relatively small $U$. The average queue backlog is reduced by $23.1\%$ when $U = 8$, while the computation time is only $1.17$ times of the traditional non-ML LQ policy. Therefore, our algorithm is computationally efficient in practice.




\subsection{Server Allocation (Ten Nodes)} \label{Sec:SA_10}

We then consider a dynamic server allocation problem with greater scale: exogenous packets arrive to ten nodes according to Bernoulli process with rate $\lambda_i, \ i = 1, \cdots, 10$ respectively. At each time slot, the central server selects one node to serve. The head of line job in the selected queue $i$ completes the required service and leaves the system with probability $p_i$. The system model and parameters are illustrated in Figure~\ref{Fig:SA_10_Model}.

\begin{figure}[h]
\centerline{\includegraphics[width=0.9\linewidth]{SA_10_Model}}
\caption{System model of the dynamic server allocation model with ten nodes.}
\label{Fig:SA_10_Model}
\end{figure}

As discussed in Section \ref{Sec:Outline}, if we directly apply classical methods to solve the estimated MDP, the computational complexity is at least $\Theta(U^N)$, which grows exponentially with the number of nodes $N$. Therefore, for the ten-node dynamic server allocation model, we prefer to utilize the special structures of the MDP or apply approximation methods to reduce the computational cost.

It has been shown in \cite{buyukkoc1985cmu} that when all nodes are always connected to the central server, the policy to minimize the average job delay is to serve the node with the largest service rate $p_i$ among all non-empty nodes. Therefore, for each exploitation episode, instead of solving the estimated MDP $\tilde{\mathcal{M}}_k$ with classical methods, we directly obtain the optimal solution: serving the node with the largest estimated $p_i$ among all non-empty nodes. We emphasize that our algorithm does not depend on how the estimated MDPs are solved. The above computational trick utilizes the special structure of the server allocation model to speed up the computation, and other techniques for efficiently solving the MDPs can also potentially be employed.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.5\linewidth]{SA_10_Result}}
\caption{Simulation results for the server allocation model with ten nodes.}
\label{Fig:SA_10_Result}
\end{figure}

The results are shown in Figure~\ref{Fig:SA_10_Result}. From the figure, we observe that RL-QN outperforms LQ. The LQ reaches an average queue backlog of $7.6$, while RL-QN reaches $7$ when $U=5$. When $U$ is beyond $10$, RL-QN reaches an average queue backlog of $6$, which is close to the optimal result.


\subsection{Routing} \label{Sec:Routing_2}

We consider a simple routing problem: exogenous packets arrive at the source node $s$ according to Bernoulli process with rate $\lambda = 0.85$. Node $1$ and node $2$ are two intermediate nodes and can serve at most one packet during each time slot, with probability $p_1$ and $p_2$ respectively. Node $d$ is the destination node. At each time slot, node $s$ has to choose between routes $s \to 1 \to d$ and $s \to 2 \to d$ to transit new exogenous packets. Specifically, the system model and parameters are shown in Figure~\ref{Fig:Routing_2_Model}.

The parameters $(p_1,p_2)$ are queue-dependent here:
$$
(p_1, p_2) = 
\begin{cases}
(0.9, 0.1), & Q_2(t) \leqslant 5 \\
(0.1, 0.9), & Q_2(t) > 5
\end{cases}
.
$$

For each $\lambda < 0.9$, an intuitive stabilizing policy is to always use the fixed path $s \to 1 \to d$, while never choose $s \to 2 \to d$. Therefore, we can use the policy that always routes through $s \to 1 \to d$ as $\pi_0$. However, it is possible that we could split the external arrivals into the two routes to fully utilize the service capacities of both node $1$ and $2$, and achieve better performance.

We simulate RL-QN for $U=10$. The results are plotted in Figure~\ref{Fig:Routing_2_Result}, which shows that RL-QN outperforms the fixed path stabilizing policy and converges to the optimum quickly.

\begin{figure}[H]
\centering
\begin{minipage}{.4\textwidth}
  \centerline{\includegraphics[width=0.9\linewidth]{Routing_2_Model}}
  \caption{System model of the routing network with two nodes.}
  \label{Fig:Routing_2_Model}
\end{minipage}%
\begin{minipage}{.5\textwidth}
  \centerline{\includegraphics[width=0.9\linewidth]{Routing_2_Result}}
  \caption{Simulation results of the routing network with two nodes.}
  \label{Fig:Routing_2_Result}
\end{minipage}
\end{figure}


\subsection{Switching} \label{Sec:Switching}

We consider a $2 \times 2$ input-queued switch as illustrated in Figure~\ref{Fig:Switching_Model}. Data packets arriving at input $i$ destined for output $j$ are stored at input port $i,$ in queue $Q_{i,j}$, thus there are $4$ queues in total. 
We consider the case where the new data packets are arriving
at queue $(i,j)$ at rate $\lambda_{i,j}$ for $1\leq i,j \leq 2$, according to a Bernoulli
process. That is, for each time slot, the number of packets arriving at queue $Q_{i,j}$ is
a Bernoulli random variable with mean $\lambda_{i,j}$.
The server then selects a matching between the inputs and outputs to transmit packets. If input $i$ is connected with output $j$, then a buffered packet is removed from the input queue $Q_{i,j}$ and sent to output $j$.
\begin{figure}[htbp]
\centerline{\includegraphics[width=0.5\linewidth]{Switching_Model}}
\caption{System model of a $2 \times 2$ input-queued cell-switch.}
\label{Fig:Switching_Model}
\end{figure}

According to \cite{mckeown1999achieving}, whenever the condition that $\sum_{i} \lambda_{i, j} < 1, \sum_{j} \lambda_{i, j} < 1$ is satisfied, the Maximum Matching algorithm which selects the matching that maximizes the total queue backlog of the connected channels is stabilizing, hence we use it as $\pi_0$.

\begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.45\textwidth} 
        \centering
        \includegraphics[width=\linewidth]{Switching_asy_Result}
        \caption{$\lambda_{1, 1} = \lambda_{1, 2} = 0.4$, $\lambda_{2, 1} = 0.2$, $\lambda_{2, 2} = 0.1$}
        \label{Fig:Switching_asy_Result}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.45\textwidth} 
        \centering
        \includegraphics[width=\linewidth]{Switching_heavy_Result}
        \caption{$\lambda_{1, 1} = \lambda_{1, 2} = \lambda_{2, 1} = \lambda_{2, 2} = 0.49$}
        \label{Fig:Switching_heavy_Result}
    \end{subfigure}
    \caption{Simulation results for the input-queued switch model.}
    \label{Fig:Switching_Result}
\end{figure}

We implement the simulation for $U=5$ under two different settings of arrival traffic. The results are shown in Figure~\ref{Fig:Switching_Result}. When $\lambda_{1, 1} = \lambda_{1, 2} = 0.4$, $\lambda_{2, 1} = 0.2$, $\lambda_{2, 2} = 0.1$, we can see that our algorithm outperforms the stabilizing policy $\pi_0$ and converges to $\tilde{\pi}^* + \pi_0$. When $\lambda_{1, 1} = \lambda_{1, 2} = \lambda_{2, 1} = \lambda_{2, 2} = 0.49$, the system is under heavy traffic. In \cite{maguluri2016heavy}, it was shown that the Maximum Matching policy $\pi_0$  is close to the optimal policy that minimizes the average queue backlog. Our simulation result (Figure~\ref{Fig:Switching_heavy_Result}) indicates that our algorithm achieves comparable performance as the near-optimal policy $\pi_0.$






\section{Conclusion} \label{Sec:Conclusion}

In this work, we apply a model-based reinforcement learning framework to general queueing networks with unbounded state spaces. We propose the RL-QN algorithm, which applies an $\epsilon$-greedy exploration scheme. By leveraging Lyapunov analysis, we show that the average queue backlog of the proposed approach can get arbitrarily close to the optimal average queue backlog under the optimal policy. The proposed RL-QN algorithm requires the knowledge of a stable policy. An interesting future direction is to investigate this problem when such information is not available.



