# Transfer learning applications for anomaly detection in wind turbines

Cyriana M.A. Roelofs<EMAIL_ADDRESS>Christian Gück
<EMAIL_ADDRESS>Stefan Faulstich
<EMAIL_ADDRESS>Fraunhofer IEE, Joseph-Beuys-Straße 8,
34117 Kassel, Germany

###### Abstract

Anomaly detection in wind turbines typically involves using normal behaviour
models to detect faults early. However, training autoencoder models for each
turbine is time-consuming and resource intensive. Thus, transfer learning
becomes essential for wind turbines with limited data or applications with
limited computational resources. This study examines how cross-turbine
transfer learning can be applied to autoencoder-based anomaly detection. Here,
autoencoders are combined with constant thresholds for the reconstruction
error to determine if input data contains an anomaly. The models are initially
trained on one year’s worth of data from one or more source wind turbines.
They are then fine-tuned using smaller amounts of data from another turbine.
Three methods for fine-tuning are investigated: adjusting the entire
autoencoder, only the decoder, or only the threshold of the model. The
performance of the transfer learning models is compared to baseline models
that were trained on one year’s worth of data from the target wind turbine.
The results of the tests conducted in this study indicate that models trained
on data of multiple wind turbines do not improve the anomaly detection
capability compared to models trained on data of one source wind turbine. In
addition, modifying the model’s threshold can lead to comparable or even
superior performance compared to the baseline, whereas fine-tuning the decoder
or autoencoder further enhances the models’ performance.

###### keywords:

condition monitoring , transfer learning , autoencoder , anomaly detection ,
wind turbines , fault detection

††journal: Energy & AI

##  1 Introduction

††Nomenclature:

TL

    Transfer Learning
AE

    autoencoder
ANN

    artificial neural network
ARCANA

    Anomaly Root Cause Analysis
ML

    machine learning
MSE

    mean square error
RMSE

    root mean square error
NBM

    normal behaviour model
O&M

    Operation and maintenance
OP-mode

    operational mode
RE

    AE reconstruction error
SCADA

    Supervisory Control and Data Acquisition

Wind power plays a crucial role in addressing climate change, offering
significant promise among renewable energy sources. However, wind farms are
often situated in regions with challenging weather conditions, leading to
complex operating conditions and potential failures in wind turbines. This
results in increased maintenance costs and decreased efficiency. To mitigate
these challenges, there has been a growing focus on condition monitoring and
early fault detection using machine learning models for wind turbines [1].
Numerous approaches concentrate on training a normal behaviour model (NBM)
using wind turbine Supervisory Control and Data Acquisition (SCADA) data, to
identify anomalous behaviour or detect specific failures in wind turbines.

Training an anomaly detection model based on NBM requires enough data
representing normal operation of the wind turbine. However, such data is not
available for new installations or after maintenance actions, such as an oil
refill or component replacements resulting in changed normal behaviour. To
overcome the scarcity of training data, the application of Transfer Learning
(TL) can be beneficial. TL involves transferring learned behaviour from one
model to another, similar to how humans learn new tasks by applying knowledge
from related tasks. In the context of wind turbines, this involves using a
pretrained model for anomaly detection in one wind turbine (or a group of wind
turbines) and repurposing it for anomaly detection in another wind turbine.
This approach saves significant time and resources that would otherwise be
required for training a new model [2, 3].

In this study, our focus is on exploring the applications of TL for
autoencoder-based anomaly detection in wind turbine data. We first explore TL
within one wind farm, to test whether anomaly detection is still possible for
an wind turbine for which only a small amount of training data is available.
Specifically, a model is trained on data of one wind turbine and transferred
to another wind turbine within the same wind farm. We examine the impact of
different amounts of tuning data for different wind turbine types, using
various TL methods. Additionally, we use data of multiple wind turbines to
train a NBM and fine-tune it on data of a single target wind turbine which is
not in the training dataset. The hypothesis is that this model can generalise
better and be better tuned to the target task. In both cases, we validate the
transferred models’ ability to detect anomalies by evaluating three case
studies, demonstrating that early fault detection in wind turbines is indeed
possible even with a limited amount of training data.

The paper is structured as follows. First we present an overview of related
work in section 2. In section 3 we describe the data and models used and the
TL methods which we applied. Next, the results are presented and discussed in
4. The results of the transferred models are first compared to their
corresponding baseline models, to show they are capable of distinguishing
normal from anomalous behaviour almost or as well as the baseline model.
Afterwards, we demonstrate the model ability to detect failures early using
three case studies. Finally, the paper is summarised in section 5.

##  2 Related work

Research on TL applied to unsupervised or semi-supervised (only using
‘healthy’ data) machine learning, such as anomaly detection, is limited [4,
5]. Most work focuses on dimensionality reduction and clustering, which is
often not applicable to anomaly detection in predictive maintenance settings.
Due to limited occurrences of failures, anomaly detection models are often
hard to train. Here, Transfer Learning approaches may help overcome this
problem [6]. In addition, the training of wind turbine-specific models is not
always feasible. For example, if a wind turbine is newly installed or not all
historical data is informative due to re-calibration or component
replacement(s). Also, training ML models on only one turbine greatly reduces
the amount of training time for the whole wind farm (containing up to hundreds
of turbines) [5].

Zgraggen et al. present three TL approaches for anomaly detection in wind
turbines using deep learning methods [5]. The anomaly detection model used is
a normal behaviour model (NBM) of a regression CNN that has 4 SCADA parameters
as input and predicts the gearbox bearing temperature. Three TL methods are
tested for two scenarios: cross-turbine TL within the same wind farm and
cross-turbine TL across different wind farms. For the first scenario the
simplest method, where an additional linear regression (LR) is added to the
pre-trained CNN, works almost as well as CNN trained from scratch on the
target turbine. For the cross-wind-farm case, there are more non-linear domain
shifts to consider as well. Here, they added both the LR model as well as a
CNN model that predicts the error between the input and the output of the LR.
Simply fine-tuning the original CNN after pretraining did not yield better
results, as it can easily overfit.

Another application of TL for anomaly detection in wind turbines is presented
by Dimitrov et al. [7]. The authors use TL to create a physics-informed
machine-learning model, by transferring knowledge from a model trained on
simulations to a target model trained on real data. The two models tested are
a five-layer dense neural network (ANN) and a five-layer denoising autoencoder
(DAE). Both use six SCADA parameters, of which the electrical power is used as
a health indicator. The prediction performance of the ANN increases
significantly with TL, even when only one month of SCADA data is used for
fine-tuning. The power prediction of the DAE did not improve a lot, as it
already had very high performance. The anomaly detection ability of the DAE
was however not tested.

To enable the use of structurally incomplete data (where only some sensor
values are available due to different sample rates), Chen et al. use TL to
train fault classification models for different types of incomplete data and a
model based on only the complete data [6]. The information is then transferred
between these models to optimally use all available data, instead of only a
subset (the structurally complete data). The proposed method has a higher
accuracy than the model solely trained on structurally complete data and the
resulting models can be applied to all data.

Wen et al. use TL to transfer information from univariate time series models
to multivariate time series models for anomaly detection [8]. The model is
based on an encoder-decoder network (U-Net) and extended for multivariate time
series. The model for multivariate time series has encoders for each input
feature, to which the information from the models for univariate time series
is transferred. The resulting TLmodels all outperformed the models without TL.

Recently, Simani et al. study TL for anomaly detection in SCADA data of wind
turbines [9]. They focus improving the accuracy of the NBM’s output signals
and show that this indeed is possible between different wind turbines in the
same wind farm. However, the application to fault diagnosis or anomaly
detection is not investigated.

##  3 Method

In this section, the data, model and transfer learning methods are described.

###  3.1 Data

This study utilises data from two offshore wind farms, WF A and WF B, located
in the German North Sea. These wind farms feature different types of turbines.
The dataset used in this study is derived from SCADA data and spans several
years. It consists of ten-minute-average values from 240 sensors for WF A and
30 sensors for WF B. The sensors capture various measurements, including power
production, wind speed and direction, as well as temperature measurements of
bearings.

In addition to the SCADA data, the dataset includes information about the
operational mode (OP-mode) for each timestamp. Examples of OP-modes include
normal operation, derated operation, service, downtime, and others. For the
purposes of this study, only ’normal operation’ is considered as normal
behaviour, while all other OP-modes are categorised as ’anomalous’.

###  3.2 Anomaly detection model

For the anomaly detection model, we use a similar setup as presented in [10,
11].

First data representing normal behaviour is selected based on the OP-mode. In
addition, only data with power $P>0$. Due to communication issues the OP-mode
is not always properly closed, resulting in some downtimes being missed and
logged as, for example, ’normal operation’. Furthermore, we filter out any
values where the power is high but the wind speed is outside the cut-in and
cut-off values of the wind turbine. The data filtered out is labeled as
’anomalous’.

