Translating Empirical State-Dependent Service Times Into Queueing Models

Abstract

Recent empirical studies suggest that human behavior in queues causes workload-dependent service times. We investigate the translation of empirical service times into state-dependent queueing models. To this end, we identify two types of state-dependent models, static and dynamic, and two types of corresponding behavioral mechanisms. For example, we view customer early task initiation as a static mechanism and social speedup pressure as a dynamic mechanism. For each model type, we discuss behavioral mechanisms consistent with the model assumptions and indicate how empirical service times can be translated into model input parameters. We illustrate how translating service times into dynamic models can result in invalid service rates, which provides evidence against dynamic mechanisms. For dynamic models, we find that mean service times are in general not the inverse of service rates, the directional change in service rates is not always the opposite of the directional change in mean service times, and workload measurement timing can drastically impact mean service time patterns. We provide closed-form equations to convert service times into service rates and vice versa, and find conditions under which monotonic mean service times imply monotonic service rates and vice versa. Our results provide guidelines for researchers to select and specify an appropriate state-dependent queueing model from service time data, and expand the scope of previously published analytical results.

1. Introduction

Recent empirical work has shown that human customers and servers, for a broad spectrum of service operations, behave adaptively. The performance (speed and work content; see Sections 2 and 3.4) of servers who respond to changes in the number of customers in the system (workload) is typically studied through estimated workload-dependent mean service times. It is not always obvious, however, how to correctly interpret patterns in mean service times. We show that, in some situations, it is incorrect to estimate server performance directly from mean service times. For example, mean service times moving in one direction do not necessarily imply that server speed moves in the opposite direction, holding the work content constant.

Server behavioral response to workload can be static or dynamic. By static behavior, we mean that server performance remains constant during each service. By dynamic behavior, we mean that server performance changes during a service, as the workload changes. We provide examples of static and dynamic mechanisms that have been investigated empirically in Section 3.4. As a preview, for example, we will argue that customer early task initiation is a static mechanism because it only impacts server performance before the beginning of service; in contrast, social speedup pressure is a dynamic mechanism because the server may speed up or slow down during the service.

Figure 1.

Numerical examples (see EC.2 of E-Companion for details): (a) monotonicity is preserved and (b) monotonicity is not preserved.

Static and dynamic behaviors are incorporated in static and dynamic queueing models. We limit our scope to models that quantify workload as the number of customers in the system, $n$ , which we refer to as the workload or the state. In a static model, a service duration depends on the state at a single point in time, typically at the beginning of service, and is independent of state changes during the service. A static model is specified through a set of state-dependent service time distributions. In a dynamic model, a service duration can depend on state changes during the service. A dynamic model is specified through a set of state-dependent service rates. We limit our attention to static and dynamic models that make certain Markovian assumptions, which are stated precisely in Section 3.1. For a static model, mean service times equal the inverse of server speed, holding work content constant. This is not true for a dynamic model, as we demonstrate in Section 4.

Ideally, empirical work informs and inspires analytical researchers to develop more useful mathematical models, which in turn inspires further empirical research, in a positive feedback loop (Fisher, 2007; Fisher et al., 2020). This loop requires empiricists to identify an appropriate model and translate findings correctly, and analysts to formulate and validate models using empirical data. We focus on connecting these two streams of work. Our work aims to provide methods for the proper translation of empirical results into analytical models. In particular, we investigate the choice of an appropriate state-dependent model and how to translate mean service times from empirical data into the inputs needed to specify the model.

Although translating mean service times into static models is relatively straightforward, this is not the case for dynamic models. Our results on dynamic models contribute to the empirical literature in three ways. First, we identify a potential pitfall in estimating service rates from mean service times. In classical queueing models, such as $M/M/c$ , the service rate is the inverse of the mean service time: $\mu =1/\text{E}(S)$ . It is tempting, but incorrect, to generalize this formula to:

$${\mu }_n=\frac{1}{\mathrm{E}(S_n)},n\ge 1,\leftarrow \text{Faulty generalization 1},$$

(1)

for state $n$ of a dynamic $M/M_n/c$ model (see Section 3.1 for definitions of ${\mu }_n$ and $S_n$ ). We refer to (1) as a faulty generalization. The intuitive reason that (1) does not hold is that during a customer’s service time, the system state can change, and therefore the service rate can change. We derive the correct relationship between ${\mu }_n$ and $\text{E}(S_n)$ (see (6)) and illustrate the magnitude of the errors that result from using (1) as an approximation for ${\mu }_n$ (see Section 6.3). Figure 1(a) shows a first example, in which the average approximation error is $24.0\mathrm{\%}$ . We expect this result to help empirical researchers better estimate the magnitude of the effect of workload on service rate.

Second, we highlight the possibility that mean service time data analyses might not reveal the correct directional change in service rates. Suppose that the variation in service rate is of interest, but the empirical analysis focuses on whether the amount by which mean service times change, if $n$ increases by one unit, is positive, negative, or zero. In this setting, it is tempting, but incorrect, to assume that service rates and mean service times move in opposite directions:

$$\begin{array}{rl}\mathrm{sgn}({\mu }_{n+1}-{\mu }_n)& =-\mathrm{sgn}(\text{E}(S_{n+1})-\text{E}(S_n)),n\ge 1,\\ & \leftarrow \text{Faulty generalization 2},\end{array}$$

(2)

where $\mathrm{sgn}(\cdot )$ is the sign function. We show that (2) does not hold in general for dynamic models, as Figure 1(b) illustrates. In this example, even though E $(S_n)$ decreases, ${\mu }_n$ first increases and then decreases, following an inverted-U shape. We list conditions under which weaker versions of (2) hold in Section 4.3. For example, in Corollary 1, we show that if mean service times are convex decreasing and an additional condition holds, then service rates are increasing. We expect our monotonicity results to help empirical researchers identify hidden changes in service rates that may correspond to new mechanisms.

Third, we show that for single-server dynamic models, the formula for the translation of mean service times to service rates differs greatly if the workload is measured at the end instead of the beginning of the service. Furthermore, for given service rates, a change in measurement timing can reverse the monotonicity pattern of the mean service times. This finding indicates that service times for which the state is measured at different times are not directly comparable. In general, the measurement timing choice should depend on data availability and hypothesized mechanisms. Some papers (such as Chan et al., 2012) propose models consistent with specific workload measurement timing.

Moreover, our results for dynamic models contribute to the analytical literature and expand the scope of previously published managerial insights. Recent work that is grounded in empirical findings formulates dynamic models to provide operational suggestions, and either assumes known service rates (Do et al., 2018) or estimates service rates directly (Cho et al., 2019). In contrast, we demonstrate how to translate empirical mean service times into service rates. Our results facilitate operational improvements (such as selecting staffing levels or choosing between one pooled queue versus parallel dedicated queues) through applications of analytical models calibrated with empirical data. Given that empirical findings are typically expressed in mean service times, the ability to translate service time data into service rate estimates provides the missing puzzle piece. To summarize, we make the following key contributions:

1. (1)

We categorize and name two major state-dependent queueing model types, static and dynamic models, and connect them to corresponding behavioral mechanisms. The choice of model should depend on both the hypothesized mechanisms and whether the model is consistent with data.

(2)

Using numerical examples, we show that although translating empirical service times into static models is always possible, the translation process may fail for dynamic models, thus providing evidence against dynamic mechanisms. In cases where translation is possible for both model types, performance measures such as mean waiting times can differ drastically.

(3)

For dynamic models, we provide closed-form expressions to compute service rates from mean service times. These expressions demonstrate that (1) does not hold in general. Also, we provide exact formulas to compute mean service times from service rates and show that service times follow phase-type (PH) distributions.

(4)

Contrary to what appears to be implicitly assumed in some of the empirical literature, for dynamic models, we show that monotone mean service times do not guarantee monotone service rates. We provide conditions that, combined with monotonicity for one sequence (service rate or mean service time), guarantee monotonicity of the other sequence, thus showing that (2) holds in certain situations, but not in general.

(5)

In dynamic models, with given service rates, we show that whether the workload is measured at the beginning or end of service leads to drastically different patterns of mean service times, both in terms of the magnitude and the directional changes.

The remainder of the paper is organized as follows: Section 2 reviews related literature; Section 3 formulates and compares the models of interest, and states assumptions; Section 4 provides results for dynamic models; Section 5 provides results for the case in which workload is measured at the end of service; Section 6 presents case studies and numerical examples; and Section 7 concludes. An electronic companion provides a table of notation, all proofs, and other supplemental materials.

2. Related Literature

Our work is related to both empirical work on how service times depend on workload and analytical work on the formulation and analysis of state-dependent queueing models. We review these two streams separately.

We begin by defining general terms that are needed to explain how “workload” is defined in different streams of the literature:

