Differentiable Discrete Event Simulation for Queuing Network Control
Ethan Che Jing Dong Hongseok Namkoong



=8pt plus0pt minus3pt
=8pt plus0pt minus3pt












  Differentiable Discrete Event Simulation 

  
 for Queuing Network Control  

  
 Ethan Che     Jing Dong     Hongseok Namkoong 

  

  Columbia Business School 

  

  







§ ABSTRACT 
  Queuing network control is essential for managing congestion in job-processing systems 
such as service systems, communication networks, and manufacturing processes. Despite growing interest in applying reinforcement learning (RL) techniques,  queueing network control

poses distinct challenges, including high stochasticity, large state and action spaces, and lack of stability. 
To tackle these challenges, we propose a scalable framework for policy optimization based on differentiable discrete event simulation. 
Our main insight is that by implementing a well-designed smoothing technique for discrete event dynamics, we can compute  policy gradients for large-scale queueing networks using auto-differentiation software (e.g., Tensorflow, PyTorch) and GPU parallelization. Through extensive empirical experiments, we observe that our policy gradient estimators are several orders of magnitude more accurate than typical -based  estimators. 
In addition, we propose a new policy architecture, which drastically improves stability while maintaining the flexibility of neural-network policies.
In a wide variety of scheduling and admission control tasks, we demonstrate that training control policies with pathwise gradients leads to a 50-1000x improvement in sample efficiency over state-of-the-art RL methods.
Unlike prior tailored approaches to queueing, our methods can flexibly handle realistic scenarios, including systems operating in non-stationary environments and those with non-exponential interarrival/service times.




















§ INTRODUCTION





Queuing models are a powerful modeling tool to conduct performance analysis and optimize operational policies in diverse applications such as service systems (e.g., call centers <cit.>, healthcare delivery systems <cit.>, ride-sharing platforms <cit.>, etc), computer and communication systems <cit.>, manufacturing systems <cit.>, and financial systems (e.g., limit order books <cit.>). Standard tools for queuing control analysis involve establishing structural properties of the underlying Markov decision process (MDP) or leveraging analytically more tractable approximations such as fluid <cit.> or diffusion approximations <cit.>. These analytical results often give rise to simple control policies that are easy to implement and interpret. However, these policies only work under restrictive modeling assumptions and can be highly sub-optimal outside of these settings. Moreover, deriving a good policy for a given queuing network model requires substantial queuing expertise and can be theoretically challenging. 



Recent advances in reinforcement learning (RL) have spurred growing interest in applying learning methodologies to solve queuing control problems, which benefit from increased data and computational resources  <cit.>. These algorithms hold significant potential for generating effective controls for complex, industrial-scale networks encountered in real-world applications, which typically fall outside the scope of theoretical analysis. 
However, standard model-free RL algorithms <cit.> often under-perform in queuing control, even when compared to simple queuing policies <cit.>, unless proper modifications are made. This under-performance is primarily due to the unique challenges posed by queuing networks, including (1) high stochasticity of the trajectories, (2) large state and action spaces, and (3) lack of stability guarantees under sub-optimal policies <cit.>.
For example, when applying policy gradient methods, typical policy gradient estimators based on coarse feedback from the environment (observed costs) suffer prohibitive error due to high variability (see, e.g., the  estimators in Figure <ref>). 

To tackle the challenges in applying off-the-shelf RL solutions for queuing control, we propose a new scalable framework for policy optimization that incorporates domain-specific queuing knowledge. Our main algorithmic insight is that queueing networks possess key structural properties that allow for several orders of magnitude more accurate gradient estimation. 
By leveraging the fact that the dynamics of discrete event simulations of queuing networks are governed by observed exogenous randomness (interarrival and service times), we propose a differentiable discrete event simulation framework. This framework enables the computation of a  gradient of a performance objective (e.g., cumulative holding cost) with respect to actions. 

Our proposed gradient estimator, denoted as the  estimator, can then be used to efficiently optimize the parameters of a control policy through stochastic gradient descent (SGD). By utilizing the known structure of queuing network dynamics, our approach provides finer-grained feedback on the sensitivity of the performance objective to any action taken along the sample path. This offers an infinitesimal counterfactual analysis: how the performance metric would change if the scheduling action were slightly perturbed.  
Rather than relying on analytic prowess to compute these gradients, we utilize the rapid advancements in scalable auto-differentiation libraries such as PyTorch <cit.> to efficiently compute gradients over a single sample path or a batch of sample paths. 
Our proposed approach supports very general control policies, including neural network policies, which have the potential to improve with more data and computational resources. 
Notably, our method seamlessly handles large-scale queuing networks and large batches of data via GPU parallelization. Unlike off-the-shelf RL solutions whose performance is exceedingly sensitive to implementation details <cit.>, our method is easy to implement (see e.g., Figure <ref>) and requires minimal effort for parameter tuning.

Across a range of queueing networks, we empirically observe that our  estimator substantially improves the sample efficiency and stability of learning algorithms for queuing network control while preserving the flexibility of learning approaches.  In Figure <ref>, we preview our main empirical findings which show that  gradients lead to a 50-1000x improvement in sample efficiency over model-free policy gradient estimators (e.g.,  <cit.>). 

Buoyed by the promising empirical results, we
provide several theoretical insights explaining the observed efficiency gains. 




Our proposed approach draws inspiration from gradient estimation strategies developed in the stochastic modeling and simulation literature, particularly infinitesimal perturbation analysis (IPA) <cit.>. While IPA has been shown to provide efficient gradient estimators for specific small-scale queuing models (e.g., the G/G/1 queue), it is well-known that unbiased IPA estimates cannot be obtained for general multi-class queuing networks due to non-differentiability of the sample path  <cit.>. Our framework overcomes this limitation by proposing a novel smoothing technique based on insights from fluid models/approximations for queues and tools from the machine learning (ML) literature. To the best of our knowledge, our method is the first to provide a gradient estimation framework capable of handling very general and large-scale queuing networks and various control policies. Our modeling approach is based on discrete-event simulation models, and as a result, it can accommodate non-stationary and non-Markovian inter-arrival and service times, requiring only samples instead of knowledge of the underlying distributions. 

Our second contribution is a simple yet powerful modification to the control policy architecture. It has been widely observed that training a standard RL algorithm, such as proximal policy optimization <cit.> (PPO), may fail to converge due to instabilities arising from training with random initialization. To address this issue, researchers have proposed either switching to a stabilizing policy when instability occurs <cit.> or imitating (behavior cloning) a stabilizing policy at the beginning <cit.>. However, both methods limit policy flexibility and introduce additional complexity in the training process. We identify a key source of the problem: generic policy parameterizations (e.g., neural network policies) do not enforce work conservation, leading to scenarios where even optimized policies often assign servers to empty queues. 
To address this, we propose a modification to standard  policy parameterizations in deep reinforcement learning,

which we refer to as the `work-conserving softmax'. This modification is compatible with standard reinforcement learning algorithms and automatically guarantees work conservation. Although work conservation does not always guarantee stability, we empirically observe across many scenarios that it effectively eliminates instability in the training process, even when starting from a randomly initialized neural network policy.


This modification not only complements our gradient estimator but is also compatible with other model-free RL approaches. We find that while PPO without any modifications fails to stabilize large queuing networks and leads to runaway queue lengths, PPO with the work-conserving softmax remains stable from random initialization and can learn better scheduling policies than traditional queuing policies.

Since rigorous empirical validation forms the basis of algorithmic progress, we provide a thorough empirical validation of the effectiveness of the differentiable discrete event simulator for queuing network control. We construct a wide variety of benchmark control problems, ranging from learning the cμ-rule in a simple multi-class queue to scheduling and admission control in large-scale networks. Across the board, we find that our proposed  gradient estimator achieves significant improvements in sample efficiency over model-free alternatives, which translate to downstream improvements in optimization performance.


  * In a careful empirical study across 10,800 parameter settings, we find that for 94.5% of these settings our proposed  gradient estimator computed along a single sample path achieves greater estimation quality than  with 1000x more data (see section <ref>).

  * In a scheduling task in multi-class queues, gradient descent with  gradient estimator better approximates the optimal policy (the cμ-rule) and achieves a smaller average cost than  with a value function baseline and 1000x more data (see section <ref>).

  * In an admission control task, optimizing the buffer sizes with  gradient estimator achieves smaller costs than randomized finite differences (SPSA <cit.>) with 1000x more data, particularly for higher-dimensional problem instances (see section <ref>).

  * For large-scale scheduling problems, policy gradient with  gradient estimator and work-conserving softmax policy architecture achieves a smaller long-run average holding cost than traditional queuing policies
and state-of-the-art RL methods such as , which use 50x more data (see section <ref>). Performance gains are greater for larger networks with non-exponential noise.

These order-of-magnitude improvements in sample efficiency translate to improved computational efficiency when drawing trajectories from a simulator and improved data efficiency if samples of event times are collected from a real-world system.

Overall, these results indicate that one can achieve significant improvements in sample efficiency by incorporating the specific structure of queuing networks, which is under-utilized by model-free reinforcement learning methods. In section <ref>, we investigate the M/M/1 queue as a theoretical case study and show that even with an optimal baseline,  has a sub-optimally large variance under heavy traffic compared to a pathwise policy gradient estimator. This analysis identifies some of the statistical limitations of , and illustrates that a better understanding of the transition dynamics, rather than narrowly estimating the value-function or Q-function, can deliver large improvements in statistical efficiency. Given the scarcity of theoretical results comparing the statistical efficiency of different policy gradient estimators, this result may be of broader interest.

Our broad aim with this work is to illustrate a new paradigm for combining the deep, structural knowledge of queuing networks developed in the stochastic modeling literature with learning and data-driven approaches.
Rather than either choosing traditional queuing policies, which can be effective for certain queueing control problems but do not improve with data, or choosing model-free reinforcement learning methods, which learn from data but do not leverage known structure, our framework offers a favorable midpoint: we leverage structural insights to extract much more informative feedback from the environment, which can nonetheless be used to optimize black-box policies and improve reliability. Beyond queuing networks, our algorithmic insight provides a general-purpose tool for computing gradients in general discrete-event dynamical systems. Considering the widespread use of discrete-event simulators with popular modeling tools such as AnyLogic <cit.> or Simio <cit.> and open-source alternatives such as SimPy <cit.>, the tools developed in this work can potentially be applied to policy optimization problems in broader industrial contexts.

The organization of this paper is as follows. In section <ref>, we discuss connections with related work. In section <ref>, we introduce the discrete-event dynamical system model for queuing networks. In section <ref>, we introduce our framework for gradient estimation. In section <ref>, we perform a careful empirical study of our proposed gradient estimator, across estimation and optimization tasks.
In section <ref>, we discuss the instability issue in queuing control problems and our proposed modification to the policy architecture to address this. In section <ref>, we empirically investigate the performance of our proposed pathwise gradient estimation and work-conserving policy architecture in optimizing scheduling policies for large-scale networks. In section <ref>, we discuss the M/M/1 queue as a theoretical case study concerning the statistical efficiency of  compared to  estimators. Finally, section <ref> concludes the paper and discusses extensions.












§ RELATED WORK



We discuss connections to related work in queuing theory, reinforcement learning, and gradient estimation in machine learning and operations research.



  
Scheduling in Queuing Networks

Scheduling is a long-studied control task in the queuing literature for managing queues with multiple classes of jobs  <cit.>. Standard policies developed in the literature include static priority policies such as the cμ-rule <cit.>, threshold policies <cit.>, policies derived from fluid approximations <cit.>, including discrete review policies <cit.>, policies that have good stability properties such as MaxWeight  <cit.> and MaxPressure  <cit.>. Many of these policies satisfy desirable properties such as throughput optimality <cit.>, or cost minimization <cit.> for certain networks and/or in certain asymptotic regimes. In our work, we aim to leverage some of the theoretical insights developed in this literature to design reinforcement learning algorithms that can learn faster and with less data than model-free RL alternatives. We also use some of the standard policies as benchmark policies when validating the performance of our  policy gradient algorithm.



  
Reinforcement Learning in Queueing Network Control
Our research connects with the literature on developing reinforcement learning algorithms for queuing network control problems <cit.>. These works apply standard model-free RL techniques (e.g. Q-learning, , value iteration, etc.) but introduce novel modifications to address the unique challenges in queuing network control problems. Our work differs in that we propose an entirely new methodology for learning from the environment based on differentiable discrete event simulation, which is distinct from all model-free RL methods. The works <cit.> observe that RL algorithms tend to be unstable and propose fixes to address this, such as introducing a Lyapunov function into the rewards, or behavior cloning of a stable policy for initialization. In our work, we propose a simple modification to the policy network architecture, denoted as the work-conserving softmax as it is designed to ensure work-conservation. We find empirically that work-conserving softmax ensures stability with even randomly initialized neural network policies. In our empirical experiments, we primarily compare our methodology with the  algorithm developed in <cit.>. In particular, we construct a  baseline with the same hyper-parameters, neural network architecture, and variance reduction techniques as in <cit.>, although with our policy architecture modification that improves stability. 




  
Differentiable Simulation in RL and Operations Research

While differentiable simulation is a well-studied paradigm for control problems in physics and robotics <cit.>, it has only recently been explored for large-scale operations research problems. For instance, <cit.> study inventory control problems and train a neural network using direct back-propagation of the cost, as sample paths of the inventory levels are continuous and differentiable in the actions. In our work, we study control problems for queuing networks, which are discrete and non-differentiable, preventing the direct application of such methods. To address this, we develop a novel framework for computing pathwise derivatives for these non-differentiable systems, which proves highly effective for training control policies. Another line of work, including <cit.>, proposes differentiable agent-based simulators based on differentiable relaxations. While these relaxations have shown strong performance in optimization tasks, they also introduce unpredictable discrepancies with the original dynamics. We introduce tailored differentiable relaxations in the back-propagation process only, ensuring that the forward simulation remains true to the original dynamics. 



  
Gradient Estimation in Machine Learning

Gradient estimation <cit.> is an important sub-field of the machine learning literature, with applications in probabilistic modeling <cit.> and reinforcement learning <cit.>. There are two standard strategies for computing stochastic gradients <cit.>. The first is the score-function estimator or  <cit.>, which only requires the ability to compute the gradient of log-likelihood but can have high variance <cit.>. Another strategy is the reparameterization trick <cit.>, which involves decomposing the random variable into the stochasticity and the parameter of interest, and then taking a pathwise derivative under the realization of the stochasticity. Gradient estimators based on the reparameterization trick can have much smaller variance <cit.>, but can only be applied in special cases (e.g. Gaussian random variables) that enable this decomposition. Our methodology makes a novel observation that for queuing networks, the structure of discrete-event dynamical systems gives rise to the reparameterization trick. Nevertheless, the function of interest is non-differentiable, so standard methods cannot be applied. As a result, our framework also connects with the literature on gradient estimation for discrete random variables <cit.>. In particular, to properly smooth the non-differentiability of the event selection mechanism, we employ the straight-through trick <cit.>, which has been previously used in applications such as discrete representation learning <cit.>. Our work involves a novel application of this technique for discrete-event systems, and we find that this is crucial for reducing bias when smoothing over long time horizons.




  
Gradient Estimation in Operations Research
There is extensive literature on gradient estimation for stochastic systems <cit.>, some with direct application to queuing optimization  <cit.>.
Infinitesimal Perturbation Analysis (IPA) <cit.> is a standard framework for constructing pathwise gradient estimators, which takes derivatives through stochastic recursions that represent the dynamics of the system. While IPA has been applied successfully to some specific queuing networks and discrete-event environments more broadly <cit.>, standard IPA techniques cannot be applied to general queuing networks control problems, as has been observed in <cit.>. There has been much research on outlining sufficient conditions under which IPA is valid, such as the commuting condition in <cit.> or the perturbation conditions in <cit.>, but these conditions do not hold in general. Several extensions to IPA have been proposed, but these alternatives require knowing the exact characteristics of the sampling distributions and bespoke analysis of event paths <cit.>. Generalized likelihood-ratio estimation <cit.> is another popular gradient estimation framework, which leverages an explicit Markovian formulation of state transitions to estimate parameter sensitivities. However, this requires knowledge of the distributions of stochastic inputs, and even with this knowledge, it may be difficult to characterize the exact Markov transition kernel of the system. Finally, finite differences <cit.> and finite perturbation analysis <cit.> are powerful methods, particularly when aided with common random numbers <cit.>, as it requires minimal knowledge about the system. However, it has been observed that performance can scale poorly with problem dimension <cit.>, and we also observe this in an admission control task (see Section <ref>).

Our contribution is proposing a novel, general-purpose framework for computing pathwise gradients through careful smoothing, which only requires samples of random input (e.g., interarrival times and service times) rather than knowledge of their distributions.

Given the negative results about the applicability of IPA for general queuing network control problems (e.g., general queuing network model and scheduling policies), we introduce bias through smoothing to achieve generality. It has been observed in <cit.> that biased IPA surrogates can be surprisingly effective in simulation optimization tasks such as ambulance base location selection. Our extensive empirical results confirm this observation and illustrate that while there is some bias, it is very small in practice, even over long time horizons (>10^5 steps).







§ DISCRETE-EVENT DYNAMICAL SYSTEM MODEL FOR QUEUING NETWORKS






































































































































We describe multi-class queuing networks as discrete-event dynamical systems. This is different from the standard Markov chain representation, which is only applicable when inter-arrival and service times are exponentially distributed.
To accommodate more general event-time distributions, the system description not only involves the queue lengths, but also auxiliary information such as residual inter-arrival times and workloads. Surprisingly, this more detailed system description leads to a novel gradient estimation strategy (discussed in Section <ref>) for policy optimization.

We first provide a brief overview of the basic scheduling problem. We then describe the discrete-event dynamics of multi-class queuing networks in detail and illustrate with a couple of well-known examples. While queuing networks have been treated as members of a more general class of Generalized Semi-Markov Processes (GSMPs) that reflect the discrete-event structure of these systems <cit.>, we introduce a new set of notations tailored for queuing networks to elaborate on some of their special structures.
In particular, we represent the discrete event dynamics 
via matrix-vector notation that maps directly to its implementation in auto-differentiation frameworks, allowing for the differentiable simulation of large-scale queueing networks through GPU parallelization. 



 §.§  The Scheduling Problem



A multi-class queuing network consists of n queues and m servers. The core state variable is the queue lengths associated with each queue, denoted as  x(t) ∈_+^n, which evolves over continuous time. As a discrete-event dynamical system, the state also includes auxiliary data denoted as (t)— consisting of residual inter-arrival times and workloads at time t—which determines state transitions but are typically not visible to the controller.


The goal of the controller is to route jobs to servers, represented by an assignment matrix u ∈{0, 1 }^m × n, to manage congestion. More concretely, the problem is to derive a policy π(x), which only depends on the observed queue lengths and selects scheduling actions, to minimize the integral of some instantaneous costs c(x,u). A typical instantaneous cost is a linear holding/waiting cost:

    c(x,u)=h^⊤x






for some vector h∈_+^n. The objective is to find a
policy π that minimizes the cumulative cost over a time horizon:

    min_π𝔼[∫_0^Tc(x(t),π(x(t)))dt].


Optimizing a continuous time objective can be difficult and may require an expensive discretization procedure. However, discrete-event dynamical systems are more structured in that x(t) is piecewise constant and is only updated when an event occurs. For the multi-class queuing networks, events are either arrivals to the network or job completions, i.e., a server finishes processing a job.

It is then sufficient to sample the system only when an event occurs, and we can approximate the continuous-time objective with a performance objective in the discrete-event system over N events,

    min_π{J_N(π) := 𝔼[∑_k=0^N-1c(x_k,π(x_k))τ_k+1^*]
    }

where x_k is the queue lengths after the kth event update. τ^*_k+1 is an inter-event time that measures the time between the kth and (k+1)th event, and N is chosen such that the time of the Nth event, a random variable denoted as t_N, is “close" to T.





The dynamics of queuing networks are highly stochastic, with large variations across trajectories. Randomness in the system is driven by the random arrival times of jobs and the random workloads (service requirements) of these jobs.  We let ξ_1:N={ξ_i}_i=1^N denote a single realization, or `trace', of these random variables over the horizon of N events.