Next, data is prepared for the autoencoder (AE) model. First, we remove any
features that represent counter values (such as the number of cable windings
or the energy production since commissioning of the wind turbine) or set
points. Then, we eliminate features that are constant or have up to 3 unique
values (indicating low variance). The data does not contain missing values, as
they are replaced with zero values. Consequently, we also remove features that
have more than 80% zero values. For features representing angles, we transform
them into sine and cosine values to ensure their continuity. Finally, we scale
the data between 0 and 1 using min-max-scaling.

The AE model used is a fully connected feed forward neural network with a
small bottleneck layer and employs the PReLU [12] activation function. For WF
A the network consists of 7 layers, where the hidden layers have 200, 100, 50,
100 and 200 units respectively. For WF B, there are 3 hidden layers, with 25,
10, and 25 units. The network’s parameters are optimised using the Adam
optimisation algorithm [13] with a learning rate of 0.001.

To determine whether the model detected an anomaly, we calculate the root mean
square error (RMSE) of the reconstruction errors for each timestamp in the
training data. This score is referred to as the Anomaly Score $S$. Next, we
set a threshold for the training data to maximise the $F_{1/2}$-score. The
anomaly labels used for optimisation are based on the OP-mode and power values
as described before for normal data selection. We use the $F_{1/2}$-score
instead of the $F_{1}$-score because we aim to prioritise precision over
recall for our model.

###  3.3 Transfer learning

In this paper, we use two types of source models: models trained on 12 months
of data of one wind turbine and models trained on 12 months of data of
multiple wind turbines. These models are then transferred to the target wind
turbine using the following three methods:

  * 1.

Threshold: the AE of the source model is applied as is. Only the anomaly score
threshold is updated.

  * 2.

Decoder: In addition to updating the threshold, the decoder of the AE is fine-
tuned while keeping the encoder fixed. We use a learning rate of 0.001 and
train for 10 epochs.

  * 3.

AE: In this case the complete AE network is fine-tuned using a learning rate
of 0.0001 (1/10 of original learning rate) and 10 epochs. Afterwards, the
threshold is updated.

From this point on, we will refer to the two types of TL models as asset-to-
asset and multi-asset models.

###  3.4 Evaluation

The anomaly detection capability of the model is evaluated using the
$F_{1/2}$-score on one month of data immediately following the training or
tuning data. We determine whether anomalies were correctly detected by using
the labels ”anomaly” and ”normal” based on the OP-mode, wind, and power
measurements, as explained in Section 3.2.

We compare the performance of tuned models with a baseline model trained on 12
months of data. Furthermore, we conduct three case studies to demonstrate that
the transferred models are still effective in detecting anomalous behaviour
due to failures. Here, we use the criticality $C$ measure as defined in [11]
to evaluate the anomalies detected. The criticality value increases by one if
an anomaly is detected during normal operation, decreases by one if no anomaly
is detected during normal operation, and remains constant if the OP-mode is
not normal. The range of the criticality value is limited to [0, 1000].

##  4 Results and discussion

First, baseline models for 10 wind turbines of WF A and 6 wind turbines of WF
B were created to compare the TL models.

Both the baseline models as well as the source models were trained 1 year of
SCADA data representing normal behaviour. The average performance of these
baseline models by amount of training data is shown in figure 1. For WF A wind
turbines 1 and 40 are on opposite sides of the wind farm, whereas 28 and 29
are in the middle. For WF B wind turbines 42 and 70 are opposite sides, wind
turbines 55 and 56 are in the middle of the wind farm.

(a) WF A

(b) WF B

Figure 1: Performance of the baseline models.

Next, the source models were transferred using the methods described in
section. 3.3. In Section 4.1 and Section 4.2 the results for the asset-to-
asset and multi-asset models are described respectively. Finally, in Section
4.3 the case studies are described and analysed.

###  4.1 Asset to asset transfer learning

Figure 2: Transfer Learning results for WF A. The difference in $F_{1/2}$
score with respect to the baseline model for each wind turbine is shown.
Figure 3: Transfer Learning results for WF B. The difference in $F_{1/2}$
score with respect to the baseline model for each wind turbine is shown.

For WF A, the source models were based on data from wind turbines 1 and 40,
which are located on opposite sides of the wind farm. The purpose was to
investigate whether the location has an influence on the model’s performance
on the target wind turbines. We will refer to these wind turbines as ”source
wind turbines” from now on. For WF B, the same setup was used, using source
wind turbines 42 and 70.

Results are shown in Figure 2 and Figure 3 for TL models based on only a small
amount of tuning data, ranging from 1 to 3 months. Notably, these results show
that it is indeed possible to create an anomaly detection model with a
relatively small amount of data. The F1/2 scores of the TL models are
comparable or slightly better than the baseline models, which were trained on
a full year of SCADA data. Whether 1 month or 3 months of tuning data is used,
does not seem to have a significant effect.

Interestingly, the location of the chosen source wind turbines does not appear
to have a significant impact on the model’s performance. The TL models trained
using data from wind turbine 1 fit the data of wind turbine 40, on the
opposite side of the wind farm, quite well, even with just one month of tuning
data. The same is observed for WF B. There is however not enough data to state
that one can simply use any of the wind turbines as source, since the data
quality of the source wind turbine overshadows this effect. Therefore, we use
one of the middle assets in the WFs with a good baseline model performance as
source wind turbine for the case studies.

The transferred models show promising results, performing on par or slightly
better than the baseline models. It is worth noting that the autoencoder of
the source model is capable of reconstructing the data of the target wind
turbines, even without tuning the autoencoders. These models can yield good
model scores, even though the autoencoder was trained on data from a different
wind turbine, showing that the wind turbines within the same wind farm behave
very similarly, since the normal behaviour primarily depends on the
environment.

In some cases, the transferred models even outperform the baseline models. The
improved performance of the transferred models may be attributed to their
ability to fill knowledge gaps in the baseline models, particularly when the
training data of the baseline models lacks normal behaviour within a period
(e.g. winter or summer) similar to the prediction period.

Fine-tuning the decoder or the complete autoencoder further enhances the
models’ performance in terms of F1/2 score compared to the threshold method,
as depicted in Figures 2 and 3. Specifically, fine-tuning the decoder alone
outperforms the other TL methods in most cases. However, we find that caution
must be exercised when tuning the autoencoder, as it could lead to overfitting
and negative transfer.

###  4.2 Multi-asset model transfer learning

In this section, we discuss the results for the multi-asset TL models. Here,
the source models were trained on 9 wind turbines of WA and fine-tuned on the
target wind turbine, which is not part of the source model’s training data.
For WF B the source models were trained on data from 5 wind turbines and fine-
tuned on the target wind turbine. Results are shown in Figure 4.

(a) WF A

(b) WF B

Figure 4: Transfer Learning results using data of multiple assets to train a
source model and fine-tuning on the target wind turbine. The difference in
$F_{1/2}$ score with respect to the baseline model for each wind turbine is
shown.

Compared to the results of the asset-to-asset models, these models do not
perform as well, particularly for WF A. Since the wind turbines within one
wind farm behave very similar, this could be caused by overfitting, as the
model sees almost the same data multiple times. In case of WF A, the TL models
for all wind turbines except 40 perform significantly better with 3 months of
tuning data. For WF B the models’ performance does not seem to be
significantly influenced by the amount of tuning data.

###  4.3 Case Studies

For the case studies we use the same set up as before, where a baseline model
is trained on one year of data of the target wind turbine to compare with TL
models. The source models are trained with data of a wind turbine in the
middle of the wind farm, for which the baseline model performed well. For WF A
this is wind turbine 29, for WF B wind turbine 55. Next, the models’ decoders
were fine-tuned with 1 to 2 months of data of the target wind turbines.

The first two case studies are for wind turbines of WF A, the third case study
contains a failure in a wind turbine in WF B.

####  4.3.1 Case study 1

Anomalies are detected from 2022-12-22 onward, which were caused by problems
in the rotor brake and hydraulic system. The anomalies are caused by several
related sensor measurements suddenly behaving very differently. Both TL models
are able to detect this anomaly easily, as seen in Figure 5. The anomaly score
during abnormal operation on 2022-12-27 is lower for the TL models compared
the anomaly score of the baseline model, indicating that the TL models are
less sensitive than the baseline model. This is however expected, since only a
small amount of target data was used for tuning.