$$\begin{array}{rl}N(t)& =\text{number of customers in the system at time }t,\\ c(t)& =\text{number of servers at time }t,\\ B(t)& =min(N(t),c(t))=\text{number of busy servers at time }t,\\ U(t)& =\text{unfinished work at time }t,\\ W(t)& =\text{virtual waiting time at time }t.\end{array}$$

The analytical models that we discuss assume a constant number of servers $c$ , but some of the empirical work involves settings where the number of servers varies with time. The unfinished work is the sum of the service times of all waiting customers and the residual service times of all customers currently in service. The virtual waiting time is the time that a customer arriving at time $t$ would have to wait until commencing service. The unfinished work equals the virtual waiting time ( $U(t)=W(t)$ ) for single-server but not for multi-server queues (Sonderman, 1979).

Empirical Studies: We focus our discussion on three issues: (a) how the workload is defined, (b) when the workload is measured, and (c) what pattern in the mean service times is revealed. Table 1 compares a sample of recent empirical papers in terms of these three issues.

Table 1.

Empirical studies that involve the impact of workload on service time.

Citation	Workload at $t$	Time $t_i$	Service time pattern
Kc and Terwiesch 2009 ${}^{(1)}$	$B(t)/c(t)$	Entry to system	$\searrow$
	$B(t)/c(t)$	Average during service	$\searrow$
Kc and Terwiesch 2012 ${}^{(2)}$	$f(B(t)/c(t))$	End of service	$\searrow$
Tan and Netessine 2014	$N(t)$	Average during service	$\nearrow \searrow$
Batt and Terwiesch 2017	$N(t)-B(t)$	Beginning of service	$\nearrow \searrow$
	$B(t)$	Beginning of service	$\nearrow \searrow$
Berry Jaeker and Tucker 2017 ${}^{(3)}$	$f(B(t))$	Beginning and end of service	$\nearrow \searrow \nearrow$
Wang and Zhou 2018	$N(t)-B(t)$	Beginning of service	$\searrow$
Deo and Jain 2019 ${}^{(4)}$	$f(N(t))$	Beginning of service	$\nearrow$
Rastpour et al. 2020	$B(t)$	Beginning of service	$\nearrow \searrow$

${}^{(1)}$ This study also reports the impact of “overwork” as $\nearrow$ , but we focus only on the impact of workload.

${}^{(2)}f(\cdot )$ : Indicator function for whether occupancy is above a threshold.

${}^{(3)}f(\cdot )$ : Quadratic and spline functions.

${}^{(4)}f(\cdot )$ : Anticipated number of customers in the next period.

• (a)

How is the workload defined? In the empirical studies of interest to us, regression models are used to estimate mean service times, $\mathrm{E}(S_n)$ (or $\ln \mathrm{E}(S_n)$ ), as a function of workload and other variables. Possible ways to define workload include the proportion of servers that are busy ( $B(t)/c(t)$ ), the number of customers in the system ( $N(t)$ ), the number of customers in the queue ( $N(t)-B(t)$ ), or functions involving one or more of the aforementioned choices, and all of these definitions have been used in empirical studies (see Table 1, Column 2). The appropriate way to define workload in an empirical study depends on data availability, context, and the research question. It is important to consider which aspects of the workload are “visible” to servers, accurately and in real-time. For example, in a hospital intensive care unit (ICU), the number of busy servers (i.e., the number of occupied beds) is visible, but the queue length (i.e., the number of patients waiting to be admitted to the ICU) is ambiguous and might not be visible to those working in the ICU (Berry Jaeker and Tucker, 2017; Kc and Terwiesch, 2009, 2012). Therefore, servers can only react to changes in the number of busy servers, but not to changes in the queue length. In contrast, in a grocery store checkout lane, the number of waiting customers is visible to the server and can have an important impact on their behavior (Wang and Zhou, 2018).

(b)

When is the workload measured? Empirical researchers need to decide the time $t_i$ at which the workload relevant to the service time of customer $i$ is measured. Possibilities for $t_i$ include the time of entry to the system, the time at which service began, the time at which service ended, or an average of the workload from the beginning to the end of service (see Table 1, Column 3). The appropriate choice depends on the data availability and the research question. Arguments justifying the choice of $t_i$ are not always provided. Kc and Terwiesch (2012) provide an example of a well-justified choice for $t_i$ : workload should be measured at or near the end of an ICU patient’s length of stay (LOS) because the way in which the LOS might be impacted by workload is through the discharge of the “least-sick patient” in order to make room for a new admission. We show in Section 6.5 that for a given pattern of service rates, the pattern of corresponding mean service times can vary greatly depending on the time at which the workload is measured.

(c)

What pattern in the mean service times is revealed? Hypotheses for empirical studies of the impact of workload on server performance are typically expressed in terms of high-level features: for example, the relationship is increasing, decreasing, U-shaped, or inverted-U-shaped. The dependent variable is typically the service time for a particular customer. We show the patterns in Table 1 (Column 4), with $\nearrow$ denoting that mean service times increase with workload, $\searrow$ denoting that mean service times decrease with workload, and a sequence of these symbols denoting a non-monotone pattern. Early studies anticipated finding generalizable patterns, but as the field has matured, researchers have come to recognize that workload can impact service times through a wide variety of mechanisms and the research focus has shifted to investigating individual mechanisms (Delasay et al., 2019).

Analytical Studies: Queueing theorists have studied state-dependent queueing models since the 1960s. Such models can potentially capture the empirically observed ways in which server performance depends on workload. As we discussed in Section 1, many of these models fall into one of the two categories: static and dynamic models.

Table 2 summarizes a sample of papers on state-dependent queueing models. Some models do not fit neatly into our static and dynamic model categorization. For example, in Abouee-Mehrizi and Baron (2016), a static initial service time requirement distribution and a dynamic varying server speed determine service times. Similarly, Bekker et al. (2004), Bekker and Borst (2006), Bekker and Boxma (2007), and Sakuma et al. (2021) present models with both static and dynamic components.

Table 2.

Analytical studies of state-dependent queueing models.

Citation	Reported service process	Number of servers	Service time/rate depends on	Model type
Conway and Maxwell (1962)	$M$	1	$N(t)$	Dynamic
Jackson (1963)	$M$	$c$	$N(t)$	Dynamic
Hillier et al. (1964)	$M$	$c$	$N(t)$	Dynamic
Gebhard (1967)	$M$	1	$N(t)$	Dynamic
Hadidi and Conolly (1969)	$M$	1	$N(t)$	Dynamic
Mandelbaum and Pats (1995) ${}^{(1)}$	$M$	1, $c$ , $\mathrm{\infty }$	$N(t)$	Dynamic
George and Harrison (2001)	$M$	1	$N(t)$	Dynamic
Chan et al. (2014)	$M$	$c$	$N(t)$	Dynamic
Dong et al. (2015)	$M$	$c$	$N(t)$	Dynamic
Do et al. (2018)	$M$	$c$	$N(t)$	Dynamic
Chan et al. (2017)	$M$	$c$	$N(t)$	Static
D’Auria et al. (2022)	$M$	$c$	$W(t)$	Static
Harris (1967)	$G$	1	$N(t)$	Static
Schäl (1971)	$G$	1	$N(t)$	Static
Shanthikumar (1979)	$G$	1	$N(t)$	Static
Gray et al. (1992)	$G$	1	$N(t)$	Static
Gupta and Srinivasa (1998)	$G$	1	$N(t)$	Static
Bekker et al. (2004)	$G$	1	$U(t)$
Bekker and Borst (2006)	$G$	1	$U(t)$
Bekker and Boxma (2007)	$G$	1	$U(t)$
Abouee-Mehrizi and Baron (2016) ${}^{(2)}$	$G$	1	$N(t)$	Static
Oz et al. (2017)	$G$	1	$N(t)$	Static
Baron et al. (2018)	$G$	1	$N(t)$	Static
Sakuma et al. (2021)	$G$	1	$U(t)$

Note: For some papers, we describe a special case of a more general model analyzed in the paper.

${}^{(1)}$ This study investigates multiple models that differ only in terms of the number of servers.

${}^{(2)}$ The model is a static model only if the rescaling parameters are set to 1 for all states.

Model Properties: In all of the dynamic models in the literature, the process $\{N(t),t\ge 0\}$ is Markovian. Therefore, the model is relatively easy to solve. However, the Markovian process does not imply that service time distributions are exponential, as noted by Harris (1967). As we show in Section 4.2.2, service times in dynamic models follow infinite-state PH distributions. In static models, $N(t)$ is not Markovian, even for exponential service times, such as in Chan et al. (2017) and D’Auria et al. (2022). Therefore, analyses of static models are more complicated, may involve other state variables ( $U(t)$ or $W(t)$ ), and are typically limited to single-server models. The two multi-server static models listed in Table 2 (Chan et al., 2017; D’Auria et al., 2022) assume exponential service times as well as other limiting assumptions.