We can then view the expected cost (<ref>) more explicitly as a policy cost averaged over traces. 
In addition, we focus on a parameterized family of policies {π_θ:θ∈Θ}, for some Θ⊆^d, in order to optimize (<ref>) efficiently. In this case, we utilize the following shorthand J_N(θ;ξ_1:N) for the
policy cost over a single trace and J_N(θ) for the average policy
cost under π_θ,
which leads to the parameterized control problem:

    min_θ{J_N(θ) := 𝔼[J(θ;ξ_1:N)]
    := 
    𝔼[∑_k=0^N-1c(x_k,π_θ(x_k))τ^*_k+1]
    }.

We now turn to describe the structure of the transition dynamics of multi-class queuing networks, to elaborate how scheduling actions affect the queue lengths.






 §.§ System Description

Recall that the multi-class queuing network consists of n queues and m servers, where each queue is associated with a job class, and different servers can be of different compatibilities with various job classes. Recall that x(t)∈ℕ_+^n denotes the lengths of the queues at time t ∈_+.














The queue lengths x(t) are updated by one of two types of events: job arrivals or job completions. 
Although the process evolves in continuous time, it is sufficient to track the system only when an event occurs. We let k∈_+ count the kth event in the system, and let t_k denote the time immediately after the kth event occurs. By doing so, we arrive at a discrete-time representation of the system. Given that we do not assume event times are exponential, the queue lengths x_k alone are not a Markovian descriptor of the system. Instead, we must consider an  augmented state s_k = (x_k, _k), where x_k ∈^n_+ is the vector of queue lengths and _k = (τ_k^A, w_k),∈^2n_+ is an auxiliary state vector that includes residual inter-arrival times τ_k^A = {}_j = 1^n∈^n_+ and residual workloads w_k = {}_j = 1^n∈^n_+ of the `top-of-queue' jobs in each queue. The auxiliary state variables determine the sequence of events.

More explicitly, for each queue j∈[n], the residual inter-arrival time  keeps track of the time remaining until the next arrival to queue j occurs. Immediately after an arrival to queue j occurs, the next inter-arrival time is drawn from a probability distribution F_j^A. When a job arrives to queue j, it comes with a workload (service requirement) drawn from a distribution F_j^S. We allow the distributions F_j^A's and F_j^S's to vary with time, i.e., the interarrival times and service requirements can be time-varying. For notational simplicity, we will not explicitly denote the time dependence here. We refer to the residual workload at time t_k of the top-of-queue job in queue j as , which specifies how much work must be done before the job completion. A job is only processed if it is routed to a server i∈[m], in which case the server processes the job at a constant service rate μ_ij∈_+. We refer to μ∈_+^m × n  as the matrix of service rates.
Under this scheduling decision, the residual processing time, i.e., the amount of time required
to process the job, is τ^S_k,j= /μ_ij. 






The augmented state s_k is a valid Markovian descriptor of the system and we now describe the corresponding  transition function f such that

    s_k+1 = f(s_k, u_k, ξ_k+1),

where u_k is an action taken by the controller and ξ_k+1 contains external randomness arising from new inter-arrival times or workloads drawn from F^A_j's or F^S_j's depending on the event type.

The transition is based on the next event, which is the event with the minimum residual time. The controller influences the transitions through the processing times, by deciding which jobs get routed to which servers. We focus on scheduling problems where the space of controls 𝒰 are feasible assignments of servers to queues. Let _n∈^n denote an n-dimensional vector consisting of all ones. The action space is, 

    𝒰:={ u∈{0,1}^m× n:u_n = _m, _m^⊤ u=_n,u≤ M},

where M∈{0,1}^m× n is the topology of the network, which indicates which job class can be served by which server. Following existing works on scheduling in queuing networks <cit.>, we consider networks for which each job class has exactly 1 compatible server.
 

For every queue j, there is 1 compatible server,
i.e., ∑_i=1^mM_ij=1.





Given an action u, the residual processing time is /μ_ij when u_ij=1 and ∞ 
when u_ij=0. This can be written compactly as 





    τ^S_k,j≡/∑_i=1^mu_ijμ_ij = /_j,j,

where 𝐝𝐢𝐚𝐠(u^⊤μ) ∈^n × n extracts the diagonal entries of the matrix u^⊤μ∈^n × n.

As a result, at time t_k the residual event times τ_k∈_+^2n consists of the residual inter-arrival and processing times,








    τ_k≡ (τ_k^A, τ_k^S) = (τ_k^A, ^-1 )

We emphasize that τ_k depends on the action u. The core operation in the transition dynamics is the event selection mechanism. The next event is the one with the minimum residual time in τ_k.

We define _k+1∈{0,1}^2n to be a one-hot vector representing the  of τ_k – the position of the minimum in τ_k:

    _k+1(_k, u_k)    ≡ (τ_k)  ∈{0,1}^2nEvent Select

_k+1(_k, u_k) indicates the type of the (k+1)th event.
In particular, if the minimum residual event time is a residual inter-arrival time, then the next event is an arrival to the system. If it is a residual job processing time, then the next event is a job completion.
We denote τ^*_k+1 to be the inter-event time, which is equal to the minimum residual time:

    τ^*_k+1(_k, u_k) 
    = min{τ_k}Event Time

τ^*_k+1(_k, u_k) is the time between the kth and (k+1)th event, i.e. t_k+1 - t_k.

After the job is processed by a server, it either leaves the system or proceeds
to another queue. Let R∈ℝ^n× n denote the routing matrix, where the jth column, R_j details the change in the queue lengths when a job in class j finishes service. For example, for a tandem queue with two queues, the routing
matrix is

    R=[[ -1  0;  1 -1 ]]

indicating that when a job in the first queue completes service, it leaves its
own queue and joins the second queue. When a job in the second queue
completes service, it leaves the system. 



We define the event matrix D as a  block matrix of the form

    D=[[ I_n   R ]],

where I_n is the n× n identity matrix. The event matrix determines the update to the queue lengths, depending on which event took place. In particular, when the (k+1)th event occurs, the update to the queue lengths is

    x_k+1   =x_k+D _k+1(_k, u_k)
    Queue Update

Intuitively, the queue length of queue j increases by 1 when the next event is a class j job arrival; the queue lengths update according to R_j when the next event is a queue j job completion.

The updates to the auxiliary state _k = (τ_k^A, w_k) ∈_+^2n is typically given by







    [
    [ τ_k+1^A;   w_k+1 ]]
       = 
    [
    [ τ_k^A;   w_k ]]
    - 
    τ_k+1^*[
    [ _n;    ]]
    _reduce residual times + 
    [
    [ T_k+1; W_k+1 ]] ⊙ e_k+1_draw new times / workloadsAux Update

where ⊙ is the element-wise product and
T_k+1 ={T_(k+1),j}_j=1^n∈^n are new inter-arrival times T_(k+1),j∼ F_j^A 
and W_k+1 = {W_k+1,j}_j=1^n∈^n are workloads W_k+1,j∼ F_j^S.

Intuitively, after an event occurs, we reduce the residual inter-arrival times by the inter-event time. We reduce workloads by the amount of work applied to the job, i.e., the inter-event time multiplied by the service rate of the allocated server. Finally, if an arrival occurred we draw a new inter-arrival time; if a job was completed, we draw a new workload for the top-of-queue job (if the queue is non-empty).
r0.4
  

  
    < g r a p h i c s >

  

  M/M/1 queue. 


There are two boundary cases that make the update slightly different from (<ref>). First, if a new job arrives at an empty queue j (either an external arrival or a transition from a job completion), we also need to update w_k+1,j to W_k+1,j. Second, if a queue j job completion leaves an empty queue behind, we set w_k+1, j = ∞, indicating that no completions can occur for an empty queue.




Let ξ_k denote exogenous noise in the environment, which consists of the sampled inter-arrival times and workloads for resetting the time of a completed event,

    ξ_k = (T_k, W_k) ∈_+^2n.

We finally arrive at the stated goal of describing the transition dynamics of s_k = (x_k, _k) in terms of a function f(s_k, u_k, ξ_k+1). Notably, all the stochasticity is captured by ξ_k's, which are independent of the states and actions. 






































































































It is worth mentioning a few features of this discrete-event representation.

    
  * While auxiliary data _k = (τ^A_k, w_k) is necessary for s_k = (x_k, _k) to be a valid Markovian system descriptor, this information is typically not available to the controller. We assume the controller only observes the queue lengths, i.e., the control policy π only depends on x_k.
    
  * The representation can flexibly accommodate non-stationary and non-exponential event-time distributions, i.e., F_j^A's and F_j^S's can be general and time-varying.
    
  * This model enables purely data-driven simulation, as it only requires samples of the event times ξ_k. One does not need to know the event time distributions F_j^A's and F_j^S's to simulate the system if data of these event times are available.
    
  * The matrix-vector representation enables GPU parallelism, which can greatly speed up the simulation of large-scale networks. 
    
    
  * As we will explain later, this representation enables new gradient estimation strategies.





  
Queuing Network Examples
As a concrete illustration, we show how a few well-known queuing networks are described as discrete-event dynamical systems.








The M/M/1 queue (see Figure <ref>) with arrival rate λ>0 and
service rate μ≥λ features a single queue n=1 and a single server m=1, and exponentially distributed inter-arrival times and workloads, i.e., T_k ∼𝖤𝗑𝗉(λ) and
W_k ∼𝖤𝗑𝗉(1) respectively. The network topology is M=[1], the service rate is μ, and the routing matrix is R =[-1], indicating that jobs leave the system after service completion.
The scheduling policy is work-conserving, the server always serves the queue when it is non-empty, i.e. u_k=1{x_k>0}. The state update is,

    x_k+1   =x_k+[[  1 -1 ]]^⊤_k+1
    _k+1   =min{τ_k^A,τ_k^S} =min{τ_k^A,/μ· 1{x_k>0 }}∈{0,1}^2
    τ^*_k+1   = min{τ_k^A,τ_k^S}
    [
    [ τ_k+1^A;   w_k+1 ]]
       = 
    [
    [ τ_k^A;   w_k ]]
    - τ_k+1^*[
    [             1; μ 1{ x_k > 0} ]] +
    [
    [ T_k+1; W_k+1 ]] ⊙ e_k+1.









The multi-class singer-server queue features an n queues and a single server m=1 (see Figure <ref>). While the inter-arrival times and workloads, i.e., (T_k, W_k)'s are usually exponentially distributed, they can also follow other distributions. The network topology is M=[1,...,1] ∈^n, the service rates are μ = [μ_1,...,μ_n], and the routing matrix is R =[-1,...,-1]∈^n, indicating that jobs leave the system after service completion. 
A well-known scheduling policy for this system is the cμ-rule, a static priority rule.

Let h = (h_1,...,h_n) ∈^n denote the holding costs. The cμ-rule sets

    u_k=_j∈ [n]{ h_jμ_j1{x_j>0}}∈{0,1}^n.


The state update is,

    x_k+1   =x_k+[[  _n -_n ]]^⊤_k+1
    _k+1   =min{τ_k,1^A,...,τ_k,n^A,τ_k,1^S,...,τ_k,n^S} 
    =min{τ_k,1^A,...,τ_k,n^A,
    w_k,1/μ_1 u_k,1,...
    w_k,n/μ_n u_k,n}
    τ^*_k+1   = min{τ_k,1^A,...,τ_k,n^A,τ_k,1^S,...,τ_k,n^S}
    [
    [ τ_k+1^A;   w_k+1 ]]
       = 
    [
    [ τ_k^A;   w_k ]]
    - τ_k+1^*[
    [     _n; μ⊙ u_k ]] +
    [
    [ T_k+1; W_k+1 ]] ⊙ e_k+1.







The criss-cross network <cit.> features n=3 queues and m=2 servers (see Figure <ref>). External jobs arrive to queues 1 and 3. The first server can serve queues 1 and 3 with service rates μ_11 and μ_13 respectively, while the second server is dedicated to serving queue 2 with service rate μ_22. After jobs from queue 1 are processed, they are routed to queue 2; jobs from queues 2 and 3 exit the system after service completion. The inter-arrival times and workloads, i.e., (T_k, W_k)'s, can follow general distributions.
The network topology M, service rate matrix, and the routing matrix R are:

    M = [ [ 1 0 1; 0 1 0;   ]],     μ = [ [ μ_11    0 μ_13;    0 μ_22    0;      ]],     
    R = [ [ -1  0  0;  1 -1  0;  0  0 -1 ]]


<cit.> develop a work-conserving threshold policy for this system. For a threshold a ∈_+, server 1 prioritizes jobs in queue 1 if the number of jobs in queue 2 is below a. Otherwise, it prioritizes queue 3. This gives the scheduling action

    u_k,11 = 1{x_k,2≤ a },   
    u_k,22 = 1{x_k,2 > 0 },    
    u_k,13 = (1 - u_k,11)1{x_k,3 > 0},

and the transition dynamics 

    x_k+1   =x_k+[[ I_3   R ]]
    _k+1
    _k+1   =min{τ_k,1^A, ∞, τ_k,3^A, 
    τ_k,1^S,
    τ_k,2^S,
    τ_k,3^S}  
    =min{τ_k,1^A, ∞, τ_k,3^A, 
    w_k,1/μ_11u_k,11,
    w_k,2/μ_22u_k,22, 
    w_k,3/μ_13u_k,13}
    τ^*_k+1   = min{τ_k,1^A, ∞, τ_k,3^A, 
    τ_k,1^S,
    τ_k,2^S,
    τ_k,3^S}
    [
    [ τ_k+1^A;   w_k+1 ]]
       = 
    [
    [ τ_k^A;   w_k ]]
    - τ_k+1^*[
    [ _3;    ]] +
    [
    [ T_k+1; W_k+1 ]] ⊙ e_k+1.

Here, τ_k,2^A = ∞ since queue 2 has no external arrivals. 

































































§ GRADIENT ESTIMATION





In this section, we introduce our proposed approach for estimating the gradient of the objective (<ref>), ∇ J_N(θ). We start with a brief discussion of existing methods for gradient estimation, including their advantages and limitations. We then outline the main challenges for computing pathwise derivatives in multi-class queuing networks, and introduce our strategy for overcoming these challenges. Finally, we formally define our gradient estimation framework and discuss its computational and statistical properties. Later in section <ref>, we perform a comprehensive empirical study and find that our gradient estimation framework is able to overcome many of the limitations of existing methods in that (1) it is capable of estimating gradients for general queuing networks, (2) it provides stable gradient estimations over very long horizons (>10^5 steps), (3) it provides greater estimation accuracy than model-free policy gradient methods with 1000x less data, 
and (4) when applying to policy optimization, it drastically improves the performance of the policy gradient algorithm for various scheduling and admission control tasks. 

L0.3
  

  
    < g r a p h i c s >

  

  Multi-class, single-server queue. 


Our goal is to optimize the parameterized control problem (<ref>). A standard optimization algorithm is (stochastic) gradient descent, which has been considered for policy optimization and reinforcement learning <cit.>. The core challenge for estimating policy gradient ∇ J_N(θ) = ∇[J(θ;ξ_1:N)] from sample paths of the queuing network is that the sample path cost J(θ, ξ_1:N) is in general not differentiable in θ. As a consequence, one cannot change the order of differentiation and expectation, i.e.,

    ∇ J_N(θ)=∇[J_N(θ;ξ_1:N)]𝔼[∇ J_N(θ;ξ_1:N)],

where ∇ J_N(θ;ξ_1:N) is not even well-defined.
The non-differentiability of these discrete-event dynamical systems emerges from two sources. First, actions u_1:N are discrete scheduling decisions, and small perturbations in the policy can result in large changes in the scheduling decisions produced by the policy. Second, the actions affect the dynamics through the event times. The ordering of events is based on the `𝖺𝗋𝗀𝗆𝗂𝗇' of the residual event times, which is not differentiable.