Figure 5: Case Study 1: Rotor brake and hydraulic problems for a wind turbine
of WF A cause a change point in some sensor measurements. The anomalous
behaviour is detected correctly by TL models using only 1 to 2 months of data
for fine-tuning. The upper graph depicts the anomaly score, threshold and
resulting criticality measure for the baseline model, the lower graph compares
the criticality for the three models.

####  4.3.2 Case study 2

Problems in the converter cooling system causes a sensor drift, which in turn
causes a positive trend in the anomaly score of the baseline model, s. Figure
6. After maintenance on 2022-10-07 the problem was fixed and the anomaly score
normalises. Both tuned models are able to detect the problem, although they
detect the anomaly later than the baseline model. Interestingly, the model
fine-tuned with 2 months of data is less sensitive to this behaviour and
detects the anomaly a week later. Nevertheless, the anomaly could be detected
weeks before maintenance takes place. Also, the anomaly score shows a positive
trend, before the threshold is exceeded.

Figure 6: Case Study 2: Sensor drift caused by problems in the water cooling
system of a wind turbine of WF A cause a change point in some sensor
measurements. The anomalous behaviour is detected correctly by TL models using
only 1 to 2 months of data for fine-tuning. The upper graph depicts the
anomaly score, threshold and resulting criticality measure for the baseline
model, the lower graph compares the criticality for the three models.

####  4.3.3 Case study 3

Anomalies are detected between 2023-02-28 and 2023-03-27 for a wind turbine in
WF B, as shown in Figure 7. The cause of the anomalies was a failure in the
gearbox of the wind turbine. Both TL model are able to detect some anomalous
behaviour on 2023-03-01, after which the anomaly score normalises for a short
time. Next, the wind turbine is regularly in derated operation. Afterwards the
turbine was down for maintenance and the gearbox was replaced.

Figure 7: Case Study 3: Anomalies cause by a failure in the gearbox of a wind
turbine of WF B. Using only 1 to 2 months of data for fine-tuning, the anomaly
detection model is capable of detecting the failure correctly. The upper graph
depicts the anomaly score, threshold and resulting criticality measure for the
baseline model, the lower graph compares the criticality for the three models.

##  5 Conclusions

To investigate the feasibility of autoencoder-based anomaly detection with
limited training data for wind turbines, two transfer learning applications
using tree methods were conducted. The study focused on anomaly detection
within a single wind farm, using different amounts of tuning data and wind
turbine types from two wind farms.

Two types of TL models were studied, asset-to-asset and multi-asset TL models.
In case of asset-to-asset the source models were trained on a single wind
turbine and fine-tuned on data of the target wind turbine. For the multi-asset
TL models, the source models were trained on data of multiple wind turbines
within the same wind farm and fine-tuned on data of the target wind turbine.

In both cases, the TL models fine-tuned using 1 to 3 months of data performed
as well or slightly better than the baseline models, which were trained on a
full year of data. The asset-to-asset TL models showed superior performance
compared to the multi-asset TL models. Additionally, three TL methods were
investigated: only updating the threshold, fine-tuning the decoder and fine-
tuning the autoencoder, of which the decoder-method showed the best results.

To demonstrate that early fault detection in wind turbines is indeed possible
even with a limited amount of training data three case studies were conducted.
In all three cases the failure was detected by TL models tuned on only 1 to 2
months of data, highlighting the potential of TL for improving anomaly
detection in wind turbines, reducing the amount of data and resources needed
for effective detection.

##### Acknowledgements

The development of methods presented was funded by the German Federal Ministry
for Economic Affairs and Climate Action (BMWK) through the research project
“ADWENTURE” (FKZ 03EE2030).

## References

  * [1] A. Stetco, F. Dinmohammadi, X. Zhao, V. Robu, D. Flynn, M. Barnes, J. Keane, G. Nenadic, Machine learning methods for wind turbine condition monitoring: A review, Renewable Energy 133 (2019) 620–635.  doi:10.1016/j.renene.2018.10.047.
URL https://www.sciencedirect.com/science/article/pii/S096014811831231X

  * [2] K. Weiss, T. M. Khoshgoftaar, D. Wang, A survey of transfer learning, Journal of Big Data 3 (1) (2016).  doi:10.1186/s40537-016-0043-6.
  * [3] F. Zhuang, Z. Qi, K. Duan, D. Xi, Y. Zhu, H. Zhu, H. Xiong, Q. He, A comprehensive survey on transfer learning, Proceedings of the IEEE 109 (1) (2021) 43–76.  doi:10.1109/JPROC.2020.3004555.
  * [4] G. Michau, O. Fink, Unsupervised transfer learning for anomaly detection: Application to complementary operating condition transfer, Knowledge-Based Systems 216 (2021) 106816.  doi:10.1016/j.knosys.2021.106816.
  * [5] J. Zgraggen, M. Ulmer, E. Jarlskog, G. Pizza, L. Goren Huber, Transfer learning approaches for wind turbine fault detection using deep learning, PHM Society European Conference 6 (1) (2021) 12.  doi:10.36001/phme.2021.v6i1.2835.
URL http://papers.phmsociety.org/index.php/phme/article/view/2835

  * [6] D. Chen, S. Yang, F. Zhou, Transfer learning based fault diagnosis with missing data due to multi-rate sampling, Sensors (Basel, Switzerland) 19 (8) (2019).  doi:10.3390/s19081826.
  * [7] L. Schröder, N. K. Dimitrov, D. R. Verelst, J. A. Sørensen, Using transfer learning to build physics-informed machine learning models for improved wind farm monitoring, Energies 15 (2) (2022).  doi:10.3390/en15020558.
URL https://www.mdpi.com/1996-1073/15/2/558

  * [8] T. Wen, R. Keyes, Time series anomaly detection using convolutional neural networks and transfer learning (2019).
URL http://arxiv.org/pdf/1905.13628v1

  * [9] S. Simani, S. Farsoni, P. Castaldi, Supervisory control and data acquisition for fault diagnosis of wind turbines via deep transfer learning, Energies 16 (9) (2023).  doi:10.3390/en16093644.
URL https://www.mdpi.com/1996-1073/16/9/3644

  * [10] C. M. Roelofs, M.-A. Lutz, S. Faulstich, S. Vogt, Autoencoder-based anomaly root cause analysis for wind turbines, Energy and AI 4 (2021) 100065.  doi:10.1016/j.egyai.2021.100065.
  * [11] Lutz, A., et al., Evaluation of anomaly detection of an autoencoder based on maintenace information and scada-data, Energies (2020).
URL https://www.mdpi.com/1996-1073/13/5/1063/htm

  * [12] K. He, X. Zhang, S. Ren, J. Sun, Delving deep into rectifiers: Surpassing human-level performance on imagenet classification, in: Proceedings of the IEEE international conference on computer vision, 2015, pp. 1026–1034.
  * [13] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization (2014).
URL http://arxiv.org/pdf/1412.6980v9

  *[NBM]: normal behaviour model
  *[SCADA]: Supervisory Control and Data Acquisition
  *[TL]: Transfer Learning
  *[OP-mode]: operational mode
  *[OP-modes]: operational mode
  *[AE]: autoencoder
  *[RMSE]: root mean square error


# Overlap Times in the $GI^{B}/GI/\infty$ Queue

Sergio Palomo
Systems Engineering
Cornell University
<EMAIL_ADDRESS>Jamol Pender 111Corresponding Author
School of Operations Research and Information Engineering
Cornell University
<EMAIL_ADDRESS>
###### Abstract

Overlap times have been studied as a way of understanding the time of
interaction between customers in a service facility. Most of the previous
analysis relies on the single jump assumption for arrivals, which implies the
queue increases by one for each arrival epoch. In this paper, we relax the
single arrival assumption and explore the impact of having batch arrivals.
Unfortunately, with batch arrivals it is not clear how one measures an overlap
time between batches of customers. Thus, we develop two ways of capturing the
notion of an overlap time in a batch setting and derive exact results in the
infinite server queue with batch arrivals. Finally, we derive new results for
analyzing overlap times of more than two batches.

##  1 Introduction

The study of overlap times has emerged as a way of understanding how customers
interact in service systems. Recently, overlap times have been analyzed in
infinite server queues, single server queues, and multi-server queues in [7,
13, 14, 9] respectively. This has been done mostly in Markovian cases since
the analysis is a bit easier. In [7] the overlap times were used as a way of
computing a new $R_{0}$ value for understanding infection rates in
compartmentalized epidemic models in relation to COVID-19 transmission. Palomo
and Pender [13] proves that the overlap distribution is exponential for the
M/M/1 queue and [14] studies the overlap distribution for infinite server
queues that have renewal arrivals and i.i.d service times from a renewal
process. Finally, Ko and Xu [9] uses fluid limits to approximate the
overlapping time of customers in multi-server queues with time varying arrival
rates.