Translation: For static models, empirical service time data can be used directly to estimate state-dependent service time distributions, and thus specify the service process parameters. For dynamic models, a customer’s service time is the sum of a random number of sojourn times in system states visited during the service. To specify the service process in such models, one needs to translate empirical service times to sojourn time parameters.

None of the dynamic model papers in Table 2 investigate service time distributions. In empirical studies, as illustrated by the papers in Table 1, service times are typically the variables that are used to evaluate server performance. This creates a disconnect, which we address by developing methods to translate empirical service times into dynamic models, thus enabling empiricists to compare and evaluate both types of models. Another disconnect is that many empirical papers define workload in terms of $B(t)$ , whereas analytical papers allow the service process to depend on $N(t)$ , $U(t)$ , or $W(t)$ . We focus on models with a dependence on $N(t)$ , because $B(t)=min(N(t),c)$ is a simple function of $N(t)$ . Consequently, translating empirical findings that involve $B(t)$ into a model that allows $N(t)$ dependence presents no difficulties.

3. Models and Notation

In this section, we start by outlining notational conventions. In Sections 3.1 to 3.4, we compare specific static and dynamic models in terms of their formulations, performance evaluation, service time definitions, and corresponding mechanisms. We state assumptions in Section 3.5.

Figure 2.

Transition rate diagrams: (a) $M/M_n/c$ underlying chain, $N(t)$ and (b) $M/M_n/c$ absorbing chain, $X(t)$ .

We use uppercase bold letters to represent matrices and lowercase bold letters to represent vectors. The $(i,j)$ th element of matrix $\mathbf{A}$ is $A(i,j)$ and the $i$ th element of vector $\mathbf{a}$ is $a(i)$ . $\mathbf{I}$ is an identity matrix, $\mathbf{e}$ is a vector of ones, $\mathbf{0}$ is a vector of zeros, and ${\mathbf{u}}_i$ is the $i$ th unit vector (i.e., $u_i(i)=1,u_i(j)=0$ for $j\ne i$ ), all of appropriate size. $I_{(C)}$ is an indicator function for condition $C$ . The first-order and second-order forward differences for the sequence $\{x_n\}$ are $\mathrm{\Delta }{x}_n=x_{n+1}-x_n$ and ${\mathrm{\Delta }}^2{x}_n=\mathrm{\Delta }(\mathrm{\Delta }{x}_n)=\mathrm{\Delta }{x}_{n+1}-\mathrm{\Delta }{x}_n$ , $n\ge 0$ , respectively. We call a sequence $\{x_n\}$ convex if it is discretely convex (see Definition 1 in Zacharias and Pinedo, 2017). We use increasing and decreasing in their non-strict sense unless noted otherwise.

We assume throughout that there is a single queue, a single customer class, and a constant number $c$ of identical parallel servers that serve customers (who never abandon) from the queue in first-come-first-served order. Servers do not idle if there are customers waiting.

3.1. Model Formulations

In both the static and dynamic models, let $S_n$ be the service time random variable for a focal customer for whom the system state is $n$ immediately after service begins, that is, the workload is measured at the beginning of service.

The static model we focus on assumes exponential service times ( $M/M_n*/c$ ). It has a Poisson arrival process with rate $\lambda$ . The service time $S_n\sim \exp ({\mu }_n)$ for a customer whose service commences at time $t$ where $N(t^{+})=n$ . Previous literature uses the Kendall notation $M/M_n/c$ for this model, but we will show that static and dynamic models differ in many aspects, and using the same notation can be misleading. Therefore, we add an asterisk to distinguish the static model and denote it by $M/M_n*/c$ .

The dynamic model we focus on assumes exponential service rates ( $M/M_n/c$ ). Let the arrival rate be $\lambda$ and set $N(t)=n$ . We use $b_n=min(n,c)$ to denote the number of busy servers in state $n$ . The probability that a particular server completes service in state $n$ , $n\ge 1$ , in the infinitesimal interval $[t,t+\delta )$ is ${\mu }_n\delta I_{(\text{server is busy})}+o(\delta )$ , which implicitly defines the service rate ${\mu }_n$ in state $n$ .

Note that the natural way to specify a static model is through the parameters of the state-dependent service time distributions, for example, their means. In contrast, the natural way to specify a dynamic model is through state-dependent service rates.

3.2. Performance Evaluation

Performance evaluation is the process of evaluating performance measures (such as mean waiting time) for a queueing model with specified input parameters. Common approaches for performance evaluation include simulation and obtaining analytical solutions.

It is easy and natural to obtain performance measures in $M/M_n*/c$ using simulation because a customer’s service time only needs to be generated once, at the beginning of service. After that, the service completion time becomes available for this particular service, and the system’s next event (a customer arrival or departure) is determined. However, obtaining its analytical results is difficult. Only boundary results for multi-server cases are provided (Chan et al., 2017) and exact solutions are only available for single-server cases (such as in Oz et al., 2017). Even for single-server cases, no published papers provide ready-to-use numerical algorithms to evaluate performance from the analytical solutions. The analytical challenge is that the process $\{N(t),t\ge 0\}$ is not Markovian. To formulate the model as a Markov process, one needs to add state variables to keep track of the workload at the beginning of service for each busy server.

In contrast, it is easy and natural to obtain analytical solutions for $M/M_n/c$ : The process $\{N(t),t\ge 0\}$ is a time-homogeneous continuous-time Markov chain (CTMC), a birth-death process with the birth rate $\lambda$ and death rate $b_n{\mu }_n$ in state $n$ , as illustrated in Figure 2(a). We refer to this infinite-state CTMC as the “underlying chain” for $M/M_n/c$ . While obtaining analytical results for $M/M_n/c$ is relatively easy, simulating it is significantly harder than $M/M_n*/c$ . One natural way to simulate $M/M_n/c$ is to generate sojourn times between transitions. Once a transition occurs, we use independent random numbers to determine the type of event (an arrival or a departure) and which customer leaves the system (if the event is a departure). Parameters for sojourn times are updated after each transition. To collect service times, we need to sum up all sojourn times during a service.

3.3. Service Time Definitions

In $M/M_n*/c$ , the distribution of the service time $S_n$ is defined explicitly. Translating $S_n$ into model specifications is straightforward: $S_n\sim \exp ({\mu }_n)$ and $\mathrm{E}(S_n)=1/{\mu }_n$ .

In contrast, in $M/M_n/c$ , the distribution of $S_n$ is implicitly defined, through sojourn times distributions, and must be obtained through model analysis. $S_n$ is the time to absorption in a modified version of the underlying chain, which we call the “absorbing chain,” $\{X(t),t\ge 0\}$ , as discussed in detail in Section 4.1 and illustrated in Figure 2(b). The absorbing chain begins in state $n$ , $n\ge 1$ at $t=0$ , corresponding to the beginning of service for the focal customer, and $S_n$ is the time to reach the absorbing state $a$ , corresponding to the completion of service. Therefore, $S_n$ is a first-passage time from $n$ to $a$ in the absorbing chain, defined as:

$$S_n=\text{inf }\{t^{\mathrm{\prime }}:X(t^{\mathrm{\prime }})=a|X(0)=n\},n\ge 1.$$

(3)

$S_n$ equals the sum of a random number of independent and exponentially distributed sojourn times, that is, times from one instant at which $N(t)$ changes to the next. Thus, the total duration of service depends on stochastic events (such as arrivals or service completions of non-focal customers) that occur after the focal customer’s service has begun.

Therefore, in $M/M_n/c$ , translating $S_n$ into ${\mu }_n$ is more complicated, as shown in Section 4.1, and $S_n$ does not follow an exponential distribution. Instead, it follows an infinite-state PH distribution, as discussed in Section 4.2.2. We define $s_n\equiv \text{E}(S_n)$ , to simplify the notation in the remainder of the paper.

We revisit the dynamic $M/M_n/1$ model in Section 5, with $S_n^l$ defined as the service time for a focal customer for whom the system state is $n$ immediately before service ends, with $s_n^l\equiv \text{E}(S_n^l)$ . That is, the workload is measured at the end of service.

3.4. Models and Mechanisms

$M/M_n*/c$ and $M/M_n/c$ imply different server behavior. To explain the differences, we express service time as work content divided by server speed (as in Delasay et al., 2019). For example, the work content and server speed for supermarket checkout could be approximated by the number of items to scan and items scanned per time unit, respectively. Depending on how work content and speed are affected by $n$ , we identify two groups of mechanisms:

3.4.1. Dynamic Mechanisms

The remaining work content or the server speed changes during the service.

3.4.2. Static Mechanisms

