Dynamic Control of Queuing Networks
via Differentiable Discrete-Event Simulation
Ethan Che * 1 Jing Dong 1 Hongseok Namkoong 1

Abstract
Queuing network control is a problem that arises
in many applications such as manufacturing, communications networks, call centers, hospital systems, etc. Reinforcement Learning (RL) offers abroad set of tools for training controllers for general queuing networks, but standard model-free
approaches suffer from high variance of trajectories, large state and action spaces, and instability.
In this work, we develop a modeling frameworkfor queuing networks based on discrete-event simulation. This model allows us to leverage tools
from the gradient estimation literature to compute
approximate first-order gradients of sample-path
performance metrics through auto-differentiation,
despite discrete dynamics of the system. Using
this framework, we derive gradient-based RL algorithms for policy optimization and planning.
We observe that these methods improve sample
efficiency, stabilize the system even when startingfrom a random initialization, and are capable ofhandling non-stationary, large-scale instances.

Figure 1. Histogram of gradient estimators for average queue
length of a linear softmax policy in the N-model network (seeAppendix A.2), across a horizon of N = 1000 steps.

ically allocate server capacity to jobs in order to optimize
various performance metrics, such as minimizing average
waiting time of jobs in the queue, average length of queues,
maximizing job-server compatibility, etc.
Solving this control problem is known to be difficult, due
to high stochasticity, nonlinear and non-smooth dynamics,
non-stationarity, and complexity of the network topology.
While there are many effective heuristics such as the cµ-rule(Mandelbaum & Stolyar, 2004), the MaxWeight policy (Dai& Lin, 2005), or fluid policies (Chen & Yao, 1993), theirperformance can vary across different problem settings, andtheoretical guarantees are only available for specific network architectures and specific load regimes. For thesereasons, recently there has been much interest in Reinforcement Learning (RL) as a way to develop control policies
for generic queuing networks. While there has been muchdevelopment in applying model-free RL, these approaches
face challenges due to high variance of trajectories, large
state and action spaces, as well as instability of the system
under poor controllers (Meyn, 2008; Dai & Gluzman, 2022;
Liu et al., 2022; Zaki et al., 2021).

1. Introduction
Queuing networks are ubiquitous for modeling jobprocessing systems, which arise in applications related to
manufacturing systems (Perkins & Kumar, 1989), call centers (Bassamboo et al., 2006), communications networks(Maguluri et al., 2012), hospitals (Dai & Shi, 2019), andmore. In these systems, jobs arrive to queues at random
interarrival epochs and are either routed to a server to beprocessed, or wait in the queue if no servers are available.
After a random processing time, they either leave the system or are routed to further servers for further processing.
The optimal control problem is determining how to dynam*

Equal contribution 1 Decision, Risk and Operations, Columbia
Business School, New York, US. Correspondence to: Ethan Che<eche25@gsb.columbia.edu>.

In this work, we provide a model-based framework for policy optimization and planning for controlling multi-class,multi-pool queuing networks through the lens of discrete
event simulation. Rather than standard RL approaches forqueuing, which typically rely on a discrete-time Markov De-

Published at the Differentiable Almost Everything Workshop of the
40 th International Conference on Machine Learning, Honolulu,
Hawaii, USA. July 2023. Copyright 2023 by the author(s).

1Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

cision Process (MDP) representation which only holds whenarrivals and service times are stationary and are exponentially distributed, we model the system as a discrete event
dynamical system (Ho, 1987; Ho & Cao, 2012; Tsitsiklis,
1989) in continuous time, sampling the system wheneveran arrival or service occurs. This model only relies on
samples of event times and allows for non-stationary andnon-Markovian arrival and service processes, which better
matches behavior of real-world systems (Green et al., 2007).
Importantly, this modeling framework provides a novel way
to obtain approximate first-order derivatives of performence
metrics with respect to controls applied to the system. First,
we use a continuous relaxation of the action space motivated by fluid models to obtain a ‘reparameterization trick’(Kingma & Welling, 2013) of the dynamics with respect tothe actions. While these dynamics are non-differentiable,
we apply the ‘straight-through’ trick (Bengio et al., 2013)
to obtain gradients. This way, we can directly compute
the gradient of rewards or costs accumulated in a rollout
through reverse-mode autodifferentiation. This can be donein parallel with minibatches of rollouts using GPUs. Whilestraight-through gradients are known to be potentially biased, we observe the bias is very small in practice even forlong episodes, and the variance is significantly lower thanstandard zero-order estimators (see Figure 1 and Figure 2).