Beyond infectious disease transmission and design of service systems there are
other applications where overlap times are important. One application is
public transportation like buses and subways. In these transportation
settings, several passengers get on the bus or subway in batches and ride for
some time to their destination. Thus, it is important to study the impact of
batch arrivals to get a better understanding of how infectious disease spreads
in these transportation settings.

Fortunately, the study of batch arrival queues has increased significantly in
recent years to due the emerging application of cloud-based data storage and
processing. In this setting, the batches of jobs arriving to the system are
collections of jobs submitted simultaneously by a user. These jobs are then
served by each being processed individually and the results are returned to
the user. For more discussion, see works such as Lu et al. [10], Pender and
Phung-Duc [19], Xie et al. [23], Yekkehkhany et al. [24] and references
therein. Another relevant application in the context of infectious disease
modeling is modeling the number infections of COVID-19. See for example Kaplan
[8], Morozova et al. [12], Palomo et al. [15]. In this setting, the results
for patients who potentially have COVID-19 arrive in a large batch to be
processed at a facility. Moreover, the data that we observe from COVID-19 is
also of batch form as counts are made daily. Finally, an emerging application
of batch queues is in context of autonomous vehicles moving in platoons
(batches) down highways and roads, e.g. Mirzaeian et al. [11], Hampshire et
al. [5]. Such applications also serve as the inspiration for the batch arrival
queue staffing problem studied by Daw et al. [3].

In this work, we study overlap times for the $GI^{B}/GI/\infty$ queue. Our
results on the $GI^{B}/GI/\infty$ queue also serve as approximations and lower
bounds for systems of finite capacity such as the $GI^{B}/GI/C$ queue. For
example, consider the recent work on multi-server jobs, which are queueing
systems where collections of jobs arrive together and also have a requirement
that they must start together (see, e.g., Rumyantsev and Morozov [20], Pender
and Ko [18], Afanaseva et al. [1], Grosof et al. [4], Weng and Wang [22], Hong
and Wang [6], Wang et al. [21], and references therein). The simultaneous
start requirement is a salient model feature relative to batch arrival many
server queues, but if there were infinitely many servers available then these
models reduce to one another. The multi-server setting model is known to be
quite challenging to analyze, so we offer our following analysis in the
unlimited server setting to understand the challenges and provide insights
about more complicated models.

In what follows we describe the contributions of our work and how the rest of
the paper is organized.

###  1.1 Contributions of the Paper

  * •

We derive explicit expressions for the steady state distribution for the
overlap time for individual customers that are exactly $k$ batches apart.

  * •

We derive exact expressions for the tail distribution of the overlap time when
we measure overlap in the first of each batch or last of each batch setting.

  * •

We also derive explicit expressions for the tail distribution of the overlap
time when there are more than two customers or two batches.

###  1.2 Organization of Paper

In Section 2, we describe the stochastic model that we will use in this work
and the definition of overlap times. We derive an exact expression for
describing the overlap times for a pair of individual customers in the
infinite server queue with batch arrivals i.e the $GI^{B}/GI/\infty$ queue. We
use this equation to compute the steady state distribution of the overlap time
of customers that are $k$ batches apart. In Section 3, we also analyze the
overlap times from a batch perspective. In particular we analyze two ways of
comparing a batch. The first is the first of each batch to leave and the
second is the last of each batch to leave. In Section 4, we analyze the
overlap times of when there are more than two customers. Finally in Section 5,
we provide a conclusion and some future research directions.

##  2 Infinite Server Overlap Times

In this section, we study the infinite server queue with the intention of
understanding how much time do adjacent customers spend in the system
together. A similar type of analysis has been completed Kang et al. [7],
Palomo and Pender [13, 14]. Here we consider the $GI^{B}$/GI/$\infty$ queue
where we let $A_{i}$ be the arrival time of the $i^{th}$ batch of customers,
the inter-arrival time between the $i^{th}$ and $(i+1)^{th}$ batches is given
by $A_{i+1}-A_{i}$. We assume the inter-arrival times are i.i.d with
cumulative distribution function (cdf) $H(x)$. For convenience, we define
$h_{k}(x)$ as the distribution of a sum of $k$ inter-arrival times. We also
assume that $S_{(i,j)}$ be the service time of the $i^{th}$ customer in the
$j^{th}$ batch and the service times are i.i.d with cdf $G(x)$. Thus, the
departure time for the $j^{th}$ customer in the $n^{th}$ batch is given by the
following expression

$\displaystyle D_{(j,n)}$ $\displaystyle=$ $\displaystyle S_{(j,n)}+A_{n}.$
(1)

Given the arrival and departure time for each customer, it is possible to
describe the overlap time between any pair of customers. The overlap time
between the $j^{th}$ customer in the $n^{th}$ batch and $\ell^{th}$ customer
in the $(n+k)^{th}$ batch is given by

$\displaystyle O_{(j,n),(\ell,n+k)}$ $\displaystyle=$
$\displaystyle\left(\min(D_{(j,n)},D_{(\ell,n+k)})-A_{n+k}\right)^{+}$ (2)
$\displaystyle=$
$\displaystyle\left(\min(A_{n}+S_{(j,n)},A_{n+k}+S_{(\ell,n+k)})-A_{n+k}\right)^{+}$
(3) $\displaystyle=$ $\displaystyle\left((D_{(j,n)}-A_{n+k})^{+}\wedge
S_{(\ell,n+k)}\right)$ (4) $\displaystyle=$
$\displaystyle\left((S_{(j,n)}+A_{n}-A_{n+k})^{+}\wedge S_{(\ell,n+k)}\right)$
(5) $\displaystyle=$
$\displaystyle\left((S_{(j,n)}-\left(A_{n+k}-A_{n}\right))^{+}\wedge
S_{(\ell,n+k)}\right).$ (6)

######  Theorem 2.1.

Let $O_{j,k,\ell}$ have the steady state distribution of
$O_{(j,n),(\ell,n+k)}$ in the $GI^{B}/GI/\infty$ queue, $\mathcal{S}$ and
$\tilde{\mathcal{S}}$ be two independent service times with cdf $G(x)$, and
let $\mathbb{P}(B=j)=p_{j}$, then the tail distribution of
$O_{j,k,\ell}=\lim_{n\to\infty}O_{(j,n),(\ell,n+k)}$ is given by

$\mathbb{P}\left(O_{j,k,\ell}>t\right)=\left\\{\begin{array}[]{lr}\overline{G}(t),&\text{for
}k=0,j=\ell,\ell\geq 1\\\ \overline{G}(t)^{2},&\text{for
}k=0,j\neq\ell,\ell\geq 1\\\
\overline{G}(t)\int^{\infty}_{0}\overline{G}\left(t+x\right)h_{k}(x)dx,&\text{for
}k>0,j\neq\ell,\ell\geq 1\end{array}\right\\}$

where $h_{k}(x)$ is the density of the sum of $k$ i.i.d inter-arrival times.

###### Proof.

In the first condition where we have $k=0,j=\ell,\ell\geq 1$, this corresponds
to a self-overlap, which is precisely equal to the tail distribution of the
service time. The second condition where we have $k=0,j\neq\ell,\ell\geq 1$,
this corresponds to an overlap of two customers in the same batch. Thus, the
probability that they are both around by time $t$ is equal to the square of
the tail distribution of the service time. Finally, the last condition is
exactly similar to the case with no batches, which is given in Palomo and
Pender [14]. This completes the proof. ∎

Theorem 2.1 provides the distribution of the case where we are comparing the
overlap between individual customers. However, we might be interested in the
overlap of a batch itself. In the sequel, we provide an analysis of this batch
perspective.

##  3 The Batch Perspective

Unfortunately in the batch setting, there is not a unique way to define an
overlap time between batches of customers. In what follows, we describe two
ways to define an overlap between customers in a batch. Each one yields
different expressions for calculating the tail distribution of the respective
overlap time definition. We recognize that there are many more ways to define
an interaction or overlap between a batch of customers, however, these two are
quite natural. However, before we derive the overlap times for the batch
queues, we prove a simple result about the tail distribution of the maximum
and minimum of a random sequence of random variables.

######  Lemma 3.1.

Let $B$ be a discrete random variable on the positive integers i.e.
$\\{1,2,3...\\}$ with probabilities {$p_{1},p_{2},...p_{n}$} and let
$X=(X_{1},X_{2},...,X_{B})$ be a random vector with random size given by the
random variable $B$. Then, we have the following distributions for $\min(X)$
and $\max(X)$

$\displaystyle\mathbb{P}\left(\min(X)>t\right)$ $\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}p_{j}\overline{G}(t)^{j}$ (7)

and