In the stochastic simulation literature, there are two popular methods for gradient estimation: infinitesimal perturbation analysis (IPA) and generalized likelihood ratio (LR) gradient estimation.
To illustrate, consider abstractly and with a little abuse of notation a system following the dynamics s_k+1 = f(s_k, θ, ξ_k+1), where s_k ∈ is the state, θ∈ is the parameter of interest,  ξ_k is exogenous stochastic noise, and f is a differentiable function. Then, the IPA estimator computes a sample-path derivative estimator
by constructing a derivative process D_k=∂ s_k / ∂θ via the recursion:

    D_k+1 = ∂/∂θ f(s_k, θ, ξ_k+1) + ∂/∂ s_k f(s_k, θ, ξ_k+1) · D_kIPA

Likelihood-ratio gradient estimation on the other hand uses knowledge of the distribution of ξ_k to form the gradient estimator. Suppose that s_k is a Markov chain for which the transition kernel is parameterized by θ, i.e., s_k+1∼ p_θ(· |s_k). 
For a fixed θ_0, let 

    ∂/∂θ_θ[s_k] = ∂/∂θ_θ_0[s_kL_k(θ)] = _θ_0[s_k∂/∂θL_k(θ)]
        L_k(θ) := ∏_j=1^k-1 p_θ(s_j+1|s_j)/∏_j=1^k-1 p_θ_0(s_j+1|s_j).

This allows one to obtain the following gradient estimator:

    D_k = s_k∑_j=1^k-1∂/∂θp_θ(s_j+1|s_j)/p_θ_0(s_j+1|s_j)L_k(θ),  s_j+1∼ p_θ_0(·|s_j),∀ j ≤ k.
    LR


Despite their popularity, there are limitations to applying these methods to general multi-class queuing networks. While IPA has been proven efficient for simple queuing models, such as the G/G/1 queue through the Lindley recursion, it is well-known that unbiased IPA estimates cannot be obtained for general queuing networks <cit.>. The implementation of LR gradient estimation hinges on precise knowledge of the system's Markovian transition kernel  <cit.>. This requires knowledge of the inter-arrival time and workload distributions, and even with this knowledge, it is non-trivial to specify the transition kernel of the queue lengths and residual event times in generic systems. Modifications to IPA <cit.> also require precise knowledge of event time distributions and often involve analyzing specific ordering of events which must be done on a case-by-case basis. As a result, none of these methods can reliably provide gradient estimation for complex queuing networks under general scheduling policies and with possibly unknown inter-arrival and service time distributions. Yet, the ability to handle such instances is important to solve large-scale problems arising in many applications.


Due to the challenges discussed above, existing reinforcement learning (RL) approaches for queueing network control mainly rely on model-free gradient estimators, utilizing either the  estimator and/or Q-function estimation. As we will discuss shortly, these methods do not leverage the structural properties of queuing networks and may be highly sample-inefficient, e.g., requiring a prohibitively large sample for gradient estimation.

To address the challenges discussed above, we propose a novel gradient estimation framework that can handle general, large-scale multi-class queuing networks under any differentiable scheduling policy, requiring only samples of the event times rather than knowledge of their distributions. Most importantly, our approach streamlines the process of gradient estimation, leveraging auto-differentiation libraries such as PyTorch <cit.> or Jax <cit.> to automatically compute gradients, rather than constructing these gradients in a bespoke manner for each network as is required for IPA or LR. As shown in Figure <ref>, computing a gradient in our framework requires only a few lines of code. To the best of our knowledge, this is the first scalable alternative to model-free methods for gradient estimation in queuing networks.





 



 §.§ The standard approach: the  estimator







Considering the lack of differentiability in most reinforcement learning environments, the standard approach for gradient estimation developed in model-free RL is the score-function or  estimator <cit.>. This serves as the basis for modern policy gradient algorithms such as Trust-Region Policy Optimization (TRPO) <cit.> or Proximal Policy Optimization (PPO) <cit.>. As a result, it offers a useful and popular baseline to compare our proposed method with.

 The core idea behind the  estimator is to introduce a randomized policy π_θ and differentiate through the action probabilities induced by the policy. Under mild regularity conditions on π_θ and c(x_k,u_k), the following expression holds for the policy gradient:

    ∇ J_N(θ) = 𝔼[∑_t=0^N-1(∑_k=t^N-1c(x_k,u_k)τ^*_k+1) ∇_θlogπ_θ(u_t|x_t)],


which leads to the following policy gradient estimator:

    ∇^𝖱 J_N(θ; ξ_1:N) = ∑_t=0^N-1(∑_k=t^N-1c(x_k,u_k)τ^*_k+1) ∇_θlogπ_θ(u_t|x_t).


While being unbiased, the 𝖱𝖤𝖨𝖭𝖥𝖮𝖱𝖢𝖤 estimator is known to have a very high variance <cit.>. The variance arises from two sources. First, the cumulative cost ∑_k=t^N-1c(x_k,u_k)τ^*_k+1 can be very noisy, as has been observed for queuing networks <cit.>. Second, 

as the policy converges to the optimal policy, the score function ∇_θlogπ_θ(u_t|x_t) can grow large, magnifying the variance in the cost term.  Practical implementations involve many algorithmic add-ons to reduce variance, e.g., adding a `baseline' term <cit.> which is usually (an estimate of) the value function V_π_θ(x_k),

    ∇^𝖱𝖡 J_N(θ; ξ_1:N) = ∑_t=0^N-1(∑_k=t^N-1c(x_k,u_k)τ^*_k+1  
     - V_π_θ(x_k)) ∇_θlogπ_θ(u_t|x_t).
     𝖡𝖠𝖲𝖤𝖫𝖨𝖭𝖤

These algorithmic add-ons have led to the increased complexity of existing policy gradient implementations <cit.> and the outsized importance of various hyperparameters <cit.>. It has even been observed that seemingly small implementation “tricks" can have a large impact on performance, even more so than the choice of the algorithm itself <cit.>.







 §.§  Our approach: Differentiable Discrete-Event Simulation









We can view the state trajectory as a repeated composition of the transition function s_k+1 = f(s_k,u_k, ξ_k+1), which is affected by exogenous noise ξ_1:N, i.e., stochastic inter-arrival and service times. If the transition function were differentiable with respect to the actions u_k, then under any fixed trace ξ_1:N, one could compute a sample-path derivative of the cost J(θ;ξ_1:N) using auto-differentiation frameworks such as PyTorch <cit.> or Jax <cit.>. Auto-differentiation software computes gradients efficiently using the chain rule. To illustrate, given a sample path of states, actions, and noise (s_k,u_k,ξ_k+1)_k=0^N-1, we can calculate the gradient of s_3 with respect to u_1 via

    ∂ s_3/∂ u_1 = ∂ s_3/∂ s_2∂ s_2/∂ u_1 =
    ∂ f(s_2,u_2, ξ_3)/∂ s_2∂ f(s_1,u_1, ξ_2)/∂ u_1.

This computation is streamlined through a technique known as backpropagation, or reverse-mode auto-differentiation. The algorithm involves two steps. The first step, known as the forward pass, evaluates the main function or performance metric (in the example, s_i's) and records the partial derivatives of all intermediate states relative to their inputs (e.g. ∂ s_2 / ∂ u_1). 
This step constructs a computational graph, which outlines the dependencies among variables. The second step is a backward pass, which traverses the computational graph in reverse. It sequentially multiplies and accumulates partial derivatives using the chain rule, propagating these derivatives backward through the graph until the gradient concerning the initial input (in this example, u_1) is calculated. Due to this design, gradients of functions involving nested compositions can be computed in a time that is linear in the number of compositions. By systematically applying the chain rule in reverse, auto-differentiation avoids the redundancy and computational overhead typically associated with numeric differentiation methods.

However, as mentioned before, the dynamics do not have a meaningful derivative due to the non-differentiability of actions and the 𝖺𝗋𝗀𝗆𝗂𝗇 operation which selects the next event based on the minimum residual event time. Yet if we can utilize suitably differentiable surrogates, it would be possible to compute meaningful approximate sample-path derivatives using auto-differentiation.







  §.§.§ Capacity sharing relaxation
 


First, we address the non-differentiability of the action space. Recall that u_k∈{0,1}^m × n are scheduling decisions, which assign jobs to servers. Since u_k lies in a discrete space, a small change in the policy parameters can produce a jump in the actions. To alleviate this, we consider the transportation polytope as a continuous relaxation of the original action space (<ref>):






    𝒰:={ u∈ [0,1]^m× n:u_n = _m, _m^⊤ u=_n,u≤ M}.

The set of extreme points of 𝒰 coincide with the original, integral action space 𝒰. 
For a fractional action u_k∈𝒰, 
we can interpret it as servers splitting their capacity among multiple job classes motivated by the fluid approximation of queues <cit.>. As a relaxation, it allows servers to serve multiple jobs simultaneously. The effective service rate for each job class is equal to the fraction of the capacity allocated to the job class multiplied by the corresponding service rate. 


As a result, instead of considering stochastic policies over discrete actions, we approach this problem as a continuous control problem and consider
deterministic policies over continuous actions, i.e., the fractional scheduling decisions. Under this relaxation, the processing times are differentiable in the (fractional) scheduling decision. Finally, it is worth mentioning that we only use this relaxation when training policies. For policy evaluation, we enforce that actions are integral scheduling decisions in 𝒰. To do so, we treat the fractional action as a probability distribution and use it to sample a discrete action. 





Under the capacity sharing relaxation, the service rate for queue j under the routing decision u∈𝒰 is μ_j^⊤ u_j≡∑_i=1^mμ_iju_ij.
Thus, given workload w_j, the processing time of the job will be

    τ^S_j =w_j/∑_i=1^mu_ijμ_ij = w_j/_j,j.


Note that this is identical to the original definition of the processing times in (<ref>). The only difference is that we now allow fractional routing actions, under which a server can serve multiple jobs at the same time.


For a concrete example, consider a single server i compatible with two job classes 1 and 2 with service rates μ_i1 = 9 and μ_i2 = 15 respectively. Suppose it splits its capacity between job classes 1 and 2 according to u_i1 = 1/3 and u_i2 = 2/3. Then for residual workloads w_1 and w_2, the corresponding processing times are τ^S_1=w_1/3 and τ^S_2=w_2/10. If u_i1=0 and u_i2=1 instead, then the corresponding processing times are τ^S_1=w_1/0 = ∞ and τ^S_2=w_2/15. 









  §.§.§ Differentiable event selection



To determine the next event type, the 𝖺𝗋𝗀𝗆𝗂𝗇 operation selects the next event based on the minimum residual event time. This operation does not give a meaningful gradient.



  
Pitfalls of `naive' smoothing
In order to compute gradients of the sample path, we need to smooth the 𝖺𝗋𝗀𝗆𝗂𝗇 operation. There are multiple ways to do this. A naive approach is to directly replace 𝖺𝗋𝗀𝗆𝗂𝗇 with a differentiable surrogate. 
One such popular surrogate is . With some inverse temperature β > 0, _β applied to the vector of residual event times τ∈_+^2n returns a vector in _+^2n, 
which we use to replace the event selection operation _k+1: 

    _k+1 = 𝗌𝗈𝖿𝗍𝗆𝗂𝗇_β(τ_k)
          𝗌𝗈𝖿𝗍𝗆𝗂𝗇_β(τ)_j=e^-βτ_j / ∑_l=1^2n e^-βτ_l.
    Direct Smoothing