Figure 2. Normalized RMSE and bias of gradient estimators foraverage queue length of a linear softmax policy in the N-model
network (see Appendix A.2) across a horizons ranging from N =50 to N = 1000 steps. Normalized by the magnitude of theaverage gradients.

τ A (t) ∈ Rn+ , which may be depend on time t, and if they
are not served immediately they wait in queue.
2.1. Continuous action relaxation
Upon observing the state, the controller chooses an actionm×na ∈ P := {X ∈ R+
: X1 = 1, 1⊤ X = 1}, which
is an assignment of job types to servers, which must be
feasible according to M . If server i is assigned to job typej, then a job from queue j is routed first-in first-out (FIFO)
to the server, upon which the server begins processing it.
The processing time for each job-server pair is describedm×nby τ S (t) ∈ R+
, which may different across job-server
pairs and can vary across time. Yet, since the servers only
process jobs routed to them according to a, the service timeS
for server i and a job of type j is equal to τij
only if aij = 1,
and infinite otherwise.

We propose two RL algorithms based on these first-order
gradients. Straight-through Policy Gradient (StraightPG)
updates policy parameters using straight-through gradientsof cost computed via autodifferentiation along mini-batchesof sample paths. Straight-through Planning is a modelpredictive control style policy that solves a planning problem via straight-through gradients. Our preliminary results
show that despite relying on approximate gradients, thesepolicies match the performance of well-established heuristics and model-free RL algorithms on standard instances,
but with improved sample efficiency. And unlike other RL
methods, which require switching to a stabilizing policy(Liu et al., 2022; Zaki et al., 2021) or behavorial cloningof a stabilizing policy (Dai & Gluzman, 2022), the training
process is robust to instability even when starting from arandom initialization.

Drawing inspiration from fluid models (Chen & Yao, 1993),
we consider a novel continuous relaxation of the dynamics,
in which the realized service times are equal to τ S /a, withdivision being element-wise. This allows for non-integral
assignments, which represent a server splitting its capacity between several jobs or a job being served by several
servers. This way, we consider deterministic policies overa continuous relaxation of the action space, which departs
from standard RL treatments of queuing control that typically use stochastic policies over discrete actions. Note that
we use this relaxation for training and optimization, but forevaluation we restrict to integral assignments. After a job isserved it leaves the system or enters another queue.

2. ModelWe consider a continuous-time multi-class, multi-pool queuing network with m servers and n queues corresponding to
n different job types, with time indexed by t ∈ Rn+ . The
topology of the queuing network M = {0, 1}m×n describes
the compatibility of servers and jobs, with Mij = 1 if server
i can serve job type j. The job lengths (or state) x(t) ∈ Nn
are the number of jobs of each type still remaining in thesystem, including jobs waiting in queue and jobs in service. Jobs arrive according to stochastic interarrival times

2.2. Discrete event dynamical system
Although the process evolves in continuous time, we canobtain a discrete-time characterization by sampling the process when an arrival or a service occurs. We refer to this2

Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

as the discrete-event dynamical system model. Let τk ∈ R
for k ∈ N+ be denote as the kth event epoch. We define
xk := x(τk ) to be the job length process sampled at τk , andlet τkA , τkS be the arrival times and service times at epoch k.
Starting from τ0 = 0, the next event epoch is determinedas the minimum of the all the total possible event times:τk+1 = min{τkA , τkS /a}.