$\displaystyle\mathbb{P}\left(\max(X)>t\right)$ $\displaystyle=$
$\displaystyle 1-\sum^{\infty}_{j=1}p_{j}G(t)^{j}.$ (8)

###### Proof.

For the minimum we have that

$\displaystyle\mathbb{P}\left(\min(X)>t\right)$ $\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}\mathbb{P}\left(\min(X)>t|B=j\right)\cdot\mathbb{P}\left(B=j\right)$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}\prod^{j}_{i=1}\mathbb{P}\left(X_{i}>t\right)\cdot
p_{j}$ $\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}p_{j}\overline{G}(t)^{j}.$

Finally, for the maximum we have that

$\displaystyle\mathbb{P}\left(\max(X)>t\right)$ $\displaystyle=$
$\displaystyle 1-\mathbb{P}\left(\max(X)\leq t\right)$ $\displaystyle=$
$\displaystyle 1-\sum^{\infty}_{j=1}\mathbb{P}\left(\max(X)\leq
t|B=j\right)\cdot\mathbb{P}\left(B=j\right)$ $\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}\prod^{j}_{i=1}\mathbb{P}\left(X_{i}\leq
t\right)\cdot p_{j}$ $\displaystyle=$ $\displaystyle
1-\sum^{\infty}_{j=1}p_{j}G(t)^{j}.$

∎

The reader should note two things. First, when the batch size is one with
probability one, then this reduces to the directly to the tail distribution of
the random variable and this holds for either the max or min since they are
equivalent for one random variable. Second, we note that when the random
variable $B$ has finite support, then the infinite sums reduce to finite ones.
The finite support of the random variable $B$ is important from a
computational perspective since it often can be the difference between an
exact expression and an approximation of it.

###  3.1 Last of Each Batch to Leave

The first method for capturing overlap between customers in batches is called
”last of each batch to leave”. This method defines an overlap for a batch as
any time that any of customers of each batch overlap. Mathematically, if we
define $\vec{D}_{n}^{B_{n}}=(D_{n,1},D_{n,2},...,D_{n,B_{n}})$ as the vector
of departure times of the $n^{th}$ batch and
$\vec{S}_{n}^{B_{n}}=(S_{n,1},S_{n,2},...,S_{n,B_{n}})$ as the vector of
service times of the $n^{th}$ batch, then we have the following expression of
the overlap time of the ”last of each batch to leave”

$\displaystyle O_{n,n+k}$ $\displaystyle=$
$\displaystyle\left(\min\left(\max\left(\vec{D}_{n}^{B_{n}}\right),\max\left(\vec{D}_{n+k}^{B_{n+k}}\right)\right)-A_{n+k}\right)^{+}$
(9) $\displaystyle=$
$\displaystyle\left(\left(\max\left(\vec{S}_{n}^{B_{n}}\right)-\left(A_{n+k}-A_{n}\right)\right)^{+}\wedge\max\left(\vec{S}_{n+k}^{B_{n+k}}\right)\right).$
(10)

######  Theorem 3.2.

Let $O_{k}$ be the distribution of $O_{n,n+k}$ for all values of $n$ in the
$GI^{B}/GI/\infty$ queue and let $\mathcal{A}_{k}$, be the sum of k i.i.d
inter-arrival times from a renewal process with density $h_{k}(x)$, then the
tail distribution of $O_{k}$ under the last of each batch to leave overlap
process is given by

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\left(1-\mathbb{E}\left[G\left(t+\mathcal{A}_{k}\right)^{\mathcal{B}}\right]\right)\cdot\left(1-\mathbb{E}\left[G(t)^{\mathcal{B}}\right]\right).$
(11)

###### Proof.