As β→∞, 𝗌𝗈𝖿𝗍𝗆𝗂𝗇_β converges to 𝖺𝗋𝗀𝗆𝗂𝗇. Thus, one may expect that for large β, 𝗌𝗈𝖿𝗍𝗆𝗂𝗇_β would give a reliable differentiable surrogate. 

However, queuing networks involve a unique challenge for this approach: one typically considers very long trajectories when evaluating performance in queuing networks, as one is often interested in long-run average or steady-state behavior. Thus, even if one sets β to be very large to closely approximate , the smoothing nonetheless results in `unphysical', real-valued queue lengths instead of integral ones, and small discrepancies can accumulate over these long horizons and lead to entirely different sample paths. 
This can be observed concretely in the left panel of Figure <ref>, which displays the sample paths of the total queueing length processes for a criss-cross queueing network (in Example <ref>) under the original dynamic and under direct smoothing, using the same inter-arrival and service times.  We observe that when setting the inverse temperature β = 1, the sample path under direct smoothing is completely different from the original one, even though all of the stochastic inputs are the same. 

Even when setting a very high inverse temperature, i.e., β = 1000, for which 𝗌𝗈𝖿𝗍𝗆𝗂𝗇_β is almost identical to 𝖺𝗋𝗀𝗆𝗂𝗇, the trajectory veers off after only a hundred steps. 




This can greatly affect the quality of the gradient estimation. 

We observe in the right panel of Figure <ref> that across a range of inverse temperatures, the average cosine similarity between the surrogate gradient and the true gradient (defined in (<ref>)) are all somewhat low.

In the same plot, we also show the average cosine similarity between our proposed gradient estimator, which we will discuss shortly, and the true gradient. Our proposed approach substantially improves the gradient estimation accuracy, i.e., the average cosine similarity is close to 1, and as we will show later, it does so across a wide range of inverse temperatures. 






  
Our approach: `straight-through' estimation

The failure of the direct smoothing approach highlights the importance of preserving the original dynamics, as errors can quickly build up even if the differentiable surrogate is only slightly off. We propose a simple but crucial adjustment to the direct smoothing approach, which leads to huge improvements in the quality of gradient estimation. 


Instead of replacing the 𝖺𝗋𝗀𝗆𝗂𝗇 operation with 𝗌𝗈𝖿𝗍𝗆𝗂𝗇_β when generating the sample path, we preserve the original dynamics as is, and only replace the Jacobian of 𝖺𝗋𝗀𝗆𝗂𝗇 with the Jacobian of _β when we query gradients. In short, we introduce a gradient operator ∇ such that

    _k+1=𝖺𝗋𝗀𝗆𝗂𝗇(τ_k),     ∇_k+1= ∇_β(τ_k).

where ∇ is respect to the input τ.