Both the work content and the server speed are fixed during the service.

Dynamic mechanisms include social speedup pressure and task reduction. Static mechanisms include forgetting and customer early task initiation (see Table 2 in Delasay et al., 2019, for definitions.) By fitting both models to data, one might find evidence against the dynamic $M/M_n/c$ and thus against dynamic mechanisms. We illustrate this in Section 6.2.

3.5. Assumptions

We invoke the following assumptions for some of our analyses.

Assumption 1

$\lambda$ , ${\mu }_n$ , $s_n$ , and $s_n^l$ , for $n\ge 1$ , are positive and finite real numbers, and $c$ is a positive and finite integer.

The requirements that ${\mu }_n\in (0,\mathrm{\infty })$ and $\lambda \in (0,\mathrm{\infty })$ are standard for $M/M_n/c$ . Specifically, ${\mu }_n<\mathrm{\infty }$ and $\lambda <\mathrm{\infty }$ avoid instantaneous transitions, and $0<{\mu }_n$ and $0<\lambda$ ensures that $N(t)$ is irreducible. Similarly, we assume $s_n\in (0,\mathrm{\infty })$ and $s_n^l\in (0,\mathrm{\infty })$ to avoid technical difficulties.

We invoke Assumption 1 for all $M/M_n/c$ derivations and theoretical results. Some formulas we derive can be used for calculations, however, even if Assumption 1 fails to hold. We discuss this in Remark 1 in Section 4.

Assumption 2

The service rates converge at state $m$ ( $c\le m<\mathrm{\infty }$ ), that is: ${\mu }_{m+i}=\mu$ for $i\ge 0$ .

This assumption is invoked for certain $M/M_n/c$ derivations. As a special case, if $m=c$ , ${\mu }_n=\mu$ for $n\ge c$ and the service rates in the $M/M_n/c$ model depend only on the number of busy servers at time $t$ , $B(t)$ . Defining workload as $B(t)$ or a function of $B(t)$ is common in the empirical literature (see Table 1).

Assumption 3

We set $s_0=0$ .

This assumption, which is always in effect, is made to simplify the statement of certain equations.

4. Dynamic Model: $M/M_n/c$

Since translating $S_n$ into $M/M_n*/c$ is straightforward (see Sections 3.1 and 3.3), we hereafter focus on translation issues in dynamic $M/M_n/c$ . In this section, we provide a series of results: formulas to compute ${\mu }_n$ given $s_n$ in Section 4.1; formulas to compute $s_n$ given ${\mu }_n$ and a PH distribution analysis of $S_n$ for $M/M_n/1$ in Section 4.2; and monotonicity results regarding the ${\mu }_n$ and $s_n$ sequences in Section 4.3. Results for $s_n$ given ${\mu }_n$ and an infinite-state PH distribution analysis of $S_n$ in $M/M_n/c$ are in EC.3.2 and EC.3.4 of E-Companion, respectively.

4.1. Service Rates Given Mean Service Times

Empirical studies typically evaluate state-dependent server performance in terms of estimated mean service times (see Section 2). Should one wish to use a dynamic model, one needs to translate the mean service times, $s_n$ , into service rates ${\mu }_n$ . We provide a closed-form solution for this task, namely, for computing ${\mu }_n$ given $s_n$ and the arrival rate $\lambda$ , which we assume have been estimated using empirical data. We show that, in general, the equation  ${\mu }_n=1/s_n$ in (1) is a faulty generalization.

We use the absorbing chain (see Figure 2(b)) to analyze $s_n$ . As indicated in (3), $S_n$ is the time to absorption in the absorbing chain, starting in state $n,n\ge 1$ . The absorbing chain is obtained from the underlying chain (shown in Figure 2(a)) through the following modifications: (a) Remove state 0 (because service cannot begin in state 0) and add absorbing state $a$ (corresponding to service completion); (b) set the transition rate from $n$ to $a,n\ge 1$ , to the rate of service completion for the focal customer, that is, ${\mu }_n$ ; and (c) set the transition rate from $n$ to $n-1,n\ge 2$ , to the rate at which non-focal customers depart, that is, $(b_n-1){\mu }_n$ . The sojourn time $T_n$ in state $n$ is exponentially distributed with rate $b_n{\mu }_n+\lambda$ and its expected value $t_n\equiv \mathrm{E}(T_n)=1/(b_n{\mu }_n+\lambda )$ . We use first-step analysis (see Kulkarni, 2010, Section 3.1) to obtain a recursive equation for $s_n$ . This method involves conditioning on the next transition in the absorbing chain. Since we are interested in the next-transition probabilities, we use the embedded discrete-time Markov chain (DTMC) $\{X_k,k\ge 0\}$ of the absorbing chain $\{X(t),t\ge 0\}$ , where $X_k$ is the number of customers immediately after the $k$ th transition (state change) for $X(t)$ .

Before absorption, the next transition (“the first step”) for both $X(t)$ and $X_k$ must be one of the following: focal customer finishes service, new customer arrives, or non-focal customer departs. We denote these events and their probabilities as $F_n$ , $f_n$ , $A_n$ , $a_n$ , $D_n$ , and $d_n$ , respectively, where

$$\begin{array}{rl}{f}_n& =Pr(F_n)=\frac{{\mu }_n}{b_n{\mu }_n+\lambda },a_n=Pr(A_n)=\frac{\lambda }{b_n{\mu }_n+\lambda },\\ d_n& =Pr(D_n)=\frac{(b_n-1){\mu }_n}{b_n{\mu }_n+\lambda }.\end{array}$$

Therefore, using the Markov property:

$$\begin{array}{rl}{s}_n& =Pr(F_n)\mathrm{E}(S_n|F_n)+Pr(D_n)\mathrm{E}(S_n|D_n)+Pr(A_n)\mathrm{E}(S_n|A_n)\\ & =\frac{{\mu }_n}{b_n{\mu }_n+\lambda }{t}_n+\frac{(b_n-1){\mu }_n}{b_n{\mu }_n+\lambda }(t_n+s_{n-1})+\frac{\lambda }{b_n{\mu }_n+\lambda }(t_n+s_{n+1})\\ & =\frac{1}{b_n{\mu }_n+\lambda }+\frac{(b_n-1){\mu }_n}{b_n{\mu }_n+\lambda }{s}_{n-1}+\frac{\lambda }{b_n{\mu }_n+\lambda }{s}_{n+1},n\ge 1,\end{array}$$

(4)

$$\begin{array}{rl}{s}_n& =t_n+d_n{s}_{n-1}+a_n{s}_{n+1},n\ge 1,\end{array}$$

(5)

with $s_0=0$ by Assumption 3.

By solving (4) for ${\mu }_n$ , we obtain the following expression for $M/M_n/c$ :

$${\mu }_n=\frac{\lambda s_{n+1}-\lambda s_n+1}{b_n{s}_n-(b_n-1)s_{n-1}}=\frac{\lambda \mathrm{\Delta }{s}_n+1}{(b_n-1)\mathrm{\Delta }{s}_{n-1}+s_n},n\ge 1.$$

(6)

For $M/M_n/1$ , $b_n=1$ and $d_n=0$ , and therefore (5) and (6) reduce to:

$$\begin{array}{rl}{s}_n& =t_n+a_n{s}_{n+1},n\ge 1,\end{array}$$

(7)

$$\begin{array}{rl}{\mu }_n& =\frac{\lambda \mathrm{\Delta }{s}_n+1}{s_n},n\ge 1.\end{array}$$

(8)

In empirical work, it is common to assume that service times depend only on $B(t)$ (see Table 1), which implies that $s_n=s_c$ for $n\ge c$ . For $c=1$ , this implies that all service rates are equal. For $c>1$ , this implies that faulty generalization (1) holds true in the tail, that is:

$${\mu }_n=\frac{1}{s_c},n\ge c+1.$$

Remark 1

We implicitly invoked Assumption 1 in deriving (6), as we do for all of our derivations. However, equation (6) can be used to calculate ${\mu }_n$ values even if Assumption 1 does not hold. Such calculations could lead to negative or infinite ${\mu }_n$ values. It might not be known a priori whether Assumption 1 holds, for a set of $s_n$ values that have been estimated as part of an empirical study. Calculating ${\mu }_n$ values using (6) helps one to check whether the assumption holds.

Based on (6), (8), and Lemma EC.3 of E-Companion, we see that the numerator of ${\mu }_n$ is non-positive if $\mathrm{\Delta }{s}_n\le -1/\lambda$ , which implies that ${\mu }_n$ will be non-positive if we try to fit an $M/M_n/c$ model. In other words, a sharp decrease from $s_n$ to $s_{n+1}$ can cause ${\mu }_n$ to become negative, and hence, invalid—which we interpret as evidence that a data set is inconsistent with the dynamic model. The corresponding example in Section 6.2 illustrates this.