First, we need to decompose the overlap probability into two probabilities by
using a property of the minimum of two independent random variables i.e.

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\mathbb{P}\left(\left(\left(\max\left(\vec{S}_{n}^{B_{n}}\right)-\left(A_{n+k}-A_{n}\right)\right)^{+}\wedge\max\left(\vec{S}_{n+k}^{B_{n+k}}\right)\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\left(\max\left(\vec{S}_{n}^{B_{n}}\right)-\left(A_{n+k}-A_{n}\right)\right)^{+}>t\right)\cdot\mathbb{P}\left(\max\left(\vec{S}_{n+k}^{B_{n+k}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\max\left(\vec{S}_{n}^{B_{n}}\right)>t+A_{n+k}-A_{n}\right)\cdot\mathbb{P}\left(\max\left(\vec{S}_{n+k}^{B_{n+k}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[1-\sum^{\infty}_{j=1}p_{j}G\left(t+A_{n+k}-A_{n}\right)^{j}\right]\cdot\left(1-\sum^{\infty}_{j=1}p_{j}G\left(t\right)^{j}\right)$
$\displaystyle=$
$\displaystyle\left(1-\sum^{\infty}_{j=1}p_{j}\mathbb{E}\left[G\left(t+A_{n+k}-A_{n}\right)^{j}\right]\right)\cdot\left(1-\sum^{\infty}_{j=1}p_{j}G\left(t\right)^{j}\right)$
$\displaystyle=$
$\displaystyle\left(1-\sum^{\infty}_{j=1}p_{j}\left(\int^{\infty}_{0}G\left(t+x\right)^{j}h_{k}(x)dx\right)\right)\cdot\left(1-\sum^{\infty}_{j=1}p_{j}G\left(t\right)^{j}\right)$
$\displaystyle=$
$\displaystyle\left(1-\mathbb{E}\left[G\left(t+\mathcal{A}_{k}\right)^{\mathcal{B}}\right]\right)\cdot\left(1-\mathbb{E}\left[G(t)^{\mathcal{B}}\right]\right).$

This completes the proof. ∎

######  Corollary 3.3.

Let $O_{k}$ be the distribution of $O_{n,n+k}$ for all values of $n$ in the
$M^{\mathcal{B}}/M/\infty$ queue, then the tail distribution of $O_{k}$ under
the last of each batch to leave setting is given by

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\left(1-\mathbb{E}\left[\sum^{\mathcal{B}}_{m=0}{\mathcal{B}\choose
m}(-1)^{m}e^{-\mu tm}\left(\frac{\lambda}{\lambda+\mu
m}\right)^{k}\right]\right)\cdot\left(1-\mathbb{E}\left[(1-e^{-\mu
t})^{\mathcal{B}}\right]\right).$ (12)

###### Proof.

Using the previous result in Theorem 3.2, we have that

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\left(1-\sum^{\infty}_{j=1}p_{j}\int^{\infty}_{0}\left(1-e^{-\mu(t+x)}\right)^{j}\frac{(\lambda
x)^{k-1}}{(k-1)!}\lambda e^{-\lambda
x}dx\right)\cdot\left(1-\sum^{\infty}_{j=1}p_{j}(1-e^{-\mu t})^{j}\right)$
$\displaystyle=$
$\displaystyle\left(1-\sum^{\infty}_{j=1}p_{j}\int^{\infty}_{0}\left(\sum^{j}_{m=0}{j\choose
m}(-1)^{m}e^{-\mu(t+x)m}\right)\frac{(\lambda x)^{k-1}}{(k-1)!}\lambda
e^{-\lambda x}dx\right)\cdot\left(1-\sum^{\infty}_{j=1}p_{j}(1-e^{-\mu
t})^{j}\right)$ $\displaystyle=$
$\displaystyle\left(1-\sum^{\infty}_{j=1}p_{j}\sum^{j}_{m=0}{j\choose
m}(-1)^{m}\int^{\infty}_{0}e^{-\mu(t+x)m}\frac{(\lambda
x)^{k-1}}{(k-1)!}\lambda e^{-\lambda
x}dx\right)\cdot\left(1-\sum^{\infty}_{j=1}p_{j}(1-e^{-\mu t})^{j}\right)$
$\displaystyle=$
$\displaystyle\left(1-\sum^{\infty}_{j=1}p_{j}\sum^{j}_{m=0}{j\choose
m}(-1)^{m}e^{-\mu tm}\left(\frac{\lambda}{\lambda+\mu
m}\right)^{k}\right)\cdot\left(1-\sum^{\infty}_{j=1}p_{j}(1-e^{-\mu
t})^{j}\right)$ $\displaystyle=$
$\displaystyle\left(1-\mathbb{E}\left[\sum^{\mathcal{B}}_{m=0}{\mathcal{B}\choose
m}(-1)^{m}e^{-\mu tm}\left(\frac{\lambda}{\lambda+\mu
m}\right)^{k}\right]\right)\cdot\left(1-\mathbb{E}\left[(1-e^{-\mu
t})^{\mathcal{B}}\right]\right).$

This completes the proof. ∎

######  Corollary 3.4.

Let $O_{k}$ be the distribution of $O_{n,n+k}$ for all values of $n$ in the
$GI^{B}/D/\infty$ queue, then the tail distribution of $O_{k}$ under the last
of each batch to leave setting is given by

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\mathbb{F}^{(k)}\left((\Delta-t)^{+}\right)$ (13)

where $\mathbb{F}^{(k)}(\cdot)$ is the $k$-fold convolution of the inter-
arrival distribution $F(x)$.

###### Proof.

Using the previous result in Theorem 3.2, we have that

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\mathbb{P}\left(\left(\left(\max\left(\vec{S}_{n}^{B_{n}}\right)-\left(A_{n+k}-A_{n}\right)\right)^{+}\wedge\max\left(\vec{S}_{n+k}^{B_{n+k}}\right)\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\left(\Delta-\left(A_{n+k}-A_{n}\right)\right)^{+}>t\right)\cdot\mathbb{P}\left(\Delta>t\right)$
$\displaystyle=$ $\displaystyle\mathbb{F}^{(k)}\left((\Delta-t)^{+}\right).$

This completes the proof. ∎

###  3.2 First of Each Batch to Leave

The second method for capturing overlap between customers in batches is called
”first of each batch to leave”. This method defines an overlap for a batch as
any time that all of the customers of each batch overlap. Mathematically, we
have the following expression of the overlap time of the ”first of each batch
to leave”

$\displaystyle O_{n,n+k}$ $\displaystyle=$
$\displaystyle\left(\min\left(\min\left(\vec{D}_{n}^{B_{n}}\right),\min\left(\vec{D}_{n+k}^{B_{n+k}}\right)\right)-A_{n+k}\right)^{+}$
(14) $\displaystyle=$
$\displaystyle\left(\left(\min\left(\vec{S}_{n}^{B_{n}}\right)-\left(A_{n+k}-A_{n}\right)\right)^{+}\wedge\min\left(\vec{S}_{n+k}^{B_{n+k}}\right)\right).$
(15)

######  Theorem 3.5.

Let $O_{k}$ be the distribution of $O_{n,n+k}$ for all values of $n$ in the
$GI^{B}/GI/\infty$ queue and $\mathcal{A}_{k}$, be the sum of k i.i.d inter-
arrival times from a renewal process with density $h_{k}(x)$, then the tail
distribution of $O_{k}$ under the first of each batch to leave setting is
given by

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\overline{G}(t+\mathcal{A}_{k})^{\mathcal{B}}\right]\cdot\mathbb{E}\left[\overline{G}(t)^{\mathcal{B}}\right].$
(16)

###### Proof.

First, we need to decompose the overlap probability into two probabilities by
using a property of the minimum of two independent random variables i.e.

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\mathbb{P}\left(\left(\left(\min\left(\vec{S}_{n}^{B_{n}}\right)-\left(A_{n+k}-A_{n}\right)\right)^{+}\wedge\min\left(\vec{S}_{n+k}^{B_{n+k}}\right)\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\left(\min\left(\vec{S}_{n}^{B_{n}}\right)-\left(A_{n+k}-A_{n}\right)\right)^{+}>t\right)\cdot\mathbb{P}\left(\min\left(\vec{S}_{n+k}^{B_{n+k}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\min\left(\vec{S}_{n}^{B_{n}}\right)>t+A_{n+k}-A_{n}\right)\cdot\mathbb{P}\left(\min\left(\vec{S}_{n+k}^{B_{n+k}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum^{\infty}_{j=1}p_{j}\overline{G}\left(t+A_{n+k}-A_{n}\right)^{j}\right]\cdot\left(\sum^{\infty}_{j=1}p_{j}\overline{G}\left(t\right)^{j}\right)$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}p_{j}\mathbb{E}\left[\overline{G}\left(t+A_{n+k}-A_{n}\right)^{j}\right]\cdot\left(\sum^{\infty}_{j=1}p_{j}\overline{G}\left(t\right)^{j}\right)$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}p_{j}\mathbb{E}\left[\overline{G}\left(t+\mathcal{A}_{k}\right)^{j}\right]\cdot\left(\sum^{\infty}_{j=1}p_{j}\overline{G}\left(t\right)^{j}\right)$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}p_{j}\left(\int^{\infty}_{0}\overline{G}\left(t+x\right)^{j}h_{k}(x)dx\right)\cdot\left(\sum^{\infty}_{j=1}p_{j}\overline{G}\left(t\right)^{j}\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\overline{G}(t+\mathcal{A}_{k})^{\mathcal{B}}\right]\cdot\mathbb{E}\left[\overline{G}(t)^{\mathcal{B}}\right].$

This completes the proof. ∎

######  Corollary 3.6.

Let $O_{k}$ be the distribution of $O_{n,n+k}$ for all values of $n$ in the
$M^{B}/M/\infty$ queue and $\mathcal{A}_{k}$, be the sum of k i.i.d inter-
arrival times from a renewal process with density $h_{k}(x)$, then the tail
distribution of $O_{k}$ under the first of each batch to leave setting is
given by

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[e^{-\mu
t\mathcal{B}}\left(\frac{\lambda}{\lambda+\mu\mathcal{B}}\right)^{k}\right]\cdot\mathbb{E}\left[e^{-\mu
t\mathcal{B}}\right].$ (17)

###### Proof.

Using the previous result in Theorem 3.5, we have that

$\displaystyle\mathbb{P}\left(O_{k}>t\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum^{\infty}_{j=1}p_{j}e^{-\mu
j(t+A_{n+k}-A_{n})}\right]\cdot\left(\sum^{\infty}_{j=1}p_{j}e^{-\mu
jt}\right)$ $\displaystyle=$ $\displaystyle\sum^{\infty}_{j=1}p_{j}e^{-\mu
jt}\mathbb{E}\left[e^{-\mu
j(A_{n+k}-A_{n})}\right]\cdot\left(\sum^{\infty}_{j=1}p_{j}e^{-\mu jt}\right)$
$\displaystyle=$ $\displaystyle\sum^{\infty}_{j=1}p_{j}e^{-\mu
jt}\mathbb{E}\left[e^{-\mu
j\mathcal{A}}\right]^{k}\cdot\left(\sum^{\infty}_{j=1}p_{j}e^{-\mu jt}\right)$
$\displaystyle=$ $\displaystyle\mathbb{E}\left[e^{-\mu
t\mathcal{B}}\left(\frac{\lambda}{\lambda+\mu\mathcal{B}}\right)^{k}\right]\cdot\mathbb{E}\left[e^{-\mu
t\mathcal{B}}\right].$

This completes the proof. ∎

##  4 More Than A Pair Setting

Our above analysis for overlap times only characterizes the overlap for two
customers or two batches. In some settings it may be interesting to
characterize the behavior of more than a pair of customers or a pair of
batches. We analyze the setting of more than two customers in what follows
below.

###  4.1 Individual Customer Perspective

For our first result, we characterize the overlap time of an arbitrary number
of customers in the $GI/GI/\infty$ queue. In particular, the overlap time for
a k-tuple of customers in the $GI/GI/\infty$ queue, which is denoted by
$O_{n_{1},n_{2},\dots,n_{k}}$, can be computed by the following expressions

$\displaystyle O_{n_{1},n_{2},\dots,n_{k}}$ $\displaystyle=$
$\displaystyle\left(\min(D_{n_{1}},D_{n_{2}},\cdots,D_{n_{k}})-A_{n_{k}}\right)^{+}$
$\displaystyle=$
$\displaystyle\left(\min(A_{n_{1}}+S_{n_{1}},\dots,A_{n_{j}}+S_{n_{j}},A_{n_{k}}+S_{n_{k}})-A_{n_{k}}\right)^{+}$
$\displaystyle=$
$\displaystyle(S_{n_{1}}+A_{n_{1}}-A_{n_{k}})^{+}\bm{\wedge}\cdots\bm{\wedge}(S_{n_{j}}+A_{n_{j}}-A_{n_{k}})^{+}\bm{\wedge}S_{n_{k}}.$

First, note that the distance between any two pairs of customers is arbitrary.
Moreover, unlike the two customer case, we observe that the random variables
are not independent anymore. This is because the inter-arrival times are
connected in each of the customers except for the last one, which only has a
service time. Nonetheless, in the exponential setting, we can still show that
the arrival process only serves to determine the constants and not the
functional form of the overlap distribution.

######  Theorem 4.1.

Let $O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $GI/GI/\infty$ queue, then the
tail distribution of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the following
formula

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{j=1}\overline{G}\left(t+(A_{n_{k}}-A_{n_{j}})\right)\right].$
(18)

###### Proof.

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left((S_{n_{1}}+A_{n_{1}}-A_{n_{k}})^{+}\bm{\wedge}\cdots\bm{\wedge}(S_{n_{j}}+A_{n_{j}}-A_{n_{k}})^{+}\bm{\wedge}S_{n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(S_{n_{1}}+A_{n_{1}}-A_{n_{k}}>t,\dots,S_{n_{j}}+A_{n_{j}}-A_{n_{k}}>t,S_{n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{j=1}\overline{G}\left(t+(A_{n_{k}}-A_{n_{j}})\right)\right].$

This completes the proof. ∎

######  Corollary 4.2.

Let $O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $GI/M/\infty$ queue, then the
tail distribution of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the following
formula

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$ $\displaystyle e^{-\mu
kt}\mathbb{E}\left[\prod^{k-1}_{j=1}e^{-\mu j(A_{n_{j+1}}-A_{n_{j}})}\right].$
(19)

###### Proof.

Using the result from Theorem 4.1, we have

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(S_{n_{1}}+A_{n_{1}}-A_{n_{k}}>t,\dots,S_{n_{j}}+A_{n_{j}}-A_{n_{k}}>t,S_{n_{k}}>t\right)$
$\displaystyle=$ $\displaystyle e^{-\mu
kt}\mathbb{E}\left[\prod^{k-1}_{j=1}e^{-\mu(A_{n_{k}}-A_{n_{j}})}\right]$
$\displaystyle=$ $\displaystyle e^{-\mu
kt}\mathbb{E}\left[\prod^{k-1}_{j=1}e^{-\mu j(A_{n_{j+1}}-A_{n_{j}})}\right].$

This completes the proof. ∎

######  Corollary 4.3.

Let $O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $M/M/\infty$ queue, then the
tail distribution of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the following
formula

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$ $\displaystyle e^{-\mu
kt}\prod^{k-1}_{j=1}\left(\frac{\lambda}{\lambda+\mu
j}\right)^{n_{j+1}-n_{j}}.$ (20)

###### Proof.

Using the result from Theorem 4.1, we have

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$ $\displaystyle e^{-\mu
kt}\mathbb{E}\left[\prod^{k-1}_{j=1}e^{-\mu j(A_{n_{j+1}}-A_{n_{j}})}\right]$
$\displaystyle=$ $\displaystyle e^{-\mu
kt}\prod^{k-1}_{j=1}\left(\frac{\lambda}{\lambda+\mu
j}\right)^{n_{j+1}-n_{j}}.$

This completes the proof. ∎

The first thing to notice is that our results match the two customer case and
yield the previous results by Palomo and Pender [14]. We also observe that the
functional form of the result is not affected by the arrival distribution and
is completely governed by the service distribution, which is exponential in
this case. Moreover, the more customers you want to overlap, the smaller the
probability since the rate in the exponential is increased according to the
number of overlapping customers.

###  4.2 First of Each Batch

Now that we have analyzed the k-tuple setting in the individual setting i.e.
the batch size is equal to 1, we now move to the setting where batches can
have an arbitrary size. We will also perform our analysis for the two batch
settings we define in the previous sections i.e. first of each batch and last
of each batch.

######  Theorem 4.4.

Under the first of each batch to leave setting, we let
$O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $GI^{B}/GI/\infty$ queue, then
the tail distribution of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the
following formula

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{i=1}\overline{G}\left(t-(A_{n_{k}}-A_{n_{i}})\right)^{\mathcal{B}}\right].$
(21)

###### Proof.

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left((\min\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)+A_{n_{1}}-A_{n_{k}})^{+}\bm{\wedge}\cdots\bm{\wedge}(\min\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)+A_{n_{j}}-A_{n_{k}})^{+}\bm{\wedge}\min\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\min\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)+A_{n_{1}}-A_{n_{k}}>t,\dots,\min\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)+A_{n_{j}}-A_{n_{k}}>t,\min\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\min\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)>t+A_{n_{k}}-A_{n_{1}},\dots,\min\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)>t+A_{n_{k}}-A_{n_{j}},\min\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{i=1}\left(\sum^{\infty}_{j=1}p_{j}\cdot\overline{G}\left(t-(A_{n_{k}}-A_{n_{i}})\right)^{j}\right)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{i=1}\overline{G}\left(t-(A_{n_{k}}-A_{n_{i}})\right)^{\mathcal{B}}\right].$

This completes the proof. ∎

######  Corollary 4.5.

Under the first of each batch to leave setting, we let
$O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $GI^{B}/M/\infty$ queue, then
the tail distribution of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the
following formula

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$ $\displaystyle\mathbb{E}\left[e^{-\mu
kt\mathcal{B}}\prod^{k-1}_{i=1}e^{-\mu\mathcal{B}(A_{n_{k}}-A_{n_{i}})}\right].$
(22)

###### Proof.

Using the result from Theorem 4.4, we have that

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum^{\infty}_{j=1}p_{j}\left(\prod^{k-1}_{i=1}e^{-\mu
j(t+A_{n_{k}}-A_{n_{i}})}\cdot e^{-\mu jt}\right)\right]$ $\displaystyle=$
$\displaystyle\sum^{\infty}_{j=1}p_{j}\cdot e^{-\mu
jkt}\mathbb{E}\left[\prod^{k-1}_{i=1}e^{-\mu j(A_{n_{k}}-A_{n_{i}})}\right]$
$\displaystyle=$ $\displaystyle\mathbb{E}\left[e^{-\mu
kt\mathcal{B}}\prod^{k-1}_{i=1}e^{-\mu\mathcal{B}(A_{n_{k}}-A_{n_{i}})}\right].$

This completes the proof. ∎

######  Corollary 4.6.

Under the first of each batch to leave setting, we let
$O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $GI^{B}/M/\infty$ queue, then
the mean of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the following formula

$\displaystyle\mathbb{E}\left[O_{n_{1},n_{2},\dots,n_{k}}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{\mu
k}\mathbb{E}\left[\mathcal{B}^{-1}\prod^{k-1}_{i=1}e^{-\mu\mathcal{B}(A_{n_{k}}-A_{n_{i}})}\right].$
(23)

###### Proof.

Using the result from Theorem 4.4, we have that

$\displaystyle\mathbb{E}\left[O_{n_{1},n_{2},\dots,n_{k}}\right]$
$\displaystyle=$
$\displaystyle\int^{\infty}_{0}\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)dt$
$\displaystyle=$ $\displaystyle\int^{\infty}_{0}\mathbb{E}\left[e^{-\mu
kt\mathcal{B}}\prod^{k-1}_{i=1}e^{-\mu\mathcal{B}(A_{n_{k}}-A_{n_{i}})}\right]dt$
$\displaystyle=$ $\displaystyle\frac{1}{\mu
k}\mathbb{E}\left[\mathcal{B}^{-1}\prod^{k-1}_{i=1}e^{-\mu\mathcal{B}(A_{n_{k}}-A_{n_{i}})}\right].$

This completes the proof. ∎

###  4.3 Last of Each Batch

######  Theorem 4.7.

Under the last of each batch to leave setting, we let
$O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $GI^{B}/GI/\infty$ queue, then
the tail distribution of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the
following formula

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{i=1}\left(1-G\left(t-(A_{n_{k}}-A_{n_{i}})\right)^{\mathcal{B}}\right)\right].$
(24)

###### Proof.

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left((\max\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)+A_{n_{1}}-A_{n_{k}})^{+}\bm{\wedge}\cdots\bm{\wedge}(\max\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)+A_{n_{j}}-A_{n_{k}})^{+}\bm{\wedge}\max\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\max\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)+A_{n_{1}}-A_{n_{k}}>t,\dots,\max\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)+A_{n_{j}}-A_{n_{k}}>t,\max\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\max\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)>t+A_{n_{k}}-A_{n_{1}},\dots,\max\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)>t+A_{n_{k}}-A_{n_{j}},\max\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{i=1}\left(1-\sum^{\infty}_{j=1}p_{j}\left((1-\overline{G}(t-(A_{n_{k}}-A_{n_{i}}))\right)^{j}\right)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod^{k}_{i=1}\left(1-G\left(t-(A_{n_{k}}-A_{n_{i}})\right)^{\mathcal{B}}\right)\right].$