This is known as the `straight-through' trick in the machine learning literature and is a standard approach for computing approximate gradients in discrete environments <cit.>. To the best of our knowledge, this is the first application of this gradient estimation strategy for discrete-event dynamical systems. Using this strategy, we can use the chain rule to compute gradients of performance metrics that depend on the event selection. 

Consider any differentiable function g of e_k+1, 


    ∇ g(_k+1) = ∂ g(_k+1)/∂_k+1∇_k+1∇τ_k= ∂ g(_k+1)/∂_k+1∇_β(τ_k) ∇τ_k

In contrast, direct smoothing involves the derivative ∂ g(ẽ_k+1)/∂ẽ_k+1 where ẽ_k+1 = _β(τ_k).  Evaluating the gradient of g at g(ẽ_k+1) is a cause of additional bias.



With these relaxations, the transition function of the system is differentiable. We can now compute a gradient of the sample path cost J_N(θ; ξ_1:N), using the chain rule on the transition functions. Given a sample path of states s_k = (x_k, _k), actions u_k = π_θ(x_k), and ξ_k = (T_k, W_k), the pathwise gradient of the sample path cost is J_N(θ; ξ_1:N) with respect to an action u_k is,

    ∇_u_k J_N(θ; ξ_1:N)
    = ∇_u_kc(x_k,u_k)_current cost + ∑_t=k+1^N[  ∇_xc(x_t, u_t) + ∇_u c(x_t, u_t)∇_x_tπ_θ(x_t)_future costs]
    ∇_u_kx_t

The gradient consists of the sensitivity of the current cost with respect to the action as well the sensitivity of future costs via the current action's impact on future states. The policy gradient with respect to θ can then be computed as

    ∇_θ J_N(θ; ξ_1:N) = ∑_k=1^N. ∇_u_k J_N(θ; ξ_1:N) |_u_k = π_θ(x_k)∇_θπ_θ(x_k).

As a result of the straight-through trick, we do not alter the event selection operation _k's and thus the state trajectory {x_k}_k=1^N is unchanged.



We refer to the gradient estimator (<ref>) as the 𝖯𝖠𝖳𝖧𝖶𝖨𝖲𝖤 policy gradient estimator. Although this formula involves iterated products of several gradient expressions, these can be computed efficiently through reverse-mode auto-differentiation using libraries such as PyTorch <cit.> or Jax <cit.> with O(N) time complexity in the time horizon N. This time complexity is of the same order as the forward pass, i.e., generating the sample path itself, and is equivalent to the time complexity of 𝖱𝖤𝖨𝖭𝖥𝖮𝖱𝖢𝖤. The policy gradient algorithm with 𝖯𝖠𝖳𝖧𝖶𝖨𝖲𝖤 gradient is summarized in Algorithm <ref>. In Section <ref>, we perform a careful empirical comparison of  and , and find that  can lead to orders of magnitude improvements in sample efficiency. 













It is important to re-emphasize that when computing the gradient, we evaluate the policy differently than we would for . Instead of drawing a random discrete action u_k∼π_θ(x), we use the probabilities output by the policy directly as a fractional routing matrix in 𝒰,

    :      u_k∼π_θ(x_k),     u_k∈𝒰
    :      u_k = π_θ(x_k),     u_k∈𝒰.

However, this is only for gradient computation. When we evaluate the policy, we draw u_k∼π_θ(x_k).





  
Bias-variance trade-off for inverse temperature: One-step analysis
The inverse temperature β is a key hyperparameter that determines the fidelity of the _β approximation to . The choice of β poses a bias-variance trade-off, with a higher β leading to a smaller bias but a higher variance and a smaller β incurring a higher bias but a lower variance. 

In general, it is difficult to assess the bias since we  often do not know the true gradient. However, for some simple examples, we can evaluate the true gradient explicitly. We next analyze the gradient of the one-step transition of the M/M/1 queue with respect to the service rate μ,

    ∇_μ[x_k+1 - x_k]=  ∇_μ[D e_k+1] =∇_μλ - μ/λ + μ = -2λ/(λ + μ)^2.

This permits an exact calculation of the mean and variance of our proposed pathwise gradient estimator. Although we can derive analytical expressions for these quantities, we present the leading order asymptotics as β→∞ for conciseness of presentation. While it is straightforward to see that almost surely _β(τ)→(τ) as β→∞, it is much less clear whether the gradient converges, i.e., whether [∇_β(τ)] →∇[(τ)]. Since  has a gradient of zero almost everywhere, the expectation and gradient operators cannot be interchanged. 

Instead, we analyze the expectations directly using properties of the exponential distribution.


Let ∇_μ (x_k+1 - x_k) =  ∇_μ De_k+1 denote the 𝖯𝖠𝖳𝖧𝖶𝖨𝖲𝖤 gradient estimator of the one-step transition of the M/M/1 queue with respect to μ.  For x_k≥ 1, as β→∞,

    [∇_μ De_k+1]
    -∇_μ[ De_k+1]
       = β^-2·π^2λ (μ^2 - λ^2 + 2μλ)/6(λ + μ)^2
    + o ( β^-2), 
    (∇_μ De_k+1)
       = β·4λ/μ (λ + μ)^2 + o(β).


See section <ref> for the proof.
As β→∞, the bias is O(1/β^2) while the variance is O(β). This means that one can significantly reduce the bias with only a moderate size of β. The left panel of Figure <ref> shows the bias-variance trade-off of ∇_μ (x_k+1 - x_k) for various inverse temperatures β. The blue line is based on the analytical expression for the bias-variance trade-off curve. We observe that with the inverse temperature β∈ [0.5, 2], both the bias and the variance are reasonably small.







Suppose we compute the sample average of B iid samples of ∇_μ De_k+1, which are denoted as ∇_μ De_k+1, i, i=1,…, K. In particular, the estimator takes the form 1/B∑_i=1^B∇_μ De_k+1,i. The choice of β that minimizes the mean-squared error (MSE) of the estimator 

is β^* = O(B^1/5) and MSE(β^*)=O(B^-4/5).




The  estimator provides a more statistically efficient trade-off than other alternatives. As an example, a standard gradient estimator is the finite-difference estimator in which one evaluates the one-step transition at μ - h and μ + h for some small h ∈ (0,∞), and the estimator is constructed as

    1/B∑_i=1^BDe_k+1,i(μ + h) - De_k+1,i(μ - h)/2h,

where De_k+1,i(μ + h)'s are iid samples of De_k+1(μ + h).
If we set h = 1/β, it is well-known that the bias scales as O(1/β^2) while the variance scales as O(β^2). The choice of β that minimizes the MSE is β^* = O(B^1/6) and MSE(β^*)=O(B^-1/3).


While this analysis is restricted to the one-step transition of the M/M/1 queue, these insights hold for more general systems and control problems. The right panel of Figure <ref> displays the average cosine similarity (defined in  (<ref>)) between the  gradient estimator and the true gradient for a policy gradient task in a 6-class reentrant network across different congestion levels and for different inverse temperatures. We observe that for a wide range of inverse temperatures, β∈{0.2, 0.5, 1, 2}, the estimator has near-perfect similarity with the true gradient, while a very large inverse temperature suffers due to high variance. This indicates that while there is a bias-variance trade-off, the performance of the  gradient estimator is not sensitive to the choice of the inverse temperature within a reasonable range. In our numerical experiments, we 
find that one can get good performance using the same inverse temperature across different settings without the need to tune it for each setting. 










§ EMPIRICAL EVALUATION OF THE  GRADIENTS



In the previous section, we introduced the  gradient estimator for computing gradients of queuing performance metrics with respect to routing actions or routing policy parameters. In this section, we study the statistical properties of these gradient estimators and their efficacy in downstream policy optimization tasks. We use  as the baseline gradient estimator. First, in section <ref> we empirically study the estimation quality across a range of queuing networks, traffic intensities, and policies. After that, in section <ref>, we investigate their performance in a scheduling task: learning the cμ rule in a multi-class queuing network. Finally, we demonstrate the applicability of our framework beyond scheduling: we investigate the performance of the  gradient estimator for admission control tasks in section <ref>. 





 §.§ Gradient Estimation Efficiency


In general, it is challenging to theoretically compare the statistical properties of different gradient estimators, and very few results exist for systems beyond the M/M/1 queue (see section <ref> for a theoretical comparison between  and  for the M/M/1 queue). For this reason, we focus on numerical experiments across a range of environments and queuing policies typically considered in the queuing literature. Specifically, we will be comparing the statistical properties of 
 estimator with the baseline estimator . While  introduces bias into the estimation, we find in our experiments that this bias is small in practice and remains small even over long time horizons. At the same time, the  estimator delivers dramatic reductions in variance, achieving greater accuracy with a single trajectory than  with 10^3 trajectories.





First, recall that a policy π(x) maps queue-lengths x to assignment between servers and queues, represented by an m × n matrix in 𝒰 (allowing for fractional routing matrices). We visit three classical queuing policies: priority policies <cit.>, MaxWeight <cit.>, and MaxPressure <cit.>. Each of these methods selects the routing that solves an optimization problem. This means that the routing generated by the policy is deterministic given the state and is not differentiable in the policy parameters. In order to apply either  or the  gradient estimator to compute a policy gradient, we require differentiable surrogates of these policies. To this end, we define softened and parameterized variants of these policies, denoted as soft priority (𝗌𝖯𝖱), soft MaxWeight (𝗌𝖬𝖶), and soft MaxPressure (𝗌𝖬𝖯),

    π^𝗌𝖯𝖱_θ(x)_i= (θ_j·μ_i),   π^𝗌𝖬𝖶_θ(x)_i= 
    (θ_jx_j·μ_i),   π^𝗌𝖬𝖯_θ(x)_i = ((μ⊙ R(θ x))_i)

where θ∈_+^n are a vector of costs/weights for each queue, μ denotes the matrix of service rates with μ_i∈_+^n denoting the service rates associated with server i. 
The operation ⊙ refers to element-wise multiplication and the
 operation maps a vector a ∈^n into a set of probabilities (a)_i = e^a_i / ∑_j=1^n e^a_j.






We are interested in identifying the parameter θ that minimizes long-run average holding cost where c(x,u) = h^⊤ x. 

We use the objective 
J_N(θ) = [∑_k=0^N-1c(x_k,π_θ(x_k))τ^*_k+1]
where N is a large enough number to approximate the long-run performance, and the goal of the gradient estimation is to estimate ∇ J_N(θ).

We consider the following environments, which appear throughout our computational experiments and serve as standard benchmarks for control policies in multi-class queuing networks. We describe the network structure in detail in Figure <ref>.


  * Criss-cross: The network introduced in Example <ref> (see Figure <ref> (c)).

  * Re-entrant 1 (n classes): We consider a family of multi-class re-entrant networks with a varying number of classes, which was studied in <cit.>. The network is composed of several layers and each layer has 3 queues. Jobs processed in one layer are sent to the next layer. Arrivals to the system come to queues 1 and 3 in the first layer while queue 2 receives re-entered jobs from the last layer (see Figure <ref> (a) for a two-layer example).

  * Re-entrant 2 (n classes): We consider another family of re-entrant network architecture that was studied in <cit.>. It also consists of multiple layers with 3 queues in each layer. It differs from the Re-entrant 1 environment in that only queue 1 receive external arrivals while queues 2 and 3 receive re-entered jobs from the last layer (see Figure <ref> (b) for a two-layer example).


For a gradient estimator ĝ, the main performance metric we evaluate is (ĝ), which is the expected cosine similarity with the ground-truth gradient,

    (ĝ) ≡[𝖼𝗈𝗌( ĝ, ∇ J_N(θ) )] = [ ⟨ĝ, ∇ J_N(θ) ⟩/ĝ∇ J_N(θ)] ∈ [-1,1]

where the expectation  is over randomness in ĝ. The higher the similarity is, the more aligned ĝ is to the direction of ∇ J_N(θ). This metric incorporates both bias and variance of the gradient estimator. If the gradient estimator is unbiased but has a high variance, then each individual realization of ĝ is likely to have low correlation with the true gradient, so the average cosine similarity will be small even if [ĝ]=∇ J_N(θ). At the same time, if the gradient estimator has a low variance but a high bias, then the (ĝ) could still be small if 𝖼𝗈𝗌([ĝ], ∇ J_N(θ) ) is small. We focus on this metric, because it directly determines how informative the gradient estimates are when applying various gradient descent algorithms.  
For our experiments, we evaluate (a close approximation of) the ground-truth gradient ∇ J_N(θ)
by using the unbiased  gradient estimator over exceedingly many trajectories (in our case, 10^6 trajectories). 

We compare the similarity of  with that of . 

We denote B as the number of trajectories we use to calculate each  or  gradient estimator.

    ∇̂_θ J_N(θ; ξ^(1)_1:N)_ with B=1    ∇̂^𝖱_θ J_N,B(θ; ξ_1:N) := 1/B∑_b = 1^B∇̂^𝖱_θ J_N(θ; ξ^(b)_1:N)_ with B trajectories

We compute the  gradient with only B = 1 trajectory, while  gradient is calculated using B=10^3 trajectories.

For each policy and setting, we compute these gradients for 100 different randomly generated values of θ, which are drawn from a 𝖫𝗈𝗀𝗇𝗈𝗋𝗆𝖺𝗅(0,1) distribution (as the parameters must be positive in these policies). In total, we compare the gradients in 10,080 unique parameter settings, and each gradient estimator is computed 100 times to evaluate the average cosine similarity. When computing the policy gradient, we consider a time horizon of N = 10^3 steps. 





Figure <ref> compares the  estimator with B = 1 trajectory with the  estimator averaged over B = 10^3 trajectories. For the  estimator, costs are computed with a discount factor γ = 0.999, as using a lower discount rate introduced significant bias in the estimation. For , we use an inverse temperature β = 1 for the  relaxation across all settings. 

Each cell in Figure <ref> corresponds to a (policy, network, traffic-intensity) and the cell value is the average expected cosine similarity of the estimator averaged across the 100 randomly drawn θ values. We observe that across these diverse settings, the  estimator consistently has a much higher average cosine similarity with the true gradient despite using only a single trajectory. In fact, for 94.5% of the 10,800 parameter settings,  has a higher average cosine similarity with 99% confidence than  with B = 1000 trajectories. In most cases, the cosine similarity of  is close to 1, indicating almost perfect alignment with the true gradient even under high congestion.  on the other hand suffers greatly from high variance. 

Overall, this demonstrates that  is able to deliver greater estimation accuracy with an order of magnitude fewer samples.






    












 §.§ Learning the cμ rule


Given the strong improvements in estimation efficiency, we turn to evaluate how these translate to a downstream optimization task.
In single-server multi-class queues, it is well-known that the cμ-rule minimizes the long-run average holding cost <cit.>. We assess whether gradient descent with the  or  gradients is capable of converging to the cμ-rule, without knowing the holding costs h or μ and only using feedback from the environment. Despite its simplicity, it has been observed in prior work that this is a difficult learning task, particularly under heavy traffic <cit.>.

We revisit the soft priority policy mentioned before, but with only the parameters θ∈^n, i.e., 

    π^𝗌𝖯𝖱_θ(x)_i= (θ_i)

We also modify the policy to ensure that it is work-conserving, i.e., not assigning the server to an empty queue (see section <ref> for further discussion).

We consider a family of multi-class single-server queues with n queues. Holding costs are identically h_j = 1. Inter-arrival and service times are exponentially distributed, the service rates are μ_1j = 1 + ϵ j, for some ϵ>0, and the arrival rates are identical λ_j = λ and λ are set such that the traffic intensity ∑_j=1^nλ/μ_1j = ρ

for some ρ∈ (0,1). Note that in this case, the cμ-rule prioritizes queues with higher indices j. We consider a grid of gap sizes ϵ∈{1.0, 0.5, 0.1, 0.05, 0.01 } to adjust the difficulty of the problem; the smaller ϵ is, the harder it is to learn.

We compare  with B=1 trajectory and  with B = 100 trajectories for trajectories of N = 1000 steps. In order to isolate the effect of the gradient estimator from the optimization scheme, for both estimators we use an identical stochastic gradient descent scheme with normalized gradients (as these two estimators may differ by a scale factor). That is, for gradient estimator ĝ, the update under step-size α is

    θ_t+1 = θ_t - αĝ/ĝ.

We run T gradient descent steps for each gradient estimator. To allow for the fact that different estimators may have different performances across different step sizes, we consider a grid of step sizes α∈{0.01, 0.1, 0.5, 1.0 }. Gradient normalization may prevent convergence, so we use the averaged iterate θ̅_T for T. We then evaluate the long-run average holding cost under a strict priority policy determined by θ̅_T, i.e., π_θ̅_T(x)_i = θ̅_T,j. 







    








The left panel of Figure <ref> displays the values of θ̅_T after T=50 gradient iterates for  and  with n=5, ϵ = 0.1, and ρ = 0.99. We observe that while  sorts the queues in the correct order (it should be increasing with the queue index),  even with B = 100 trajectories fails to prioritize queues with a higher c μ index. Remarkably, we observe in the right panel of the same figure that  with just a single trajectory achieves a lower average holding cost than  uniformly across various step sizes and difficulty levels, whereas the performance of  varies greatly depending on the step size. This indicates that the improvements in gradient efficiency/accuracy of  make it more robust to the step-size hyper-parameter. 
It is also worth mentioning that when gap size ϵ becomes smaller, it is more difficult to learn. At the same time, since μ_1j's are more similar to each other, the cost difference between different priority rules also diminishes.







 §.§ Admission Control



While we focus mainly on scheduling tasks in this work, our gradient estimation framework can also be applied to admission control, which is another fundamental queuing control task <cit.>. To manage congestion, the queuing network may reject new arrivals to the network if the queue lengths are above certain thresholds. The admission or buffer control problem is to select these thresholds to balance the trade-off between managing congestion and ensuring sufficient resource utilization.

Under fixed buffer sizes = {_j}_j=1^n, new arrivals to queue j are blocked if x_j = _j. As a result, the state update is modified as follows,

    x_k+1 = min{x_k + De_k+1, }.


While a small  can greatly reduce congestion, it can impede the system throughput. To account for this, we introduce a cost for rejecting an arrival to the network. Let o_k∈{0,1}^n denote whether an arrival is overflowed, i.e., an arrival is blocked because the buffer is full,

    o_k+1 = De_k+1· 1{
    x_k + De_k+1 > }.

Given a fixed routing policy, the control task is to choose the buffer sizes  to minimize the holding and overflow costs:

    J_N(;ξ_1:N)    =  ∑_k=0^N-1  (h^⊤x_k)τ_k+1^* + b^⊤o_k.

Similar to the routing control problem, despite the fact that overflow is discrete, our gradient estimation framework is capable of computing a  gradient of the cost with respect to the buffer sizes, which we denote as ∇_ J_N(;ξ_1:N), i.e., we can evaluate gradients at integral values of the buffer size and use this to perform updates. Since the buffer sizes must be integral, 

we update the buffer sizes via sign gradient descent to preserve integrality:

    _t+1 = _t - 𝗌𝗂𝗀𝗇( ∇_ J_N(;ξ_1:N) ).






Learning for admission control has been studied in the queuing and simulation literature <cit.>. While exact gradient methods are possible in fluid models <cit.>, the standard approach for discrete queuing models is finite perturbation analysis <cit.>, given the discrete nature of the buffer sizes.
Randomized finite-differences, which is also known as Simultaneous Perturbation Stochastic Approximation (SPSA) <cit.>, is a popular optimization method for discrete search problems. This method forms a finite-differences gradient through a random perturbation. Let η∼𝖱𝖺𝖽𝖾𝗆𝖺𝖼𝗁𝖾𝗋(n, 1/2) ∈{-1,1 }^n be a random n-dimensional vector where each component is an independent 𝖱𝖺𝖽𝖾𝗆𝖺𝖼𝗁𝖾𝗋 random variable, taking values in {-1,1 } with equal probability. For each perturbation η, we evaluate the objective at ±η, i.e., J_N( + η; ξ_1:N) and J_N( - η; ξ_1:N), using the same sample path for both evaluations to reduce variance. For improved performance, we average the gradient across a batch of B perturbations, i.e., η^(b) for b=1,...,B, drawing a new sample path ξ^(b)_1:N for each perturbation. The batch SPSA gradient is 

    ∇_^𝖲𝖯𝖲𝖠,B J_N() = 1/B∑_b=1^B1/2( J_N( + η^(b); ξ^(b)_1:N) - J_N( - η^(b); ξ^(b)_1:N) ) η^(b).

We update the buffer sizes  according to the same sign gradient descent algorithm as in (<ref>). 


In comparison with existing works in the queuing literature (e.g. <cit.>), which derive analytical results for simple single-class or multi-class queues, we consider admission control tasks for large, re-entrant networks with multiple job classes. Each job class has its own buffer, resulting in a high-dimensional optimization problem in large networks. Moreover, the buffer size for one job class affects downstream congestion due to the re-entrant nature of the networks. For our experiments, we fix the scheduling policy to be the soft priority policy π^𝗌𝖯𝖱_θ(x) in (<ref>) due to its simplicity and strong performance in our environments. We emphasize however that our framework can be applied to buffer control tasks under any differentiable routing policy, including neural network policies. For each gradient estimator, we perform T=100 iterations of sign gradient descent, and each gradient estimator is computed from trajectories of length N = 1000. For SPSA, we consider batch sizes of B = {10 , 100, 1000} whereas we compute  with only B=1 trajectory. When evaluating the performance, we calculate the long-run average cost with the buffer size determined by the last iterate with a longer horizon N = 10^4 and over 100 trajectories. We also average the results across 50 runs of sign gradient descent.



The left panel of Figure <ref> displays iterates of the sign gradient descent algorithm with   for the M/M/1 queue with holding cost h = 1 and overflow cost b = 100. We observe that sign gradient descent with  (computed over a horizon of N=1000 steps) quickly reaches the optimal buffer size of L^* = 14 and remains there, oscillating between L = 14 and 15. The right panel shows the iterates for a simple 2-class queue with 1 server, h = 1, and b = 20 under a soft priority policy. We again observe that sign gradient descent with   quickly converges to a near-optimal set of buffer sizes.



To see how the estimator performs in larger-scale problems, we consider the Re-entrant 1 and Re-entrant 2 networks introduced in Section <ref> with varying number of job classes (i.e., varying number of layers). 

Figure <ref> compares the last iterate performance of SPSA and  for these two families of queuing networks with instances ranging from 6-classes to 21-classes. Holding costs are h = 1 and overflow costs are b = 1000 for all queues. We observe that  with only a single trajectory is able to outperform SPSA with B = 1000 trajectories for larger networks. Sign gradient descent using SPSA with only B = 10 trajectories is much less stable, with several of the iterations reaching a sub-optimal set of buffer sizes that assign L_j = 0 to several queues. This illustrates the well-known fact that for high-dimensional control problems, zeroth-order methods like SPSA must sample many more trajectories to cover the policy space and their performance can scale sub-optimally in the dimension. Yet , which is an approximate first-order gradient estimator, exhibits much better scalability with dimension and is able to optimize the buffer sizes with much less data.







§ POLICY PARAMETERIZATION




While our gradient estimation framework 
offers a sample-efficient alternative for 
learning from the environment, there is another practical 
issue that degrades the performance of learning algorithms for 
queuing network control: instability. Standard model-free RL algorithms are based on the `tabula rasa' principle, which aims to search over a general and unstructured policy class in order to find an optimal policy. However, it has been observed that this approach may be unsuitable for queuing network control. Due to the lack of structure, the policies visited by the algorithm often fail to stabilize the network,
which prevents the algorithm from learning and improving.
As a result, researchers have proposed structural modifications to ensure stability, 
including behavior cloning of a stabilizing policy to find a good initialization <cit.>, 
switching to a stabilizing policy 
if the queue lengths exceed some finite thresholds <cit.>, 
or modifying the costs to be stability-aware <cit.>. 