Algorithm 1 Straight-through Policy Gradient (StraightPG)
Input: Initial state x0 , learning rate α, epochs m, episodelength N , initial policy ρθ0for i = 0 to m − 1 do
Draw a sample
PN path and compute cumulative cost
V (x0 ) = k=1 c(xk , a(ρθi (xk )))(τk+1 − τk ).
Compute ∇θi V (x0 ) with straight-through trick.
Update θi+1 ← θi − α∇θi V (x0 ).
end forDefine e(τkA , τkS /a) ∈ {0, 1}2n as the one-hot vector denoting which event had the minimum time. The first n
entries refer to an arrival for type j, the latter n entries referto service of type j. The model is also described by anevent matrix ∆ ∈ R2n×n , where each row i corresponds to
how the job lengths x(t) are updated if the ith event is theminimum event time. With this, the dynamics are:
τk+1 = min{τkA , τkS /a}xk+1 = xk + ∆⊤e(τkA , τkS /a)action a. To do so, we employ the straight-through trick
to differentiate through e, a trick which has been used effectively in discrete representation learning (Van Den Oordet al., 2017) and model-based RL (Hafner et al., 2020). In
other words, we replace e with

(1)
(2)

ê(τkA , τkS /a) := e(τkA , τkS /a) + softmin(τkA , τkS /a)(3)− sg[softmin(τkA , τkS /a)]where the first equation describes the event epoch and thesecond equation is the update for the state. The arrival andservice times of other jobs are reduced by τk+1 and are reset
if the event was the minimum time event.

where sg[·] is the stop gradient operator. This will produce
identical trajectories as in the update (1) in the forward pass,
but behaves as a softmin in the backward pass.