Figure 3.

The absorbing chain and the Markov jump process, for $c=1$ , under Assumption 2: (a) absorbing chain and (b) absorbing Markov jump process.

Figure 4.

Simulated mean waiting times with 95% confidence intervals.

4.2. Service Times Given Service Rates

Some empirical papers (such as Cho et al., 2019) estimate state-dependent service rates directly from data. Therefore, computing mean service times ( $s_n$ ) from estimated service rates ( ${\mu }_n$ ) is of interest, and we provide formulas for this purpose in this subsection. An estimate of the arrival rate ( $\lambda$ ) is needed to use the formulas. In addition to the mean, we characterize the distribution of the service time random variable. We provide these results for the $M/M_n/1$ model in the main text and analogous but more complicated results for the multi-server $M/M_n/c$ model in EC.3.2 and EC.3.4 of E-Companion.

4.2.1. Mean Service Times Given Service Rates

If $s_{n-1}$ and $s_{n+1}$ are known, then (7) can be used to obtain $s_n$ . This is not sufficient, however—we also need to obtain an initial value, such as $s_1$ , which is a challenging task. The correct solution to (7) is the minimal non-negative solution (see Lemma EC.1 of E-Companion). Our analysis provides the correct value for $s_1$ , which together with (7) can be used to obtain the remaining $s_n$ values in $M/M_n/1$ . We begin by defining the following sequence:

$${\alpha }_n=\prod _{i=1}^{n-1}{a}_i,n\ge 2,{\alpha }_1=1.$$

(9)

Theorem 1

If $c=1$ then:

$$s_n=\sum _{i=n}^{\mathrm{\infty }}\frac{{\alpha }_i}{{\alpha }_n}{t}_i,n\ge 1.$$

(10)

Furthermore, under Assumption 2,