We investigate the source of instability in various queuing scheduling problems and find a possible explanation. 
We note that many policies obtained by model-free RL algorithms are not work-conserving and often allocate servers to empty queues. A scheduling policy is work-conserving if it always keeps the server(s) busy when there are compatible jobs waiting to be served.

Standard policies such as the cμ-rule, MaxPressure, and MaxWeight are all work-conserving, which partly explains their success in stabilizing complex networks. We treat work conservation as an `inductive bias'

and consider a simple modification to the policy architecture that guarantees this property without sacrificing the flexibility of the policy class. 

The de-facto approach for parameterizing policies in deep reinforcement learning is to consider a function ν_θ(x),
which belongs to a general function family, such as neural networks, and outputs real-valued scores.
These scores are then fed into a  layer, which converts the scores to probabilities over actions.
Naively, the number of possible routing actions can grow exponentially in the number of queues and servers.
Nonetheless, one can efficiently sample from the action space by having the output of
ν_θ(x) ∈^m× n be a matrix where row i, denote as ν_θ(x)_i, contains the scores for matching server i to different queues. Then by applying the  for row i, i.e., (ν_θ(x)_i), we obtain the probability that server i is assigned to each queue. We then sample the assignment independently for each server to obtain an action in 𝒰.
For the purpose of computing the  estimator, (ν_θ(x)_i) also gives a valid fractional routing in 𝒰. We let (ν_θ(x))∈𝒰 denote the matrix formed by applying the softmax to each row in ν_θ(x).



Under this `vanilla' softmax policy, the probability π_θ(x)_ij that server i is routed to queue j 
(or alternatively, the fractional capacity server i allocated to j) is given by

    π_θ(x)_ij = (ν_θ(x))_ij = e^ν_θ(x)_ij/∑_j=1^n e^ν_θ(x)_ij. Vanilla Softmax

Many of the policies mentioned earlier can be defined in this way, such as the soft MaxWeight policy, ν_θ(x)_i = {θ_jx_jμ_ij}_j=1^n. This parameterization is highly flexible
and ν_θ(x)
can be the output of a neural network. However, for a general ν_θ(x),
there is no guarantee that π(θ)(x)_i,j = 0 if x_j = 0. This means that such policies
may waste service capacity by allocating capacity to empty queues even when there are non-empty queues that 
server i could work on.

We propose a simple fix, which reshapes the actions produced by the policy.
We refer to this as the work-conserving ,




    π_θ^𝖶𝖢(x)_ij = 𝗌𝗈𝖿𝗍𝗆𝖺𝗑𝖶𝖢(ν_θ(x))_ij≡e^ν_θ(x)_ij1{x_j > 0}∧ϵ/∑_l=1^n e^ν_θ(x)_il1{x_l > 0}∧ϵ,
    WC-Softmax

where ∧ is the minimum and ϵ is a small number to prevent
division by zero when the queue lengths are all zero.

This parameterization is fully compatible with deep reinforcement learning approaches.
ν_θ(x) can be a neural network and critically, the work-conserving 
preserves the differentiability of π^𝖶𝖢_θ(x) with respect to θ.
As a result,  and  estimators can both be computed under this parameterization, since ∇_θlogπ_θ^𝖶𝖢(x) and ∇_θπ_θ^𝖶𝖢(x)
both exist.

This simple modification delivers substantial improvements in performance. Figure <ref> compares the average holding cost across policy iterations for PPO without any modifications, PPO initialized with a policy trained to imitate MaxWeight, and PPO with the work-conserving . Despite its empirical success in many other reinforcement learning problems, PPO without any modifications fails to stabilize the network and incurs an exceedingly high cost. It performs much better under an initial behavioral cloning step, which achieves stability but still underperforms the cμ-rule. On the other hand, with the work-conserving , even the randomly initialized policy stabilizes the network and outperforms the cμ-rule over the course of training. This illustrates that an appropriate choice of policy architecture, motivated by queuing theory, is decisive in enabling learning-based approaches to succeed. As a result, for all of the policy optimization experiments in sections <ref>, <ref>, and <ref>, we equip the policy parameterization with the work-conserving .





















































§ SCHEDULING FOR MULTI-CLASS QUEUING NETWORKS: BENCHMARKS




We now benchmark the performance of the policies obtained by  policy gradient (Algorithm <ref>) with standard queuing policies and policies obtained using state-of-the-art model-free reinforcement learning algorithms.
We consider networks displayed in Figure <ref>, which were briefly described in section <ref> and appeared in previous works <cit.>. <cit.> used Criss-cross and Re-entrant-1 networks to show that  can outperform standard queuing policies. <cit.> consider the Re-entrant-2 network, but did not include any RL baselines. We consider networks with exponential inter-arrival times and workloads in order to compare with previous results. We also consider hyper-exponential distributions to model settings with higher coefficients of variation, as has been observed in real applications <cit.>. The hyper-exponential distribution X∼𝖧𝗒𝗉𝖾𝗋𝖤𝗑𝗉(λ_1,λ_2,p) is a mixture of exponential distributions:

    X d= Y· E_1 + (1 - Y) · E_2,

for Y∼𝖡𝖾𝗋𝗇𝗈𝗎𝗅𝗅𝗂(p), E_1∼𝖤𝗑𝗉(λ_1), E_2∼𝖤𝗑𝗉(λ_2), and all are drawn independently of each other. We calibrate the parameters of the hyper-exponential distribution to have the same mean as the corresponding exponential distribution, but with a 1.5x higher variance.

Our empirical validation goes beyond the typical settings studied in the reinforcement learning for queuing literature, and is enabled by our discrete-event simulation framework.



We now describe the standard queuing policies considered in this section, which can all be expressed in the form

    π(x)=_u∈𝒰∑_i∈[m],j∈[n]ρ_ij(x)u_ij

for some index ρ∈^m × n that differs per method.

    
  * cμ-rule <cit.>: ρ_ij = h_jμ_ij1{x_j > 0 }. Servers prioritize queues with a higher holding cost and a larger service rate.
    
    
  * MaxWeight <cit.>: ρ_ij(x) = h_jμ_ijx_j. Servers prioritize queues that are longer, with a higher holding cost, and a larger service rate.
    
    
  * MaxPressure <cit.>: ρ_ij= ∑_ℓ=1^nμ_ijR_jℓh_ℓx_ℓ. MaxPressure is a modification of MaxWeight, in the sense that it takes workload externality within the network into account through the R_jl terms, e.g., processing a class j job may generate a new class j' job. 
    
    
  * Fluid <cit.>: The scheduling policy is based on the optimal sequence of actions in the fluid relaxation, which approximates the average evolution of stochastic queue lengths by deterministic ordinary differential equations. We aim to solve the continuous-time problem
    
    
    
    
    
    
    
    min_u̅   ∫_0^T h^⊤ x(t)dt 
    s.t.    ẋ(t) = λ - R(μ⊙u̅(t)),   ∀ t∈ [0,T]
        x(t) ≥ 0,   ∀ t∈ [0,T] 
       u̅(t) ∈𝒰,   ∀ t∈ [0,T].

    For tractability, we discretize the problem with time increment Δ t > 0 and horizon H=T/Δ t, and solve as a linear program.
    We then set u_k = u(t_k). The linear program is re-solved periodically to improve fidelity with the original stochastic dynamics.

We next describe the deep reinforcement learning methods considered in this section:

    
  * -𝖣𝖦 <cit.>:  is a standard model-free policy gradient method <cit.>.  <cit.> implement  for multi-class queuing networks and show that the policies obtained outperform several standard queuing policies. 
    Their implementation includes an initial behavioral cloning for stability and a carefully designed variance-reducing policy gradient estimation. We report the results from their paper, although in our experiments we include several problem instances not evaluated in their work.
    
  * -𝖶𝖢 (ours): Given the empirical success of  with the work-conserving , we use this algorithm as the main RL benchmark. For this policy, we use the same neural network architecture, hyper-parameters, and variance reduction methods as -𝖣𝖦  <cit.>.
    
  *  policy gradient (Algorithm <ref>): Trains a neural network policy with work-conserving  using the  policy  gradient estimator. We use an inverse temperature of β = 10 for all experiments in this section.

While value-based methods have also been considered for queuing network control <cit.>, our focus in this work is on benchmarking policy gradient algorithms. 

For the reinforcement learning policies, we train each method over 100 episodes, each consisting of N = 50,000 environment steps. Following <cit.>'s implementation, -𝖶𝖢 was trained with B = 50 actors, while  was trained only with B = 1 actor, which means that -based methods used 50x more trajectories than . See Appendix <ref> for more details on the training process.






To evaluate each scheduling policy, we run 100 parallel episodes starting from empty queues x_0 = 0_n with a long horizon N to estimate the long-run average holding cost (typically N =200,000 steps). As in previous works (e.g. <cit.>), we consider holding costs h = _n in which case the holding cost is equivalent to the total queue length. To reiterate, for each policy π we estimate the following quantity 

    J_N(π) = 
      [1/t_N∑_k=0^N-1 (_n^⊤x_k)τ^*_k+1]
    =  [1/t_N∫^t_N_0∑_i=1^n x_i(t)dt ]

where t_N is the time of the Nth event. We measure the standard deviation across the 100 episodes to form 95% confidence intervals. For the reinforcement learning policies, we report the average holding cost for the best policy encountered during training.






Tables 1-5 display the results of our benchmarking across the problem instances discussed before. The column `% improve' records the relative reduction in holding cost achieved by  over the best standard queuing policy (either cμ, MaxWeight, MaxPressure, or Fluid). Our main observations on the relative performance of the standard policies and policies obtained from reinforcement learning methods are summarized as follows.

    
  *  is a strong reinforcement-learning benchmark.  with our proposed work-conserving  is able to efficiently find policies that outperform standard policies across all problem instances, as well as the  featured in <cit.> (under the same policy network and hyper-parameters). This illustrates that simply ensuring work-conservation is a powerful inductive bias that delivers stability.
    
  *  policy gradient outperforms  in larger networks, using 50x less data. We observe that  and  achieve similar performances when the number of classes is small. However, when the number of job classes gets larger,  consistently outperforms , as seen in Figure <ref>. This is likely due to the 
    
    sample efficiency gained by the  gradient, enabling the algorithm to find better policies with less data.
    
  *  achieves large performance gains over   for higher-variance problem instances, using 50x less data. We observe that for the Re-entrant-1 networks with hyper-exponential noise, the reduction of holding cost of  relative to  is often equivalent or even larger than the cost reduction of  relative to the cμ-rule. This illustrates that even among optimized RL policies, there can be significant performance differences for difficult problem instances, and the sample efficiency of  is particularly useful in noisier environments.
    
  * While standard queuing methods work well, RL methods meaningfully improve performance in hard instances. We observe in the `% improve' column that  achieves a 3-20% improvement over the best standard queuing policy in each setting.




Altogether, these results illustrate that policy gradient with  gradient estimator and work-conserving  policy architecture can learn effective queuing network control policies with substantially less data than model-free policy gradient algorithms for large networks with high-variance event times, mirroring real-world systems. 
In particular, the improved sample efficiency of the  gradient estimator and the stability brought by work-conserving  policy architecture are the keys to enabling learning in large-scale systems with realistic data requirements.









§ WHY IS  SAMPLE-INEFFICIENT? A THEORETICAL CASE STUDY FOR THE M/M/1 QUEUE



In this section, we provide a theoretical case study explaining how 
gradients are able to learn more from a single observed trajectory compared to
. We focus on the special case of the M/M/1 queue to explain how
 and its actor-critic variants utilizing baselines/advantages suffer from
sample inefficiency.
Although we are only able to analyze a substantially simpler setting than the control problems in general multi-class queueing networks we are interested in, 
our theoretical results illustrate the essential statistical benefits of pathwise gradient estimators.
We highlight that while  applies to virtually any setting by relying on random exploration, 
it fundamentally struggles to assign credit to actions, especially in noisy environments. , on the other hand, is much better at assigning credit to actions. This allows us to crystallize why we see such a large improvement in sample efficiency in sections <ref> and <ref>. We also believe our result may be of broader interest in reinforcement learning, because it illustrates a practically relevant instance where , even with an optimal baseline, is provably sub-optimal. 

We consider the M/M/1 queue with a fixed arrival rate λ under service rate control u=μ > λ. This setting permits an analytic approach to showing that the  estimator has a sub-optimally large variance, particularly for congested systems when ρ = λ/μ→ 1. On the other hand, the  estimator achieves an order of magnitude improvement in estimation efficiency. 
In the M/M/1 queue setting, the  approach
for general queuing networks 
reduces to IPA based on Lindley recursion. 




We consider a simple service-rate control problem where the cost is the steady-state average queue length in the M/M/1 queue

    Q(μ) := _∞[x(t)] = λ/μ - λ = ρ/1-ρ.

Given that the service rate is continuous, it is natural to consider a policy that randomizes over [μ - h, μ]. One such option is a Beta-distributed policy π_θ: A = μ - hY, where Y ∼Beta(θ,1) and h>0. As θ→∞, the policy frequently sets service rates close to μ and as θ→ 0, it concentrates more probability mass on service rates close to μ - h. The task is to estimate the following policy gradient

    ∇_θ J(θ) = ∇_A∼π_θ[Q(A)].

The 𝖱𝖤𝖨𝖭𝖥𝖮𝖱𝖢𝖤 estimator of ∇_θJ(θ) involves sampling a random service rate from the policy π_θ, and then estimating the steady-state queue length from a trajectory under that service rate. For a trajectory with N steps, we denote the steady-state queue length estimator as Q_N(A). Then, the  gradient estimator takes the form

    ∇^𝖱 J_N(θ; ξ_1:N)
    = Q_N(μ - hY) ∇_θlogπ_θ(Y) 
    = Q_N(μ - hY) ( log Y + 1/θ).

The standard estimator for the steady-state queue length is simply the queue length averaged over a sample path: 
Q_N(a) = 1/N∑_k=1^Nx_kτ_k+1^* when A=a. As long as μ - h > λ, it is known that as N→∞, Q_N(a) → Q(a), which implies that ∇^𝖱 J_N(θ; ξ_1:N) →∇ J(θ).


On the other hand, the  estimator utilizes the structure of the single server queue. First, by inverse transform sampling, Y d= F_θ^-1(ω), where ω∼Uniform(0,1) and F_θ^-1(ω) = ω^1/θ. Then, we can substitute A = μ - hω^1/θ. Since Q(μ) is differentiable and the derivative is integrable, we can change the order of differentiation and integration

    ∇_θ J(θ) 
    = ∇_θ_ω[Q(μ - hω^1/θ)]
    = -_ω[∇ Q(μ - hω^1/θ) · h ∇_θω^1/θ].