Example 2.1. Consider a single server M/M/1 queue witharrival rate λ and service rate µ. Interarrivals are τA ∼Exp(1/λ), services are τS ∼ Exp(1/µ), a = 1{x > 0},( (1, 0) τ A < τ S /a
1
⊤
A S
∆ =
, e(τ , τ /a) =−1(0, 1) τ A > τ S /a

As a result, we can estimate gradients of the cost of a trajetory with respect to any routing
PNaction or sequence of actions taken along the way: ∇ak k=1 c(xk , ak )(τk+1 −τk ).
Figure 1 and 2 compare the variance and bias of the straightthrough estimator with finite differences (using the continuous relaxation and under coupled random seeds) and the REINFORCE estimator for a simple linear policy. We observethat the estimator achieves significant variance reduction
over alternatives, with low bias.

This can be seen as a ‘reparameterization trick’ (Kingma& Welling, 2013) of the dynamics with respect to actions;enabling us to evaluate different actions under the same
sample path of stochastic inputs, i.e. τ A and τ S .

2.5. Sinkhorn parameterization

2.3. Control problem

Instead of optimizing over actions a ∈ P , we parameterize actions in terms of priorities ρ ∈ Rm×n which
produce assignments through linear assignment: ā(ρ) =maxa∈P Tr(aρ).

The objective of the controller is to minimize an average
cost metric, such as the average job lengths or the average time the system is busy. Let c(x, a) be the instantaneous cost given job lengths x and assignment a. A typical cost is a c(x, a) = h⊤ x where h is a vector of holding costs associated with each job type. The objective is
to minimize the average cost over a horizon of length T :RT
min 0 c(x(t), a(t))dt. In the discrete-event dynamical system, we replace thisthe average of the cost over the
PNwithTc(xk , ak )(τk+1 − τk ). where NT
event times: min k=1
is the number of events that occur before time T .

Example 2.2. The cµ-rule uses a static priority ρij =hj µij where µij is service rate of server i and job type j.
MaxWeight uses a state-dependent priority ρij = xj hj µij .
To maintain differentiability, we replace the hard linear assignment with Sinkhorn’s algorithm: a(ρ) =maxa∈P Tr(aρ) + λKL(a|11⊤ ). Thus, actions only specify
relative values of assigning a server to a job, rather thanspecifying the assignment exactly.

2.4. Straight-through gradients

3. Algorithms and Empirical Results

The key source of non-differentiability in the model ise(τkA , τkS /a), which is a onehot argmin of the event times.
If this could replaced with a differentiable surrogate, thejob-lengths would evolve as a differentiable function of theUsing the first-order gradients obtained through the straightthrough trick and Sinkhorn parameterization, we propose
two policies. First, Straight-through Policy Gradient, which
3

Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

Algorithm 2 Straight-through Planning (StraightPlan)
Input: Initial state x0 , total episode length N , planninghorizon H.
for k = 0 to N do
Solve planning problem minρk:k+H V (xk )=Pk+H
i=k c(xi , a(ρi ))(τk+1 − τk ) via gradient descenton sample-path rollouts.
Play action ak = ā(ρk ) via linear assignment.
Update states and times (xk+1 , τk+1 ) with (1).
end forTable 2. Average Queue Length in Reentrant Network
N UM C LASSES

RFP

PPO

S TRT PG

6
9
12
15
18
21

15.42
24.92
36.86
43.63
52.98
59.05

14.13
23.27
32.17
39.30
51.47
55.12

15.68
22.43
32.38
37.51
48.36
55.45

Table 3. Average Holding Cost in Parallel Server System
P OLICY

parameterizes priorities ρθ that produce actions through
Sinkhorn’s algorithm. The parameter θ is updated by gradients of the cost from sample rollouts directly through auto
differentiation. Second, Straight-through Planning, which isa model-predictive control (MPC) style policy that solves aplanning problem via gradients of sample-path rollouts. We
evaluate the performance of the policy on a few examples.

cµ- RULEF LUIDM AX W EIGHTS TRAIGHT PG
S TRAIGHT P LAN

To further evaluate the quality of the straight-through gradients, we train policy gradient to minimize the stationary
average queue-length of the criss-cross network (Harrison &Wein, 1990), which is a simple tandem queue with 2 servers
and 3 job types. We evaluate across different load levels(see Table 4 in Appendix A.3). To assess validity for large
systems, we train the method on a class of reentrant queuing networks introduced in (Bertsimas et al., 2014) acrossdifferent system sizes. We compare against benchmarks
such as the robust fluid policy (RFP) (Bertsimas et al., 2014)
and Proximal Policy Optimization (PPO) as implementedin (Dai & Gluzman, 2022), which uses several variance reduction strategies. We observe in Table 1 and 2 that despite
potential bias in the gradients, straight-through policy gradient outperforms RFP and achieves similar performance
to PPO, typically using fewer policy iterations and without
requiring an initial behavior cloning step to ensure stability.

4. Related Work
Our work directly builds on approaches for queuing in RL
as in (Zaki et al., 2021; Dai & Gluzman, 2022; Liu et al.,2022), and proposes a new gradient estimator. Our approachis related to differentiable simulation for policy optimization (Suh et al., 2022; Mora et al., 2021), as it uses exact
gradients to update policies. Our work builds on previous
literature on sensitivity analysis in discrete event systems(Ho & Cao, 2012; Ho, 1987), but methodologically approaches this from the viewpoint of gradient estimation for
discrete random variables (Tucker et al., 2017; Mohamed
et al., 2020). Of particular relevance is (Arya et al., 2022),
who propose a framework for autodifferentation (AD) withdiscrete random variables, which produces unbiased gradient estimates with stochastic coupling for variance reduction.
This offers a promising method for queuing network control.
Our approach instead constructs a differentiable ‘reparameterization trick’ of the dynamics with respect to actions forreverse-mode AD, appropriate for computing gradients ofscalar costs with respect to high-dimensional policy parameters and actions. While their framework also allows forreverse-mode AD, gradients are no longer unbiased due tonon-linear dynamics. Our estimator has low bias in practiceand significantly reduced variance than the finite-difference
estimator with coupled random seeds (in Figure 1 and 2),
which suggests that embedding discrete actions into a continuous space and reparameterizing the dynamics may obtaingreater variance reduction than stochastic coupling alone.

3.2. Parallel-server system
To test the applicability of the model to non-stationary settings, we consider a parallel-server system with 5 servers
and 5 queues (Chen et al., 2021). We consider a non-Table 1. Average Queue Length in Criss Cross Network

I MB .L IGHT
BAL .L IGHT
I MB .M ED
BAL .M ED
I MB .H EAVY
BAL .H EAVYO PTIMAL

RFP

PPOS TRT PG

0.671
0.843
2.084
2.829
9.97
15.23

0.677
0.855
2.133
2.920
10.10
15.58

0.671
0.844
2.084
2.833
10.01
15.35

0.673
0.854
2.107
2.888
10.01
15.34

107.01
83.35
88.34
77.29
78.25

stationary arrival process where queues 1 and 3 experiencea demand surge until time t = 100. We compare against
cµ-rule, MaxWeight, and the fluid policy in Table 3. We
observe that the gradient-based algorithms outperform otherpolicies and do so by planning ahead and tapering service
to the overloaded queues before the demand surge ends.

3.1. Criss-cross Network and Reentrant Network

S ETTING

AVG H OLDING C OST

4

Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

References

Ho, Y.-C. L. and Cao, X.-R. Perturbation analysis of discrete event dynamic systems, volume 145. Springer Science & Business Media, 2012.

Arya, G., Schauer, M., Schäfer, F., and Rackauckas, C.
Automatic differentiation of programs with discrete randomness. Advances in Neural Information Processing
Systems, 35:10435–10447, 2022.

Kingma, D. P. and Welling, M. Auto-encoding variational
bayes. arXiv preprint arXiv:1312.6114, 2013.

Bassamboo, A., Harrison, J. M., and Zeevi, A. Design andcontrol of a large call center: Asymptotic analysis of an
lp-based method. Operations Research, 54(3):419–435,
2006.

Liu, B., Xie, Q., and Modiano, E. Rl-qn: A reinforcement
learning framework for optimal control of queueing systems. ACM Transactions on Modeling and Performance
Evaluation of Computing Systems, 7(1):1–35, 2022.

Bengio, Y., Léonard, N., and Courville, A. Estimating or
propagating gradients through stochastic neurons for conditional computation. arXiv preprint arXiv:1308.3432,
2013.

Maguluri, S. T., Srikant, R., and Ying, L. Stochastic models
of load balancing and scheduling in cloud computing
clusters. In 2012 Proceedings IEEE Infocom, pp. 702–
710. IEEE, 2012.

Bertsimas, D., Nasrabadi, E., and Paschalidis, I. C. Robust fluid processing networks. IEEE Transactions on
Automatic Control, 60(3):715–728, 2014.

Mandelbaum, A. and Stolyar, A. L. Scheduling flexible
servers with convex delay costs: Heavy-traffic optimalityof the generalized cµ-rule. Operations Research, 52(6):836–855, 2004.

Chen, H. and Yao, D. D. Dynamic scheduling of a multiclassfluid network. Operations Research, 41(6):1104–1115,
1993.

Meyn, S. Control techniques for complex networks. Cambridge University Press, 2008.

Chen, J., Dong, J., and Shi, P. Optimal routing under demand surges: The value of future arrival rates. Available
at SSRN 3980227, 2021.

Mohamed, S., Rosca, M., Figurnov, M., and Mnih, A.
Monte carlo gradient estimation in machine learning. The
Journal of Machine Learning Research, 21(1):5183–5244,
2020.

Dai, J. G. and Gluzman, M. Queueing network controls via
deep reinforcement learning. Stochastic Systems, 12(1):30–67, 2022.

Mora, M. A. Z., Peychev, M., Ha, S., Vechev, M., and Coros,
S. Pods: Policy optimization via differentiable simulation.
In International Conference on Machine Learning, pp.
7805–7817. PMLR, 2021.

Dai, J. G. and Lin, W. Maximum pressure policies instochastic processing networks. Operations Research,
53(2):197–218, 2005.

Perkins, J. R. and Kumar, P.
Stable, distributed,
real-time scheduling of flexible manufacturing/assembly/diassembly systems. IEEE Transactions on
Automatic Control, 34(2):139–148, 1989.

Dai, J. G. and Shi, P. Inpatient overflow: An approximate dynamic programming approach. Manufacturing & Service
Operations Management, 21(4):894–911, 2019.

Suh, H. J., Simchowitz, M., Zhang, K., and Tedrake, R.
Do differentiable simulators give better policy gradients?
In International Conference on Machine Learning, pp.
20668–20696. PMLR, 2022.

Green, L. V., Kolesar, P. J., and Whitt, W. Coping with timevarying demand when setting staffing requirements for aservice system. Production and Operations Management,
16(1):13–39, 2007.

Tsitsiklis, J. N. On the control of discrete-event dynamicalsystems. Mathematics of Control, Signals and Systems, 2(2):95–107, 1989.

Hafner, D., Lillicrap, T., Norouzi, M., and Ba, J. Mastering atari with discrete world models. arXiv preprint
arXiv:2010.02193, 2020.
Harrison, J. M. and Wein, L. M. Scheduling networks ofqueues: Heavy traffic analysis of a two-station closed
network. Operations research, 38(6):1052–1064, 1990.

Tucker, G., Mnih, A., Maddison, C. J., Lawson, J., and SohlDickstein, J. Rebar: Low-variance, unbiased gradientestimates for discrete latent variable models. Advancesin Neural Information Processing Systems, 30, 2017.

Ho, Y.-C. Performance evaluation and perturbation analysis
of discrete event dynamic systems. IEEE Transactions
on Automatic Control, 32(7):563–572, 1987.

Van Den Oord, A., Vinyals, O., et al. Neural discrete representation learning. Advances in neural information
processing systems, 30, 2017.
5

Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

Zaki, M., Mohan, A., Gopalan, A., and Mannor, S. Improper reinforcement learning with gradient-based policy
optimization. arXiv preprint arXiv:2102.08201, 2021.

6

Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

A. Appendix A
A.1. Implementation
We detail the Straight-through Policy Gradient algorithm and the Straight-through planning policy in Algorithms 1 and 2
respectively.
For the criss-cross network and reentrant networks, we parameterize the policy ρθ as a multi-layer perceptron (MLP) with
input dimension equal to the number of queues, output dimension m × n and 3 hidden layers with 128 hidden units ineach layer. For the criss-cross network, each policy gradient step uses a mini-batch of 50 episodes in parallel of length
N = 5, 000 for light loaded, N = 10, 000 for medium loaded systems, and N = 20, 000 steps for heavily loaded system.
For the reentrant network each policy gradient step uses a mini-batch of 50 episodes in parallel of length N = 20, 000 steps.
For the parallel-server system, due to non-stationary, we also input time t as another input variable, with all else being thesame.
A.2. N-model
The N-model is a network with two queues and two servers. Each queue has it’s own designated server, but server 1 can also
capable of serving jobs in queue 2, at a reduced service rate. The network topology is:


1 1
M=
0 1
In Figure 1 and Figure 2, we consider a system where interarrivals of the queues are distributed Exp(1/λ1 ) and Exp(1/λ2 )with λ1 = 0.4ρ and λ2 = 1.3ρ where ρ = 0.95. Service times of server i and queue j are distributed Exp(1/µij ) withservice rates µij :


1 0.5
µ=
0 1
To assess the quality of gradients, we consider a linear softmax policy of the form:θ11 x1 θ12 x2
a(x) = softmax0θ22 x2where the softmax is performed row-wise. We take gradients of the cumulative holding cost with respect to the parameters θ,
evaluated at θ11 = θ12 = θ22 = 2 across an episode of N events:#"NX
xk,1 + xk,2∇θ Vθ (x0 ) = ∇θ
k=1

For the straight-through and finite difference estimators, we calculate the gradient with respect to the non-integral assignmenta(x). For the finite difference estimator, we use the two point estimator with σ = 0.1:
∇θ Vθ (x0 ) =1
E [Vθ+σZ (x0 ) − Vθ (x0 )]σwhere Z ∼ N (0, I) ∈ R3 is a standard normal random vector. For both rollouts, we use the same random seed, which
couples the arrival and service times. Note that the finite difference estimator is using the same continuous action relaxation
as the straight-through estimator and produces an unbiased estimate of ∇θ Vθ (x0 ). For this reason, we use this estimator tomeasure the bias of the straight-through gradients.
To compare the REINFORCE estimator, which typically considers stochastic policies over discrete actions, we look at apolicy that instead of choosing the action a(x), samples from the softmax to produce a hard assignment. Since the evaluatedpolicy is different, the gradient differs from what we obtain from the straight through and finite difference estimators. While
the sign of the REINFORCE gradients are the same, they tend to be off by a scale factor. For a comparison, in Figure 2 we
present the normalized root mean squared error (RMSE) of the different gradient estimators, where the RMSE is normalized
by the magnitude of the gradient. For the straight-through estimator, RMSE is calculated by taking the mean squared error
with the average finite difference estimator (across 100,000 estimates). For the other estimators, RMSE is the square root ofsum of variances of each component, considering both estimators to be unbiased. In Figure 1, we rescale the gradient so thatits mean is equal to that of the finite-difference estimator.
7

Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

Table 4. Load parameters in each settings
S ETTING

λ1

λ3

µ11

µ13

µ22

I MB .L IGHT.
BAL .L IGHT.
I MB .M ED .
BAL .M ED .
I MB .H EAVY.
BAL .H EAVY.

0.3
0.3
0.6
0.6
0.9
0.9

0.3
0.3
0.6
0.6
0.9
0.9

2
2
2
2
2
2

1.5
1
1.5
1
1.5
1

2
2
2
2
2
2

A.3. Criss-cross
The criss-cross network (Harrison & Wein, 1990) is a network with 3 queues and 2 servers. Queues 1 and 3 can be served byserver 1. When a queue 1 job is served, it joins queue 2 which is served by server 2. When jobs in queue 2 and 3 are served,
they leave the system. The network topology is:

M=



1 0 1
0 1 0



We consider a system where interarrivals of the queues are distributed Exp(1/λ1 ), Exp(1/λ2 ) with service times of server iand queue j are distributed Exp(1/µij ) with service rates µij . Due to the tandem nature of the queue, the event matrix ∆ is:

1
0
0
0
0
0


0
0
1


∆=
0
−1 1

 0 −1 0 
0
0 −1


where the first 3 rows describe the state change when an arrival occurs (note that no arrivals to queue 2 occur unless theyleave from queue 1) and the second 3 rows discribe the state change when a service occurs. We consider different parametersettings for the load on the system, described in Table 4.
A.4. Reentrant network
The reentrant queuing network is parameterized by L, which is the number of layers. There are L servers and 3L many jobclasses/queues. In each layer, there are 3 queues and each can be served by a single server, who only serves those 3 queues.
Job types have an index j ranging from 0, ..., 3L − 1. After a job of type j is served, the job is sent to queue i + 3. F or jobsof type 3L − 3, after the job is processed it is sent to job type 1. For jobs of type 3L − 2 and 3L − 1 after service they are
immediately processed.
Arrivals only occur to queues 0 and 2. The arrival times for both queues are distributed Exp(140/9). Servers are indexed byi ranging from 0, ..., L. Service times are exponential Exp(1/µij ) with


1/8





1/2



1/4
µij =

1/6




1/7



1

i mod 2 = 0 and j mod 3 = 0
i mod 2 = 0 and j mod 3 = 1
i mod 2 = 0 and j mod 3 = 2
i mod 2 = 1 and j mod 3 = 0
i mod 2 = 1 and j mod 3 = 1
i mod 2 = 1 and j mod 3 = 2
8

Dynamic Control of Queuing Networks via Differentiable Discrete-Event Simulation

A.5. Parallel-server
We consider a parallel-server system with 5 servers and 5 queues. The network topology is:
1 0 0 0 0
1 1 0 0 0



M =
1 0 1 0 0
0 0 1 1 0
0 0 0 1 1
The vector of holding costs on the job lengths for each queue is [2, 1, 1, 2, 1].

We consider a system where interarrivals of the queues are non-stationary and distrubted Exp(1/λi (t)), described as follows(
2.4 t ≤ 100
λ1 =0.4 t > 100(
0.6 t ≤ 100
λ2 =0.8 t > 100
λ3 = 0.8
(
1.6
λ4 =
0.8

t ≤ 100
t > 100

λ5 = 0.6
The service times are distributed Exp(1/µij ) where
µij =

(

Mij
0.8 · Miji=ji ̸= jAfter being served, the job leaves the system. We consider a horizon length of N = 3000 events, which is around T = 365
in time units.

9