This completes the proof. ∎

######  Corollary 4.8.

Under the last of each batch to leave setting, we let
$O_{n_{1},n_{2},\dots,n_{k}}$ be the distribution of the overlap time for
customers within distance $n_{k}-n_{1}$ in the $GI^{B}/M/\infty$ queue, then
the tail distribution of $O_{n_{1},n_{2},\dots,n_{k}}$ is given by the
following formula

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$ $\displaystyle
1-\mathbb{E}\left[\left(\prod^{k}_{i=1}\left(1-e^{-\mu\mathcal{B}\left(t+A_{n_{k}}-A_{n_{i}}\right)}\right)\right)\right].$
(25)

###### Proof.

Using the result from Theorem 4.7, we have that

$\displaystyle\mathbb{P}\left(O_{n_{1},n_{2},\dots,n_{k}}>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left((\max\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)+A_{n_{1}}-A_{n_{k}})^{+}\bm{\wedge}\cdots\bm{\wedge}(\max\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)+A_{n_{j}}-A_{n_{k}})^{+}\bm{\wedge}\max\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\max\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)+A_{n_{1}}-A_{n_{k}}>t,\dots,\max\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)+A_{n_{j}}-A_{n_{k}}>t,\max\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\max\left(\vec{S}_{n_{1}}^{B_{n_{1}}}\right)>t+A_{n_{k}}-A_{n_{1}},\dots,\max\left(\vec{S}_{n_{j}}^{B_{n_{j}}}\right)>t+A_{n_{k}}-A_{n_{j}},\max\left(\vec{S}_{n_{k}}^{B_{n_{k}}}\right)>t\right)$
$\displaystyle=$ $\displaystyle
1-\mathbb{E}\left[\sum^{\infty}_{j=1}p_{j}\left(\prod^{k}_{i=1}\left(1-e^{-\mu
j\left(t+A_{n_{k}}-A_{n_{i}}\right)}\right)\right)\right]$ $\displaystyle=$
$\displaystyle
1-\mathbb{E}\left[\left(\prod^{k}_{i=1}\left(1-e^{-\mu\mathcal{B}\left(t+A_{n_{k}}-A_{n_{i}}\right)}\right)\right)\right].$