$$s_n=\{\begin{array}{ll}{\displaystyle \sum _{i=n}^m\frac{{\alpha }_i}{{\alpha }_n}{t}_i+\frac{{\alpha }_{m+1}}{{\alpha }_n(1-a_m)}{t}_m,}& 1\le n\le m,\\ {\displaystyle \frac{1}{{\mu }_m},}& n\ge m+1.\end{array}$$

(11)

Wang and Zhou (2018) previously provided the infinite-series representation in (10). We show that under Assumption 2, one obtains a closed-form expression involving only a finite sum, as in (11). Theorem 1 cannot be generalized directly to $M/M_n/c$ with $c\ge 2$ . The analysis for that model is more complex. We present $M/M_n/c$ results for $s_n$ with Theorem EC.1 in EC.3.2 of E-Companion.

Figure 5.

Estimated ${\mu }_n$ from (30) using a perturbed $s_n$ vector.

4.2.2. Formulating Service Time as a PH Random Variable

Service times in $M/M_n*/1$ and non-state-dependent $M/M/1$ models are exponentially distributed; in contrast, as we show in this section, service times in state-dependent $M/M_n/1$ models are PH distributed.

As explained in Section 4.1, we formulate $S_n$ for $M/M_n/1$ as the time to absorption in the absorbing chain $X(t)$ . So far, we have only studied the first moment, but the theory of PH distributions permits us to compute higher moments and the full distribution of $S_n$ . Most of the theory of PH distributions assumes a Markov chain with a finite number of states, but the absorbing chain for $M/M_n/1$ has an infinite number of states. Under Assumption 2, we can aggregate states in the tail, resulting in a finite-state Markov chain, for which the time to absorption provides the exact distribution of $S_n$ .

The time to absorption in a finite-state CTMC with a single absorbing state and all other states transient is, by definition, a PH random variable (Latouche and Ramaswami, 1999, Definition 2.3.1). The canonical form of the infinitesimal generator matrix (“generator” for short) is:

$$\begin{array}{r}\mathbf{Q}=\begin{array}{c}a\\ T\end{array}\stackrel{aT}{\left[\begin{array}{cc}0& \mathbf{0}\\ {\mathbf{l}}^a& \mathbf{L}\end{array}\right]}\end{array}$$

(12)

where $T$ is the set of transient states, $a$ is the absorbing state, $\mathbf{L}$ is a subgenerator matrix of transition rates among the transient states, and ${\mathbf{l}}^a$ is a column vector of absorption rates from the transient states. The row vector $\mathit{\beta }$ specifies the initial probabilities for $T$ . A PH distribution is represented by the notation PH $(\mathit{\beta },\mathbf{L})$ . The parameters $\mathit{\beta }$ and $\mathbf{L}$ suffice to fully specify the absorbing chain, because ${\mathbf{l}}^a$ must equal $-\mathbf{L}\mathbf{e}$ (so that the rows of $\mathbf{Q}$ sum to zero), and the initial probability for state $a$ must equal $1-\mathit{\beta }\mathbf{e}$ (so that the initial probabilities to sum to one). EC.3.3 of E-Companion provides the nonzero elements of (12) and $\mathit{\beta }.$

Under Assumption 2, the original absorbing chain in Figure 2(b) simplifies for $c=1$ in that only two transitions are possible from state $n\ge m$ : To $n+1$ with probability $a_m$ and to $a$ with probability $f_m$ (see Figure 3(a)). Therefore, each transition is a Bernoulli trial with a constant success probability $f_m$ , and the number of transitions until service finish is geometrically distributed. This observation allows us to lump all states $n\ge m$ together, as shown for state $m$ in Figure 3(b).

We modify the absorbing chain (Figure 3(a)) by allowing a virtual transition (a self-transition) in state $m$ , resulting in the transition diagram in Figure 3(b). We follow Kao (1996: Example 5.6.3) in referring to the corresponding process as a continuous-time Markov jump process (MJP; a CTMC with virtual transitions). We model $S_n$ as the time to absorption in the absorbing MJP in Figure 3(b), where state $a$ is the absorbing state and all other states, $\{1,2,\dots ,m\}$ , are transient. The arrival and service completion probabilities in state $m$ in the absorbing chain are $a_m=\lambda /(\mu +\lambda )$ and $f_m=\mu /(\mu +\lambda )$ , respectively; the former corresponds to a virtual transition in the absorbing MJP. Hence, we obtain:

$$L(m,m)=f_{m}L(m,m)=\frac{\mu }{\mu +\lambda }\times (-(\mu +\lambda ))=-\mu ,$$

(13)

and the remaining non-zero elements are the same as in EC.3.3 of E-Companion. If $n\le m$ , then $S_n$ is distributed as PH( $\mathit{\beta },\mathbf{L}$ ) and PH distribution formulas for the cumulative distribution function and moments apply (see EC.32 to EC.34 of E-Companion). If $n>m$ , then $S_n$ has the same distribution as $S_m$ , which is exponential with rate $\mu$ , because states $n\ge m$ in the absorbing chain are aggregated into state $m$ in the absorbing MJP.

This approach only works for $c=1$ . For $c>1$ in state $n\ge m$ , in addition to transitions to $a$ and $n+1$ , a transition to $n-1$ is possible, with probability $d_m$ (compare Figures 2(b) and 3(a)). Consequently, the Bernoulli success probability can change, if the process moves below state $m$ —something that was not possible for $c=1$ . Therefore, the states above $m$ are not identical and cannot be lumped together for $c>1$ . As a result, the service time $S_n$ does not follow a standard PH distribution because the related absorbing chain has an infinite number of states. However, we show with Theorem EC.2 in EC.3.4 of E-Companion that the distribution of $S_n$ in $M/M_n/c$ with $c\ge 1$ belongs to the SPH class (an infinite-state generalization of PH; Shi et al., 2005).

4.3. Monotonicity Properties

We have seen (Figure 1(b)) that the direction of change in service rates is not always the opposite of the direction of change for mean service times, that is, (2) is indeed a faulty generalization. In this section, we present a series of results that are weaker than (2), but correct. We present a theorem with an if-and-only-if relationship between ${\mu }_n$ and $s_n$ ; and a corollary that shows how this relationship simplifies for $c=1$ . We present three additional corollaries with if–then relationships, stating that if one sequence is monotone, for $n\ge 1$ , possibly combined with additional conditions, then the other sequence is monotone in the opposite direction, for $n\ge 1$ .

Theorem 2

$${\mu }_n\text{increases}:\mathrm{\Delta }{\mu }_n\ge 0,n\ge 1,$$

(14)

if and only if:

$$s_n\text{decreases}:\mathrm{\Delta }{s}_n\le 0,n\ge 1,$$

(15)

and the following additional condition holds:

$${\mathrm{\Delta }}^2{s}_n\ge (\mathrm{\Delta }{s}_n+\frac{1}{\lambda })\left(\frac{b_{n+1}\mathrm{\Delta }{s}_n-(b_n-1)\mathrm{\Delta }{s}_{n-1}}{b_n\mathrm{\Delta }{s}_{n-1}+s_{n-1}}\right),n\ge 1.$$

(16)

An intuitive explanation for (16) is that it sets a lower bound for ${\mathrm{\Delta }}^2{s}_n$ , the rate of change for $\mathrm{\Delta }{s}_n$ . Although (15) and (16) resemble the conditions for convex decreasing $s_n$ , we show numerically in Section 6.4 that a non-convex (convex) decreasing $s_n$ may (may not) lead to increasing ${\mu }_n$ .

Corollary 1

If

$$\begin{array}{rl}& s_n\text{ decreases:}\mathrm{\Delta }{s}_n\le 0,n\ge 1,\end{array}$$

(17)

$$\begin{array}{rl}& \text{is convex: }{\mathrm{\Delta }}^2{s}_n\ge 0,n\ge 1,\end{array}$$

(18)

$$\begin{array}{rl}& \text{and this additional condition holds: }\\ & (b_{n+1}-1){\mathrm{\Delta }}^2{s}_n\le |\mathrm{\Delta }{s}_{n+1}|,n\ge 0,\end{array}$$

(19)

then

$${\mu }_n\text{ increases: }\mathrm{\Delta }{\mu }_n\ge 0,n\ge 1.$$

(20)

Condition (19) sets an upper bound on the rate of change for $\mathrm{\Delta }{s}_n$ before the monotonicity of ${\mu }_n$ is lost.

Our next corollary requires no additional conditions.

Corollary 2

If ${\mu }_n$ increases ( $\mathrm{\Delta }{\mu }_n\ge 0,n\ge 1$ ) then $s_n$ decreases ( $\mathrm{\Delta }{s}_n\le 0,n\ge 1$ ).

Theorem 2 and Corollaries 1 and 2 hold for $c\ge 1$ . In the important single-server case, the conditions in Theorem 2 and Corollary 1 simplify considerably, as we show in Corollaries 3 and 4.

Corollary 3

For $c=1$ ,

$${\mu }_n\text{ increases:}\mathrm{\Delta }{\mu }_n\ge 0,n\ge 1,$$

(21)

if and only if:

$$\begin{array}{rl}& s_n\text{ decreases:}\mathrm{\Delta }{s}_n\le 0,n\ge 1,\end{array}$$

(22)

$$\begin{array}{rl}& \text{and this additional condition holds: }\\ & {\mathrm{\Delta }}^2{s}_n\ge (\mathrm{\Delta }{s}_n+\frac{1}{\lambda })\left(\frac{\mathrm{\Delta }{s}_n}{s_n}\right),n\ge 1.\end{array}$$

(23)

From Assumption 1, (22) and Lemma EC.3 of E-Companion, the right-hand side of (23) is guaranteed to be non-positive.

Corollary 4

For $c=1$ , if $s_n$ is convex decreasing ( $\mathrm{\Delta }{s}_n\le 0$ and ${\mathrm{\Delta }}^2{s}_n\ge 0,n\ge 1$ ) then ${\mu }_n$ increases ( $\mathrm{\Delta }{\mu }_n\ge 0,n\ge 1$ ).

Corollary 3 and our other results are more widely applicable than Proposition 5 in Wang and Zhou (2018: Appendix A). We discuss the differences in detail in EC.3.6 of E-Companion.

5. Single-Server Model With Workload Measured at the End of the Service

So far, we have assumed that workload is measured at the beginning of service—an approach that makes sense if one views the workload as a cause and the service duration as an effect. However, empirical studies vary in whether workload is measured at the beginning of service, the end of service, or in some other way (see Table 1). It is therefore important to investigate whether the relationships that we derived for dynamic models change substantially if we modify the workload measurement timing. To achieve this, we assume that the researchers have measured workload at the end of service and aim to translate the results into a dynamic $M/M_n/1$ model.

We find that the relationship between mean service times and service rates changes drastically, which suggests that the choice of how to measure workload could significantly impact estimated service times.

Notably, the results in Section 4, for workload measured at the beginning of service, do not require stability. Only the results in this section, for workload measured at the end of service, require stability.

We define the system state $n$ for the focal customer as the number of customers in the system immediately before the end of service. Specifically, we suppose that the service of the focal customer begins at time $t^{\mathrm{\prime }}$ and ends at time $t$ , for which $\underset{s\to t^{-}}{lim}N(s)=n$ ; we let $S_n^l=t-t^{\mathrm{\prime }}$ denote the service time of the focal customer; and we set $s_n^l=\mathrm{E}(S_n^l)$ , as stated in Section 3.3.

If the current state is $n$ , $n\ge 1$ , then the preceding (last) transition must have been either (a) the arrival of a customer to state $n-1$ , denoted $A_{n-1}^l$ ; or (b) the service completion of the customer preceding the focal customer, in state $n+1$ , denoted $F_{n+1}^l$ (meaning, after the preceding customer departed, the focal customer started service in state $n$ ). Therefore, by conditioning on the last event (last-step analysis), we obtain:

$$\begin{array}{rl}{s}_n^l& =Pr(F_{n+1}^l)\mathrm{E}(S_n^l|F_{n+1}^l)+Pr(A_{n-1}^l)\mathrm{E}(S_n^l|A_{n-1}^l)\\ & =Pr(F_{n+1}^l)t_n+Pr(A_{n-1}^l)(t_n+s_{n-1}^l)\\ & =t_n+Pr(A_{n-1}^l)s_{n-1}^l,n\ge 1,s_0^l=0.\end{array}$$

(24)

Next, to obtain $Pr(A_{n-1}^l)$ , we employ the DTMC $\{N_k,k\ge 0\}$ embedded to the CTMC for the $M/M_n/1$ system $\{N(t),t\ge 0\}$ . If the service completion of the focal customer is the $(k+1)$ th transition in the embedded DTMC, then:

$$\begin{array}{rl}Pr(A_{n-1}^l)& =Pr(N_{k-1}=n-1|N_k=n)\\ & =\frac{Pr(N_{k-1}=n-1,N_k=n)}{Pr(N_k=n)}\\ & =\frac{Pr(N_{k-1}=n-1)Pr(N_k=n|N_{k-1}=n-1)}{Pr(N_k=n)}\\ & =\frac{{\pi }_{n-1}{a}_{n-1}}{{\pi }_n}=f_n,\end{array}$$

(25)

where ${\pi }_{n-1}$ and ${\pi }_n$ are the steady state probabilities for states $n-1$ and $n$ , respectively, $a_{n-1}$ is the arrival probability for state $n-1$ , and $f_n$ is the finish probability. The last equation follows from that the local balance equation for states $n-1$ and $n$ is ${\pi }_{n-1}{a}_{n-1}={\pi }_n{f}_n$ . Note that (25) also follows from time reversibility for birth-death processes. Using (25) in (24), we obtain a recursive equation for $s_n^l$ as:

$$s_n^l=t_n+f_n{s}_{n-1}^l=\frac{1}{{\mu }_n+\lambda }+\frac{{\mu }_n}{{\mu }_n+\lambda }{s}_{n-1}^{l}n\ge 1,s_0^l=0.$$

(26)

Solving for ${\mu }_n$ from (26), we obtain:

$${\mu }_n=\frac{1-\lambda s_n^l}{\mathrm{\Delta }{s}_{n-1}^l},n\ge 1,s_0^l=0.$$

(27)

Equation (27) was derived under Assumption 1 but can be used for calculations even if Assumption 1 does not hold; see Remark 1.

Figure 6.

Faulty generalization (1) example: (a) comparison of service times and (b) comparison of service rates.

6. Case Studies and Numerical Examples

In this section, we illustrate how empirical researchers might translate empirical findings into queueing models, as well as the implications of some of our results. First, we provide a case study based on empirical findings to illustrate how different models can lead to widely different predicted mean waiting times. Second, we perturb the case study’s empirical findings to illustrate how translation to a dynamic model could fail, which provides evidence against dynamic mechanisms. Third, focusing on dynamic models, we use simulated data to illustrate faulty generalizations (1) and (2). Fourth, we provide an example in which convex decreasing mean service times do not correspond to increasing service rates and another example in which non-convex but decreasing mean service times correspond to increasing service rates. Fifth, we illustrate how changing the timing of workload measurement can reverse the mean service time trend.

6.1. Case Study 1: Comparing Model Predictions

Wang and Zhou (2018) estimate mean service times in a supermarket checkout setting using log-linear regression. Using results from that paper and additional information provided by the authors, after averaging over the effects of explanatory variables other than $n$ (see EC.4.1 of E-Companion), we obtain:

$${\hat{s}}_n=\exp (4.453-0.114n+0.010n^2)\text{seconds,}n\ge 1.$$

(28)

To translate the estimated mean service times into a queueing model, one must make an assumption about asymptotic behavior. If one extrapolates to large $n$ values, (28) indicates that ${\hat{s}}_n$ will eventually increase without bound, but Wang and Zhou (2018) report convex decreasing ${\hat{s}}_n$ , which we take to mean that ${\hat{s}}_n$ is convex decreasing within the range of their data. Therefore, we assume that ${\hat{s}}_n$ follows (28) until it reaches its minimum (at $n=6$ ), and remains constant thereafter, resulting in:

$$\begin{array}{r}{\hat{s}}_n=(77.4,71.2,66.7,63.9,62.4,62.1,62.1,\dots )\text{ seconds.}\end{array}$$

(29)

We translate these mean service times into three models: $M/M_n/1$ (the model that Wang and Zhou, 2018, assumed in their theoretical work), $M/M_n*/1$ , and $M/L{N}_n*/1$ (consistent with their argument for a lognormal service time distribution; see their Figures 5 and 6). We compare the mean waiting times (in seconds) for the three models, with respect to the number of dedicated checkout lanes, $i$ , which we vary from 3 to 10. The paper does not provide the total arrival rate, but we assume it to be $0.045$ customers per second, and that customers select a queue randomly, resulting in $i$ independent single-server queues, each with an arrival rate of $\lambda =0.045/i$ . For the static models ( $M/M_n*/1$ and $M/L{N}_n*/1$ ), we use ${\hat{s}}_n$ to specify the mean service times and the logged variance, equal to 0.184, to fully specify the log-normal service time distributions. For the dynamic model ( $M/M_n/1$ ), we use ${\hat{s}}_n$ and $\lambda$ to compute service rates, using (8).

We use simulation to estimate utilization rates and mean waiting times (see EC.4.2 of E-Companion for details). The utilization rates for the three models are similar: $96.5\mathrm{\%}$ for $M/M_n/1$ , $96.2\mathrm{\%}$ for $M/M_n*/1$ , and $96.9\mathrm{\%}$ for $M/L{N}_n*/1$ . However, Figure 4 shows that a single sequence of estimated mean service times (from (29)) results in noticeably different estimated mean waiting times across the three models, with $M/L{N}_n*/1$ showing lower mean waiting times. $M/M_n*/1$ and $M/M_n/1$ show similar mean waiting times, but that is not always the case, as we show in Section 6.2.

Which model is most appropriate? To answer that question, one would ideally obtain empirical data on mean waiting time (or some other performance measure)—information that is not available in this case.

6.2. Case Study 2: Translation Failure as Evidence Against Dynamic Mechanisms

Next, we illustrate how the translation of empirical results into a dynamic model can fail, providing evidence against dynamic mechanisms. We do this in a fully controlled setting, in which we abstract from the empirical mean service times in (29) to a data-generating process that is a static $M/M_n*/1$ model. Recall that this means that service times, conditional on starting state $n$ , are exponentially distributed with rate ${\mu }_n=1/s_n$ . For the purposes of the illustration, we set $s_1=200$ , $s_n={\hat{s}}_n$ from (29) for $n\ge 2$ , and $\lambda =0.010$ .

We begin by simulating data from this process (see EC.4.2 of E-Companion for details). Then we use the simulated data to estimate the arrival rate ( $\hat{\lambda }=0.010$ ) and the mean service times, using the following regression model:

$$S_n={\beta }_0+{\beta }_{1}n+{\beta }_2{n}^2+{\beta }_3{I}_{(n=1)}+ϵ,$$

(30)

which results in:

$$\begin{array}{rl}{\hat{s}}_n& =73.889-2.634n+0.114n^2+131.259I_{(n=1)}\\ & =(202.628,69.077,67.013,\dots ).\end{array}$$

(31)

Suppose that an empirical researcher attempts to translate the mean service times from (31) into an $M/M_n/1$ dynamic model. This would be reasonable, because the researcher would not know the true data-generating process but might know that the $M/M_n/1$ model has been studied extensively. The researcher might be surprised to find that using (8) to compute ${\hat{\mu }}_1$ results in a negative value:

$$\begin{array}{rl}{\hat{\mu }}_1& =\frac{\hat{\lambda }({\hat{s}}_2-{\hat{s}}_1)+1}{{\hat{s}}_1}=\frac{0.010\times (69.077-202.628)+1}{202.628}\\ & =-\mathrm{0.002.}\end{array}$$

Figure 5 illustrates these findings. One may wonder whether the negative ${\hat{\mu }}_1$ is statistically significant. We bootstrap model (30), using a “bootstrapping pairs” approach (see Efron and Tibshirani 1994, Section 5.2) to assess the robustness of negative ${\hat{\mu }}_1$ . We find that 100 bootstrap samples of ${\hat{\mu }}_1$ are all negative.

The negative ${\hat{\mu }}_1$ indicates that translation into $M/M_n/1$ is not possible. Within our fully controlled setting, we can draw a strong conclusion: the translation failed, because $M/M_n/1$ is an incorrect model. The true model is $M/M_n*/1$ . An empirical researcher would not be able to draw such a strong conclusion, because the researcher would not know the true model. However, mean service times estimated in an empirical study that results in negative service rates, using (8), provide evidence against $M/M_n/1$ , as well the dynamic mechanisms that are consistent with that dynamic model—even if that evidence might not be conclusive. To see why (8) results a negative ${\mu }_1$ value, note that within an $M/M_n/1$ model, for $s_1$ to be large, ${\mu }_1$ needs to be small, and $\lambda /(\lambda +{\mu }_1)$ (the probability that $N(t)$ moves from State 1 to State 2 before the service completion) needs to be small. In this case, it is impossible to find a positive ${\mu }_1$ that is sufficiently small to make $s_1=202.628$ .

Wang and Zhou (2018) mention two mechanisms through which workload could impact server performance for the single-server queues in their setting: social speedup pressure and customer early task initiation. Hypothetically, if their mean service times had been as assumed in this subsection, then the failure to translate the mean service times into a dynamic model could be taken as evidence against social speedup pressure—a dynamic mechanism. The large mean service time for customers who arrive at an idle server ( ${\hat{s}}_1$ ) could be seen as evidence of a physical setup mechanism (a static mechanism), perhaps because cashiers need to log into the cash register before serving the first customer after a break.

6.3. Illustrations of Faulty Generalizations (1) and (2)

Next, we present two examples to illustrate faulty generalizations (1) and (2), for the translation of mean service times into dynamic models. For each example, similar to Case Study 2, we specify the true model, simulate data from the true model, estimate mean services, and use (6) to translate the estimated mean service times into service rates.

The true models are $M/M_n/4$ models with $\lambda =2(\lambda =3)$ for Example 1 (Example 2). We use the average of two sigmoid functions (see EC.4.3 of E-Companion for details) to generate ${\mu }_n$ . Weighted sums of sigmoid functions provide a flexible way to generate monotone and non-monotone service rates with finite limits and have been validated by real-world data (Cho et al., 2019).

6.3.1. Example 1 True Model Specification and Simulation

We use a U-shaped ${\mu }_n$ , shown in Figure 6(b). We simulate the true model for 100,000 s and exclude customers who arrived during the first 20,000 s as a warm-up period. We obtain service times for 160,345 customers and the largest starting state $n$ for service is 17.

6.3.2. Mean Service Time Estimation

We begin by estimating ${\hat{s}}_n$ as the average of all service times that start in state $n$ . The resulting estimates are noisy for rarely visited states (see Figure 6(a)). To reduce noise, we extrapolate ${\hat{s}}_n$ smoothly by setting ${\hat{s}}_{n+1}={\hat{s}}_n+({\hat{s}}_n-{\hat{s}}_{n-2})/2$ for $n=13$ to $17$ . Furthermore, to translate the mean service times into service rates, we must extend ${\hat{s}}_n$ beyond the range of the sample, and we do that by setting ${\hat{s}}_n={\hat{s}}_{17}$ for $n\ge 18$ . Figure 6(a) shows the resulting mean service times as ${\hat{s}}_n$ (extended).

6.3.3. Translation

Figure 6(b) shows the ${\hat{\mu }}_n$ (extended) obtained from ${\hat{s}}_n$ (extended) using (6). We also plot $1/{\hat{s}}_n$ (extended) from faulty generalization (1) and the actual service rates ${\mu }_n$ for comparison. The ${\hat{\mu }}_n$ (extended) matches the true ${\mu }_n$ pattern better with a mean absolute percentage error (MAPE) of $6.3\mathrm{\%}$ . In contrast, even though $1/{\hat{s}}_n$ (extended) captures the directional change in ${\mu }_n$ correctly, it results in a higher MAPE of $8.4\mathrm{\%}$ .

6.3.4. Example 2 True Model Specification and Simulation

We use an inverted-U-shaped ${\mu }_n$ , shown in Figure 7(b). We simulate the true model for $200,000/3$  s and exclude customers who arrived during the first $40,000/3$  s as a warm-up period. We obtain service times for $159,905$ customers and the largest starting state $n$ for service is $12$ .

6.3.5. Mean Service Time Estimation

We use the same approach as in Example 1 to obtain ${\hat{s}}_n$ ; for ${\hat{s}}_n$ (extended), we extrapolate ${\hat{s}}_n$ smoothly for $n=10$ to $12$ , and set ${\hat{s}}_n$ $\text{(extended)}={\hat{s}}_{12}$ (extended) for $n\ge 13$ . Figure 7(a) shows the resulting ${\hat{s}}_n$ and ${\hat{s}}_n$ (extended).

Figure 7.

Faulty generalization (2) example: (a) comparison of service times and (b) comparison of service rates.

6.3.6. Translation

Figure 7(b) shows that ${\hat{\mu }}_n$ (extended) successfully captures the increase–decrease ${\mu }_n$ pattern for States 1 to 7 with a MAPE of $9.4\mathrm{\%}$ . In contrast, $1/{\hat{s}}_n$ (extended) incorrectly shows an increasing pattern for State 7 with a higher MAPE of $16.5\mathrm{\%}$ , which supports our concern regarding faulty generalization (2).

As a check on our work for both examples, we compute the true $s_n$ using (EC.34 of E-Companion), and find that their values agree closely with ${\hat{s}}_n$ (extended) for frequently visited states, but diverge slightly for $n\ge 12$ ( $n\ge 8$ ) for Example 1 (Example 2).

6.4. Convex Decreasing Mean Service Times Are Not Equivalent to Increasing Service Rates

Convex decreasing $s_n$ guarantees increasing ${\mu }_n$ for single-server dynamic models (Corollary 4) but the same is not true for multi-server dynamic models. We use two examples to illustrate this.

Example 1

We set $\lambda =2$ , $c=3$ , with $s_n$ convex decreasing ( $\mathrm{\Delta }{s}_n\le 0$ ; ${\mathrm{\Delta }}^2{s}_n\ge 0$ ) (see Figure 8(a)). The resulting ${\mu }_n$ is non-monotone, demonstrating that convex decreasing $s_n$ can result in non-increasing ${\mu }_n$ .

Figure 8.

Numerical examples for Theorem 2: (a) convex decreasing $s_n$ ; non-increasing ${\mu }_n$ and (b) non-convex $s_n$ ; increasing ${\mu }_n$ .

Example 2

We set $\lambda =1.2$ , $c=4$ , with non-convex $s_n$ obtained from a sigmoid function (see Figure 8(b); see EC.4.3 of E-Companion for details). The resulting ${\mu }_n$ is increasing, demonstrating that non-convex $s_n$ can result in increasing ${\mu }_n$ . We verified the results for both examples using simulation.

6.5. Measurement Timing can Reverse the Direction of Mean Service Times

Our analysis in Sections 4 and 5 demonstrates that in a single-server dynamic model, mean service time patterns differ markedly depending on whether the workload is measured at the beginning of service ( $s_n$ in (10)) or at the end of service ( $s_n^l$ in (26)).

To illustrate, we set $\lambda =3$ , $c=1$ , and use a sigmoid function to obtain ${\mu }_n$ (see Figure 9; see EC.4.3 of E-Companion for details). We compute $s_n$ and $s_n^l$ using (10) and (26), respectively (see Figure 9). The server utilization is $\rho =97.5\mathrm{\%}$ . The ${\mu }_n$ increases, which results in decreasing $s_n$ , but increasing $s_n^l$ . We verified these results using simulation.

Figure 9.

Measurement timing example: comparison of $s_n$ and $s_n^l$ curves.

7. Conclusion

A growing body of empirical work documents the adaptive behavior of customers and servers in service operations, including a research stream that focuses on how workload impacts service times. We provide a framework and draw attention to important subtleties in the process of translating empirical findings into the selection and specification of state-dependent queueing models. Among the available models, we identify two model types (static and dynamic) and relate them to behavioral mechanisms. For dynamic models, we provide a formula to compute service rates from mean service times and caution against a pitfall in assessing server behavior (see (2)). This translation process is crucial because it connects empirical research to queueing models that are employed in the management of service and production operations. In particular, empirically grounded dynamic models have been used to support operational decisions such as system design and staffing (Dong et al., 2015; Do et al., 2018; Cho et al., 2019).

The static model that we focus on has not been analyzed exactly, and we focus on its insights into the service process from a behavioral perspective. Mathematical analysis of static models is challenging, but simulation analysis of such models is straightforward. Static and dynamic models correspond to different behavioral mechanisms. The choice of model type should depend on both the hypothesized mechanisms and the consistency with data.

We extend the analysis of dynamic models from the system level to the individual customer level. By formulating service times as first-passage times, we show that service times in dynamic models follow PH distributions. For single-server dynamic models, we show that workload measurement timing can have a drastic impact on service time patterns.

We identify two potential mistakes in translating service times into service rates, (1) and (2), for dynamic models. Using first-step analysis, we provide a closed-form solution for ${\mu }_n$ from $s_n$ , which corrects (1). We provide a series of conditions in the form of bounds on the forward differences of mean service times that guarantee (2) to hold. Our results demonstrate what mean service times imply about how service rates vary with workload, and reveal speed-ups and slow-downs hidden in mean service times. These results enable researchers to correctly assess the impact of different mechanisms on server performance.

Building on published empirical findings, we demonstrate two possible scenarios. First, it is possible that models of both types can be specified from the same mean service times. In this scenario, further investigation of system performance measures is needed to determine the more suitable model type. Second, the outcome of the translation process could indicate the inconsistency of dynamic models with mean service times, which provides evidence against this model type and thus its associated mechanisms. Despite the greater tractability of dynamic models, which makes it easier to obtain analytical results regarding, for example, admission policies, we demonstrate fitting them to data is not always possible.

Our work suggests several topics that would benefit from further investigation. First, generalizing static and dynamic models to non-Markovian settings would benefit researchers by providing them with a less restricted toolkit. In particular, analysis of dynamic models with non-Markovian sojourn times could be a fruitful direction for future research. Comparisons of such models with static models with general service times, as well as the translation of mean service times into the inputs needed for such models, could provide insights for both the analytical and empirical literature streams.

Second, the translation process from empirical findings to a mathematical model may require extrapolation beyond a finite dataset. Most queueing models assume unlimited waiting room, which means that one needs to specify how service time distributions or rates vary as the queue length goes to infinity. In Section 6, we address this issue by extrapolating ${\hat{s}}_n$ smoothly and assuming constant $s_n$ for unobserved states. However, this is only one possible approach: if a log-linear regression is used to estimate mean service times, then one natural approach is to extrapolate the curve beyond the data. This approach can be problematic because it typically results in either negative or very large mean service times as $n\to \mathrm{\infty }$ .

Third, not only should the hypothesized mechanisms inform the choice of queueing model type; but they should also inform the choice of statistical methods. For example, if the researchers suspect that server speed changes dynamically during service, survival analysis methods could be appropriate (see Ding et al., 2024 for an example).

Fourth, our findings in Sections 5 and 6.5 indicate that measurement timing could be an important factor in determining service time patterns. As support, Table 1 lists five papers that measure workload only at the beginning of service, of which three report $\nearrow \searrow$ service time patterns. One paper measures workload at the end of service and finds a $\searrow$ service time pattern. From an empirical perspective, definitions of mechanisms (such as those proposed in Table 2 of Delasay et al., 2019) should include measurement timing, which could help researchers in data collection and model selection. The names of the server early task initiation and customer early task initiation mechanisms imply that they involve a reduction in work content that occurs near the beginning of service. In contrast, the definition in Table 2 of Delasay et al. (2019) of the task reduction mechanism is silent regarding the timing of work content reduction. Kc and Terwiesch (2012) argue convincingly that task reduction occurs near the end of service in ICU settings. However, in other settings, task reduction could in principle occur earlier during the service. From an analytical perspective, a “late task reduction” mechanism is inconsistent with the static and dynamic model types that we focus on in this paper, but it might be consistent with the type of dynamic model proposed in Chan et al. (2012).

Supplemental Material

sj-pdf-1-pao-10.1177_10591478241309662 - Supplemental material for Translating Empirical State-Dependent Service Times Into Queueing Models

Supplemental material, sj-pdf-1-pao-10.1177_10591478241309662 for Translating Empirical State-Dependent Service Times Into Queueing Models by Likang Ding, Bora Kolfal and Armann Ingolfsson in Production and Operations Management