The preceding display involves the gradient of the steady-state queue-length Q(μ) with respect to the service rate μ, i.e., ∇ Q(μ). 

For the M/M/1 queue, there are consistent sample-path estimators of ∇ Q(μ). One such estimator uses the fact that by Little's law Q(μ) = _∞[x(t)] = λ_∞[w(t)] =: λ W(μ) where W(μ) is the steady-state waiting time. The waiting time process W_i, which denotes the waiting time of the ith job arriving to the system, has the following dynamics, known as the Lindley recursion:

    W_i+1 = ( W_i - T_i+1 + S_i/μ)^+,

where T_i+1iid∼Exp(1) is the inter-arrival time of between the ith job and (i+1)th job, and S_iiid∼Exp(1) is the workload of the ith job. Crucially, this stochastic recursion specifies how the service rate affects the waiting time along the sample path, which enables one to derive a pathwise derivative via the recursion:

    ∇W_i+1 = 
     ( -S_k/μ^2 + ∇ W_i)
     1{ W_i+1 > 0 },

where ∇ W_i together with W_i form a Markov chain following the above recursion. By averaging this gradient across jobs and using Little's law, we have the following gradient estimator for ∇ Q, which we denote as ∇ Q_N(μ): 


    ∇ Q_N(μ) = λ1/L_N∑_i=1^L_N∇ W_i,

where L_N is the number of arrivals that occur during a sample path with N events.
Using this, the  policy gradient estimator (a.k.a. IPA estimator) is

    ∇_θ J_N(θ; ξ_1:N) 
    = h ·∇ Q_N(μ - h Y) ( 1/θ Y log Y ).

As long as μ - h > λ, it has been established that ∇ Q_N(μ) is asymptotically unbiased, i.e., as N→∞, [∇ Q_N(μ)] →∇ Q(μ), which implies 
[∇_θ J_N(θ; ξ_1:N)] →∇ J(θ). 

Since both the  and  gradient estimators give an asymptotically unbiased estimation of ∇ J(θ), we compare them based on their variances, which determines how many samples are needed to reliably estimate the gradient. Although the variances of the estimators are not precisely known for a finite N, Q_N and ∇Q_N both satisfy the central limit theorem (CLT) with explicitly characterized asymptotic variances, which we denote as _∞(Q) and _∞(∇Q) 

respectively. This implies that the variance of Q_N is approximately _∞(Q)/N. We define _∞(∇ J_N(θ; ξ_1:N)) to be the variance of the gradient estimator when we approximate the variance of Q_N and ∇Q_N using _∞(Q)/N and _∞(∇Q)/N respectively.


It is worth reiterating that the  estimator only required an estimate of cost Q(μ), which does not require any domain knowledge, whereas the  gradient required an estimate of ∇ Q(μ) which requires a detailed understanding of how the service rate affected the sample path dynamics. Utilizing this structural information can greatly improve the efficiency of gradient estimation. Since ∇_θ J_N(θ; ξ_1:N) = O(h) almost surely,
we have (∇_θJ_N(θ; ξ_1:N)) = O(h^2). On the other hand, the variance of the  estimator can be very large even if h is small. To highlight this, consider the extreme case where h = 0 for which the policy gradient ∇ J(θ) is obviously zero since the policy deterministically sets the service rate to μ regardless of θ. Strikingly,  does not have zero variance in this case:




If h = 0, then the variance of the estimators are

    _∞(∇ J_N(θ; ξ_1:N) ) = 0,   _∞(∇^𝖱 J_N(θ; ξ_1:N) )
    = Θ( N^-1 (1 - ρ)^-4)


Note that even when h=0, the variance of the  estimator can be quite high if the queue is congested, i.e. ρ is close to 1, while the pathwise estimator gives the correct estimate of zero with zero variance.