This completes the proof. ∎

##  5 Conclusion

In this paper, we consider the overlap times for customers in the
$GI^{B}/GI/\infty$ queue. We provide an analysis from an individual and a
batch perspective. Despite our analysis, there are many avenues for additional
research. First, as we mentioned before, there are several ways to analyze the
batch perspective. One such way would be to explicitly analyze the order
statistics of each batch. In particular one might be interested in when a
certain percentage of customers have departed from each batch.

Second, one can also extend our work by considering generalizations such as
dependent arrivals or service times like in Pang and Whitt [16, 17], Daw and
Pender [2]. In particular, dependent service times would impact the tail
distribution quite significantly since it no longer splits into a product of
individual distributions. It would also be interesting to consider the
situation where there are a finite number of servers or abandonments to assess
the impact on the tail distribution of the overlap time, see for example Ko
and Xu [9]. Finally, we are also interested in designing queueing systems to
minimize overlap times under some cost criterion. We hope to consider these
important extensions in future work.

## Acknowledgements

Jamol Pender would like to acknowledge the gracious support of the National
Science Foundation DMS Award # 2206286. Sergio Palomo would like to
acknowledge the gracious support of the Sloan Foundation for supporting his
graduate studies.

## References

  * Afanaseva et al. [2020] Larisa Afanaseva, Elena Bashtova, and Svetlana Grishunina.  Stability analysis of a multi-server model with simultaneous service and a regenerative input flow.  _Methodology and Computing in Applied Probability_ , 22(4):1439–1455, 2020.
  * Daw and Pender [2018] Andrew Daw and Jamol Pender.  Queues driven by hawkes processes.  _Stochastic Systems_ , 8(3):192–229, 2018.
  * Daw et al. [2021] Andrew Daw, Robert C Hampshire, and Jamol Pender.  How to staff when customers arrive in batches.  _arXiv preprint arXiv:1907.12650_ , 2021.
  * Grosof et al. [2020] Isaac Grosof, Mor Harchol-Balter, and Alan Scheller-Wolf.  Stability for two-class multiserver-job systems.  _arXiv preprint arXiv:2010.00631_ , 2020.
  * Hampshire et al. [2020] Robert C Hampshire, Shan Bao, Walter S Lasecki, Andrew Daw, and Jamol Pender.  Beyond safety drivers: Applying air traffic control principles to support the deployment of driverless vehicles.  _PLoS one_ , 15(5):e0232837, 2020.
  * Hong and Wang [2021] Yige Hong and Weina Wang.  Sharp waiting-time bounds for multiserver jobs.  _arXiv preprint arXiv:2109.05343_ , 2021.
  * Kang et al. [2021] Kang Kang, Sherwin Doroudi, Mohammad Delasay, and Alexander Wickeham.  A queueing-theoretic framework for evaluating transmission risks in service facilities during a pandemic.  _arXiv preprint arXiv:2103.13441_ , 2021.
  * Kaplan [2020] Edward H Kaplan.  OM forum—COVID-19 scratch models to support local decisions.  _Manufacturing & Service Operations Management_, 22(4):645–655, 2020.
  * Ko and Xu [2022] Young Myoung Ko and Jin Xu.  Overlapping time of a virtual customer in time-varying many-server queues.  _arXiv preprint arXiv:2211.03962_ , 2022.
  * Lu et al. [2011] Yi Lu, Qiaomin Xie, Gabriel Kliot, Alan Geller, James R Larus, and Albert Greenberg.  Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services.  _Performance Evaluation_ , 68(11):1056–1071, 2011\.
  * Mirzaeian et al. [2021] Neda Mirzaeian, Soo-Haeng Cho, and Alan Scheller-Wolf.  A queueing model and analysis for autonomous vehicles on highways.  _Management Science_ , 67(5):2904–2923, 2021\.
  * Morozova et al. [2020] Olga Morozova, Zehang Richard Li, and Forrest W Crawford.  A model for COVID-19 transmission in connecticut.  _medRxiv_ , 2020.
  * Palomo and Pender [2021a] Sergio Palomo and Jamol Pender.  Measuring the overlap with other customers in the G/G/1 queue.  In KH Bae, B Feng, S Kim, S Lazarova-Molnar, Z Zheng, T Roeder, and R Thiesing, editors, _Submitted to the Proceedings of the 2021 Winter Simulation Conference_ , pages 2595–2605. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc., 2021a.
  * Palomo and Pender [2021b] Sergio Palomo and Jamol Pender.  Overlap times in the infinite server queue.  _Probability in the Engineering and Informational Sciences To Appear_ , 2021b.
  * Palomo et al. [2020] Sergio Palomo, Jamol Pender, William Massey, and Robert C Hampshire.  Flattening the curve: Insights from queueing theory.  _arXiv preprint arXiv:2004.09645_ , 2020.
  * Pang and Whitt [2012] Guodong Pang and Ward Whitt.  The impact of dependent service times on large-scale service systems.  _Manufacturing & Service Operations Management_, 14(2):262–278, 2012.
  * Pang and Whitt [2013] Guodong Pang and Ward Whitt.  Two-parameter heavy-traffic limits for infinite-server queues with dependent service times.  _Queueing Systems_ , 73(2):119–146, 2013.
  * Pender and Ko [2017] Jamol Pender and Young Myoung Ko.  Approximations for the queue length distributions of time-varying many-server queues.  _INFORMS Journal on Computing_ , 29(4):688–704, 2017.
  * Pender and Phung-Duc [2016] Jamol Pender and Tuan Phung-Duc.  A law of large numbers for m/m/c/delayoff-setup queues with nonstationary arrivals.  In _International conference on analytical and stochastic modeling techniques and applications_ , pages 253–268. Springer, 2016.
  * Rumyantsev and Morozov [2017] Alexander Rumyantsev and Evsey Morozov.  Stability criterion of a multiserver model with simultaneous service.  _Annals of Operations Research_ , 252(1):29–39, 2017.
  * Wang et al. [2021] Weina Wang, Qiaomin Xie, and Mor Harchol-Balter.  Zero queueing for multi-server jobs.  In _Abstract Proceedings of the 2021 ACM SIGMETRICS/International Conference on Measurement and Modeling of Computer Systems_ , pages 13–14, 2021.
  * Weng and Wang [2020] Wentao Weng and Weina Wang.  Achieving zero asymptotic queueing delay for parallel jobs.  _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ , 4(3):1–36, 2020.
  * Xie et al. [2017] Qiaomin Xie, Mayank Pundir, Yi Lu, Cristina L Abad, and Roy H Campbell.  Pandas: Robust locality-aware scheduling with stochastic delay optimality.  _IEEE/ACM Transactions on Networking (TON)_ , 25(2):662–675, 2017.
  * Yekkehkhany et al. [2018] Ali Yekkehkhany, Avesta Hojjati, and Mohammad H Hajiesmaili.  GB-PANDAS:: Throughput and heavy-traffic optimality analysis for affinity scheduling.  _ACM SIGMETRICS Performance Evaluation Review_ , 45(2):2–14, 2018.