For non-trivial values of h, we focus on the so-called `heavy-traffic' asymptotic regime with ρ=λ/μ→ 1, which is of major theoretical and practical interest in the study of queues. Estimating steady-state quantities becomes harder as the queue is more congested, so (1 - ρ)^-1 emerges as a key scaling term in the variance. 

We set h such that h/μ - λ<c and h/μ - λ→ c ∈ (0,1) as ρ→ 1. This resembles the square-root heavy-traffic regime for capacity planning, where the service rate is set to be λ + β√(λ) for some β > 0, and one considers the limit as λ→∞. In this case, if one were choosing a policy over the square-root capacity rules A∈[λ + a √(λ), λ + b √(λ)] for some b>a>0, this is equivalent to setting μ = λ + b√(λ) and h = (b - a) √(λ) = O(√(λ)). Note that if c = 0, the gradient is zero (identical to Observation <ref>), and if c ≥ 1, the queue with service rate μ-h is unstable.

Within this regime, we have the following comparison between the gradient estimators, which utilizes recent results concerning the asymptotic variance of ∇Q_N(μ) <cit.>.


Suppose h = c(μ - λ) for c∈(0,1) as ρ→ 1. Under this scaling
∇ J(θ) ∼ (1-ρ)^-1,
and

    _∞(∇ J_N(θ; ξ_1:N) ) 
       = O ( N^-1(1 - ρ)^-3_estimation noise + (1-ρ)^-2_policy randomization) 
    _∞(∇^𝖱 J_N(θ; ξ_1:N) )    = Θ( N^-1(1 -ρ)^-4
    + (1-ρ)^-2)


See Appendix <ref> for the proof.

Overall, the  estimator is much more sample efficient than the  estimator as ρ→ 1, with the variance scaling as (1-ρ)^-3 compared to (1-ρ)^-4.
The first terms in (<ref>) and (<ref>) represent the variance occurring from the Monte Carlo estimation and becomes smaller if one generates a longer sample path (larger N), and it scales as N^-1. The second terms are the variance resulting from randomness in the service rate induced by the policy. 







Theorem <ref> illustrates that large improvements in statistical efficiency can be achieved by leveraging the structure of the system dynamics. An existing strategy for incorporating domain knowledge in  is to subtract a baseline b from the cost, which preserves un-biasedness:

    ∇^𝖱𝖡 J_N
    (θ, ξ_1:N)
    = (Q_N(μ - hY) - b) ( log Y + 1/θ).


In this case, one can characterize the optimal variance-reducing baseline in closed form if one has knowledge of the true cost Q(μ). Under h = c(μ - λ),

    b^*   = [Q(h-hY) ∇_θlogπ_θ(Y)^2]/[∇_θlogπ_θ(Y)^2]
       = λ/μ - λ[F^2_1(1,θ,1+θ,c) -2θ^2Φ(c,2,θ) + 2cθ^3Φ(c,3,1+θ)] 
        = O((1-ρ)^-1)

where F^2_1 is the hypergeometric 2F1 function and Φ is the Lerch Φ transcendental. The optimal baseline b^* is of the same order as Q(μ) as ρ→ 1.



Consider the  estimator with the optimal baseline b^*. As ρ→ 1, the variance of the estimator scales as

    _∞(∇^𝖱𝖡 J_N(θ; ξ_1:N) ) = Θ( N^-1(1 -ρ)^-4 
    +(1-ρ)^-2)



The proof of Corollary <ref> is provided in Appendix <ref>.
Simply, since the optimal baseline is a deterministic input, it is unable to improve upon the (1-ρ)^-4 dependence on ρ, which is driven by the statistical properties of Q_N. 

This illustrates that the pathwise gradient estimator can offer an order of magnitude improvement in sample efficiency than the  estimator even with an optimized baseline that requires knowledge of the true cost function (and thus precludes the need to estimate the cost in the first place). 

Intuitively, the  estimator is inefficient because it is unable to leverage the fact that 
Q(μ) ≈ Q(μ + ϵ) when ϵ is small. After all, generic MDPs do not have such a structure; a slight change in the action can result in vastly different outcomes. The  estimator cannot use the estimate of Q_N(μ) to say anything about Q(μ + ϵ), and must draw a new sample path to estimate Q(μ + ϵ). Meanwhile, using a single sample path, the pathwise estimator can obtain an estimate for Q(μ + ϵ) when ϵ is small enough via Q_N(μ + ϵ)≈Q_N(μ) + ϵ∇Q_N(μ). In this sense, the pathwise estimator can be seen as a infinitesimal counterfactual of the outcome under alternative—but similar—actions.

Even though we only study a single server queue here, 
we believe the key observations may apply more broadly.


  * Higher congestion (ρ→ 1) makes it more challenging to estimate the performance of queueing networks based on the sample path. This applies to both gradient estimators and baselines that could be used to reduce variance.

  * It is important to reliably estimate the effects of small changes in the policy, as large changes can potentially cause instability. Pathwise gradient estimators provide a promising way to achieve this. For general networks with known dynamics, their dynamics are often not differentiable, which requires the development of the  estimator.






§ CONCLUSION



In this work, we introduce a new framework for policy optimization in queuing network control. This framework uses a novel approach for gradient estimation in discrete-event dynamical systems. Our proposed  policy gradient estimator is observed to be orders of magnitude more efficient than model-free RL alternatives such as  across an array of carefully designed empirical experiments. In addition, we introduce a new policy architecture, which drastically improves stability while maintaining the flexibility of neural network policies. Altogether, these illustrate how structural knowledge of queuing networks can be leveraged to accelerate reinforcement learning for queuing control problems.

We next discuss some potential extensions of our approach:

    
    
  * We consider policies with preemption. Our proposed method can also handle non-preemptive policies by keeping track of the occupied servers as part of the state.
    
    
  * We focus on scheduling and admission control problems in queuing network satisfying Assumptions <ref>, but the algorithmic ideas can be extended to more general queuing networks by utilizing a larger state space that contains the residual workloads of all jobs in the network, rather than only the top-of-queue jobs as is done in this work. A higher dimensional state descriptor is required for more general networks as multiple jobs in the same queue can be served simultaneously.
    
  * Beyond queuing network control, our methodology can be extended to control problems in other discrete-event dynamical systems. More explicitly, our methodology can handle systems that involve a state update of the form x_k+1 = g(x_k, e_k+1) where g is a differentiable function and e_k+1 is the selected event. Recall that in this work, the state update is linear in x_k and e_k+1: x_k+1 = x_k + De_k+1. We also require that e_k+1 is differentiable almost surely in the action u_k. 









abbrvnat







§ TRAINING DETAILS



 WC was trained over 100 episodes, each consisting of 50,000 environment steps parallelized over 50 actors. We closely follow the hyper-parameters and training setup as in <cit.>. We used a discount factor of 0.998, a GAE <cit.> parameter of 0.99, and set the Kullback–Leibler divergence penalty as 0.03. For the value network, we used a batch size of 2,500, while for the policy network, we used the entire rollout buffer (batch size of 50,000) to take one gradient step. We performed 3 PPO gradient updates on the same rollout data. For all the experiments, we used the Adam optimizer with a cosine decaying warming-up learning rate scheduler. The learning rates were set to 3 × 10^-4 for the value network and 9 × 10^-4 for the policy network. We used 3% of the training horizon to warm up to the maximum learning rate and then cosine decayed to 1 × 10^-5 for both networks. We used the same neural network architecture as those in <cit.>, see Appendix E of  <cit.> for more details.

For the  policy gradient, we used the same hyperparameters across all experiments. We use an inverse temperature of β = 10 for the  relaxation of the event selection. We update the policy after every episode with the Adam optimizer, using constant step-size of 5× 10^-4, momentum parameters (0.8, 0.9), and gradient clipping of 1. For the policy neural network, we used a multilayer perceptron with 3 hidden layers, each hidden layer consisting of 128 hidden units. We use the work-conserving  for the final output.








    





§ PROOFS





 §.§ Proof of Theorem <ref>


We focus on the case where x_k≥1. Then,

    x_k+1    =x_k+1{τ_k^A<w_k/μ}-1{τ_k^A>w_k/μ}
        =x_k+e^-βτ_k^A/e^-βτ_k^A+e^-β w_k/μ-e^-β w_k/μ/e^-βτ_k^A+e^-β w_k/μ
        =x_k+e^-βτ_k^A-e^-β w_k/μ/e^-βτ_k^A+e^-β w_k/μ.

Since the inter-arrival times and workloads are exponentially distributed, by the memoryless property, we have τ_k^A∼𝖤𝗑𝗉(λ) and
w_k∼𝖤𝗑𝗉(1).


The true gradient is

    d/dμ𝔼[x_k+1-x_k]=d/dμλ-μ/λ+μ=-2λ/(λ+μ)^2.

Under our 𝗌𝗈𝖿𝗍𝗆𝗂𝗇_β approximation for the event-selection, 





we have

    𝔼[d/dμe^-βτ_k^A-e^-β w_k/μ/e^-βτ_k^A+e^-β w_k/μ]
        =𝔼[-2βe^-β(τ_k^A+w_k/μ)/(e^-βτ_k^A+e^-β w_k/μ)^2w_k/μ^2]
        =-2β/μ𝔼[τ_k^S(e^β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^21{τ_k^A<τ_k^S}+e^β(τ_k^S-τ_k^A)/(e^β(τ_k^S-τ_k^A)+1)^21{τ_k^A>τ_k^S})],

for τ_k^S=w_k/μ.






Next, note that

    𝔼[τ_k^Se^β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^21{τ_k^A<τ_k^S}]
        =𝔼[.𝔼[τ_k^Se^β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^21{τ_k^A<τ_k^S}] |τ_k^A=t]
        =𝔼[.𝔼[τ_k^Se^β(t-τ_k^S)/(e^β(t-τ_k^S)+1)^2|τ_k^A=t,τ_k^S>t]ℙ(τ_k^S>t|τ_k^A=t)]
        =𝔼[𝔼[(t+S')e^-β S'/(e^-β S'+1)^2]ℙ(τ_k^S>t|τ_k^A=t)] 
        =𝔼[(tA(β,μ)+B(β,μ))ℙ(τ_k^S>t|τ_k^A=t)]
        =𝔼[(τ_k^A A(β,μ)+B(β,μ))e^-μτ_k^A]
        =λ/(λ+μ)^2A(β,μ)+λ/(λ+μ)B(β,μ)

where

    A(β,μ)     =μ(β-μ H(μ/2β)+μ H(μ/2β-1/2))/2β^2
        =μ/2β+μ^2/2β^2(H(μ/2β-1/2)-H(μ/2β)_H̃(β,μ))
    
    B(β,μ)     =μ/4β^3(2β H(μ/2β)-2β H(μ/2β-1/2)-μψ^(1)(β+μ/2β)+μψ^(1)(2β+μ/2β))
        =-μ/2β^2H̃(β,μ)+μ^2/4β^3(ψ^(1)(2β+μ/2β)-ψ^(1)(β+μ/2β)_ψ̃^(1)(β,μ))

Moreover, note that

    H(μ/2β)     =ψ^(0)(μ/2β+1)+γ
    
    H(μ/2β-1/2)     =ψ^(0)(μ/2β+1/2)+γ

Note that

    H(μ/2β)-H(μ/2β-1/2)     =log(μ/2β+1)+1/μ/2β+1
        -log(μ/2β+1/2)+1/μ/2β+1.


Similarly,

    𝔼[τ_k^Se^β(τ_k^S-τ_k^A)/(e^β(τ_k^S-τ_k^A)+1)^21{τ_k^A>τ_k^S}]
        =𝔼[𝔼[.se^-β(τ_k^A-s)/(e^-β(τ_k^A-s)+1)^21{τ_k^A>s}|τ_k^S=s]]
        =𝔼[𝔼[.se^-β(τ_k^A-s)/(e^-β(τ_k^A-s)+1)^2|τ_k^A>s,τ_k^S=s]ℙ(τ_k^A>s|τ_k^S=s)]
        =𝔼[𝔼[se^-β T'/(e^-β T'+1)^2]ℙ(τ_k^A>s|τ_k^S=s)] 
        =𝔼[τ_k^S A(β,λ)e^-λτ_k^S]
        =μ/(λ+μ)^2 A(β,λ).

Then,

    -2β/μ𝔼[τ_k^S(e^β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^21{τ_k^A<τ_k^S}+e^β(τ_k^S-τ_k^A)/(e^β(τ_k^S-τ_k^A)+1)^21{τ_k^A>τ_k^S})]
        =-2β/μ(λ/(λ+μ)^2A(β,μ)+λ/(λ+μ)B(β,μ)+μ/(λ+μ)^2A(β,λ))
        =-2β/μ(λ A(β,μ)+μ A(β,λ)/(λ+μ)^2+λ/(λ+μ)B(β,μ))
        =-2λ/(λ+μ)^2-2β/μ(λμ^2H̃(β,μ)+μλ^2H̃(β,λ)/2β^2(λ+μ)^2-λμ/2β^2(λ+μ)H̃(β,μ)+λ/(λ+μ)μ^2/4β^3ψ̃^(1)(β,μ)).

Note that as β→∞,

    lim_β→∞H̃(β,μ)     =γ+ψ^(0)(1/2),
    lim_β→∞ψ̃^(1)(β,μ)     =-π^2/3.

This means that the leading order term is O(1/β^2). In particular,

    -2β/μ(λμ^2H̃(β,μ)+μλ^2H̃(β,λ)/2β^2(λ+μ)^2-λμH̃(β,μ)/2β^2(λ+μ))    ∼π^2λ^2(μ-λ)/6β^2(λ+μ)^2

Finally, we have the second-order term

    -2β/μλ/(λ+μ)μ^2/4β^3ψ̃^(1)(β,μ)∼π^2λμ/6β^2(λ+μ).

Thus, we have the following characterization of the bias:

    𝔼[d/dμe^-βτ_k^A-e^-β w_k/μ/e^-βτ_k^A+e^-β w_k/μ]-(-2λ/(λ+μ)^2)∼1/β^2π^2λ(μ^2-λ^2+2μλ)/6(λ+μ)^2+o(1/β^2)


For variance, we have

    𝔼[(-2βe^β(τ_k^A+w_k/μ)/(e^βτ_k^A+e^β w_k/μ)^2w_k/μ^2)^2]
        =𝔼[4β^2/μ^2e^2β(τ_k^A+τ_k^S)/(e^βτ_k^A+e^βτ_k^S)^4τ_k^S,2]
        =4β^2/μ^2𝔼[τ_k^S,2(e^2β(τ_k^A+τ_k^S)/(e^βτ_k^A+e^βτ_k^S)^41{τ_k^A<τ_k^S}+e^2β(τ_k^A+τ_k^S)/(e^βτ_k^A+e^βτ_k^S)^41{τ_k^A>τ_k^S})]
        =4β^2/μ^2𝔼[τ_k^S,2(e^2β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^41{τ_k^A<τ_k^S}+e^2β(τ_k^S-τ_k^A)/(e^β(τ_k^S-τ_k^A)+1)^41{τ_k^A>τ_k^S})].

Note that

    𝔼[τ_k^S,2e^2β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^41{τ_k^A<τ_k^S}]
        =𝔼[𝔼[.τ_k^S,2e^2β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^41{t<τ_k^S}]|τ_k^A=t]
        =𝔼[𝔼[(t+S')^2e^-2β S'/(e^-β S'+1)^4]ℙ(τ_k^S>t|τ_k^A=t)] 
        =𝔼[𝔼[(t^2+2S'+S'^2)e^-2β S'/(e^-β S'+1)^4]ℙ(S_i>t|T_i=t)]
        =𝔼[(t^2Ã(β,μ)+tB̃(β,μ)+C̃(β,μ))ℙ(S_i>t|T_i=t)]
        =𝔼[(T_i^2Ã(β,μ)+T_iB̃(β,μ)+C̃(β,μ))e^-μ T_i]
        =2λ/(λ+μ)^3Ã(β,μ)+λ/(λ+μ)^2B̃(β,μ)+λ/(λ+μ)C̃(β,μ).

Similarly,

    𝔼[τ_k^S,2e^2β(τ_k^S-τ_k^A)/(e^β(τ_k^S-τ_k^A)+1)^41{τ_k^A>τ_k^S}]
        =𝔼[𝔼[.s^2e^2β(s-τ_k^A)/(e^β(s-τ_k^A)+1)^41{τ_k^A>s}|τ_k^S=s]]
        =𝔼[𝔼[s^2e^-2β(τ_k^A-s)/(e^-β(τ_k^A-s)+1)^4|τ_k^A>s,τ_k^S=s]ℙ(τ_k^A>s|τ_k^S=s)]
        =𝔼[𝔼[s^2e^-2β T'/(e^-β T'+1)^4]ℙ(τ_k^A>s|τ_k^S=s)] 
        =𝔼[Ã(β,λ)τ_k^S,2e^-λτ_k^S]
        =2μ/(λ+μ)^2Ã(β,λ).

Putting the above two parts together, we have

    4β^2/μ^2𝔼[τ_k^S,2(e^2β(τ_k^A-τ_k^S)/(e^β(τ_k^A-τ_k^S)+1)^41{τ_k^A<τ_k^S}+e^2β(τ_k^S-τ_k^A)/(e^β(τ_k^S-τ_k^A)+1)^41{τ_k^A>τ_k^S})]
        =4β^2/μ^2(2λ/(λ+μ)^3Ã(β,μ)+λ/(λ+μ)^2B̃(β,μ)+λ/(λ+μ)C̃(β,μ)+2μ/(λ+μ)^2Ã(β,λ))
       ∼4βλ/3μ(λ+μ)^2





 §.§ Proof of Theorem <ref>



By assumption, h ≤ c(μ - λ) for some c < 1. First, we develop bound for _∞ (∇̂^𝖱 J_N(θ; ξ_1:N)).

We can compute the variance by conditioning on the value of Y:

    _∞(Q̂_N(μ-hY)·(log Y-1/θ))   =𝔼[_∞(Q̂_N(μ-hy)(log y-1/θ)|Y=y)]
       +(𝔼[Q̂_N(μ-hy)(log y-1/θ)|Y=y)]).

For the first term, note that the asymptotic variance in the CLT for the ergodic estimator Q̂_N(μ) is _∞(Q̂(μ))=2ρ(1+ρ)/(1-ρ)^4. Then, we have _∞(Q̂_N(μ)) = 2ρ(1+ρ)/N(1-ρ)^4. Since h ≤μ, _∞(Q̂_N(μ))= 2ρ(1+ρ)/N(1-ρ)^4.
Then,

    [ _∞(Q̂_N(μ-hy)·(log y+1/θ)|Y = y) ] 
        =[  (log y+1/θ)^2_∞(Q̂_N(μ-hy)|Y = y) ]
       ≥[  (log Y+1/θ)^22ρ(1+ρ)/N(1-ρ)^4]
        = 1/θ^22ρ(1+ρ)/N(1-ρ)^4=Θ( (1-ρ)^-4),

where the last equality uses the fact that for Y∼Beta(θ,1),

    [log Y]=ψ(θ) - ψ(θ + 1)=-1/θ

and

    [  (log Y+1/θ)^2]=(log Y) = ψ_1(θ) - ψ_1(θ + 1) = 1/θ^2.

For the second term, we plug in the true estimand as the expectation of Q_N(μ), i.e. Q(μ) = λ/μ - λ. Then,

    (𝔼[Q̂_N(μ-hy)(log y+1/θ)|Y=y)]) 
        = (λ/μ - hY - λ( log Y + 1/θ)) 
        = [ (λ/μ - hY - λ)^2( log Y + 1/θ)^2]
    - [ (λ/μ - hY - λ)
    ( log Y + 1/θ) ]^2.

We proceed to evaluate these expectations analytically.

    [ (λ/μ - hY - λ)
    ( log Y + 1/θ) ] 
       = ρ/1-ρ(
    Γ(θ)/Γ(1+θ)
    F^2_1(1,θ,1+θ,h/μ-λ)
    -
    θΦ(h/μ - λ,2,θ)
    ) 
       =O((1-ρ)^-1),

where F^2_1 is the Hypergoemetric 2F1  function and Φ is the Lerch transcendental function.
We also have

    [ (λ/μ - hY - λ)^2( log Y + 1/θ)^2] 
        = λ^2/θ^2(μ - λ)^3(
    2(μ - λ)
    + hθ^3Φ(h/μ - λ,2,θ+1 )
    + hθ^3(θ - 1)Φ(h/μ - λ,3,θ+1 ) . 
        
    + (μ - λ) 
    ( θμ - λ/μ - h - λ - (θ - 1) F_1^2(1,θ,1+θ,h/μ-λ)
    ) 
       . +2(μ -λ) F_2^3((2,θ,θ),(1+θ,1+θ),h/μ-λ)
    ) 
        = O((1-ρ)^-2).

Taking the difference between the above two parts under the limit as (1-ρ)→ 0, we have

    (λ/μ - hY - λ( log Y + 1/θ)) =O((1-ρ)^-2).


Next, we develop a bound for _∞ (∇̂ J_N(θ; ξ_1:N)).

    _∞(
        h ·∇̂ Q_N(μ - h Y) ( 1/θ Y log Y )
        )   =
        [_∞(h ·∇̂ Q_N(μ - h Y) ( 1/θ Y log Y )|Y=y)]
       +(𝔼[h ·∇̂ Q_N(μ - h Y) ( 1/θ Y log Y )|Y = y])

For the first term, we can use the fact that |Ylog Y| ≤ 1/e almost surely since Y∈[0,1]. We also use recent results in <cit.>, which compute the asymptotic variance of the IPA estimator:

    _∞(∇̂ Q_N(μ)) = 1 + 16ρ + 27 ρ^2 + 2 ρ^3 + 6ρ^4/μ^2N(1 + ρ)(1-ρ)^5≤ 52μ^-2N^-1(1-ρ)^-5

Under μ - hy, the congestion factor 1 - λ/μ - hy = μ - hy - λ/μ - hy≥μ - cy(μ - λ) - λ/μ = (1-cy)(1-ρ) and μ-hy≥μ-h≥ (1-c)μ.
So we have the bound,

    _∞(h∇̂ Q_N(μ - h Y) ( 1/θ y log y)|Y=y) 
       ≤  h^2θ^-2(y log y)^2_∞(∇̂ Q_N(μ - h y)) 
       ≤ h^2μ^-2θ^-2e^-2 52N^-1(1-c)^-7(1-ρ)^-5
        = O(N^-1 h^2μ^-2 (1-ρ)^-5) 
       = O(N^-1  (1-ρ)^-3)

since h = O(1-ρ).

For the second term, we plug in the true estimand as the mean of ∇̂ Q_N(μ), i.e., ∇ Q(μ) = -ρ/μ(1-ρ)^2,

    [h ·∇̂ Q_N(μ - h Y) ( 1/θ Y log Y )|Y = y] 
    = - h1/θ (y log y) λ/(μ - hy - λ)^2.


We next evaluate the variance analytically,

    (
    - h1/θ (Y log Y) λ/(μ - hY - λ)^2) 
       = 
    [ 
    (- h1/θ (Y log Y) λ/(μ - hY - λ)^2)^2] 
    - [ 
    - h1/θ (Y log Y) λ/(μ - hY - λ)^2]^2

Since

    [ 
    - h1/θ (Y log Y) λ/(μ - hY - λ)^2] 
       = h/(1+θ)^2λ/(μ - λ)^2
    F^3_2((2,1+θ,1+θ), 
    (2+θ, 2+ θ), h/μ -λ) 
        = O((1-ρ)^-1)

and

    [ 
    (- h1/θ (Y log Y) λ/(μ - hY - λ)^2)^2] 
       =
    2hθ^2Γ(θ)^3Γ(2+θ)^-3λ^2/(μ - λ)^2×( F^4_3((3,1+θ,1+θ,1+θ), 
    (2+θ, 2+ θ,2+ θ), h/μ -λ) . 
       . -F^4_3((4,1+θ,1+θ,1+θ), 
    (2+θ, 2+ θ,2+ θ), h/μ -λ) ) 
        = O((1-ρ)^-2),

we have

    (
    - h1/θ (Y log Y) λ/(μ - hY - λ)^2) = O( (1-ρ)^-2).





 §.§ Proof of Corollary <ref>



First, we can explicitly characterize the optimal baseline:

    b^*   = [Q(μ -hY) ∇_θlogπ_θ(Y)^2]/[∇_θlogπ_θ(Y)^2]
       = [Q(μ -hY) ( log Y + 1/θ)^2]/[( log Y + 1/θ)^2]
       = λ/μ - λ[F^2_1(1,θ,1+θ,h/μ-λ) -2θ^2Φ(h/μ-λ,2,θ) + 2h/μ-λθ^3Φ(h/μ-λ,3,1+θ)]_b(θ)
        = O((1-ρ)^-1).

Next, we plug this into the  estimator. 

    _∞((Q̂_N(μ-hY)-b^*)·(log Y-1/θ))   =𝔼[_∞((Q̂_N(μ-hy)-b^*)(log y-1/θ)|Y=y)]
       +(𝔼[(Q̂_N(μ-hy) - b^*)(log y-1/θ)|Y=y]).

Note that since b^* is a constant, the first term, i.e., the mean of the conditional variance given Y, has the same value as in Theorem <ref>. For the second term, note that since h = c(μ - λ),


    [ (λ/μ - hY - λ  - b^*)
    ( log Y + 1/θ) ]  
       = (λ/μ - λ)^2[ (1/1- cY  - b(θ))
    ( log Y + 1/θ) ].

Since [ (1/1- cY  - b(θ))
( log Y + 1/θ) ] > 0 and doesn't depend on μ or λ, this confirms that the second term is Θ((1-ρ)^-2)
















