Sensitivity Analysis and Simulation of a Multiserver Queueing System with Mixed Service Time Distribution

Abstract

The motivation of mixing distributions in communication/queueing systems modeling is that some input data (e.g., service time in queueing models) may follow several distinct distributions in a single input flow. In this paper, we study the sensitivity of performance measures on proximity of the service time distributions of a multiserver system model with two-component Pareto mixture distribution of service times. The theoretical results are illustrated by numerical simulation of the $M/G/c$ systems while using the perfect sampling approach.

1. Introduction

Mixtures of distributions arise in complex stochastic systems and they are extensively used for statistical analysis in many real fields, such as lifetime modeling, ageing or failure processes, engineering reliability [1], and survival theory [2], where data are assumed to be heterogeneous. The application of the mixture of distributions in the modeling of queueing systems is often induced by diverse structure of the customers in the system, e.g., by various service time requirements of multiple classes of customers that arrive into the system (for instance, the transmission time of IP datagrams with different lengths), or by the noisy/biased measurements that induce the so-called contaminated distributions. Ignoring such a diversity at the modeling phase may lead to significant deviation of system performance at practical implementation phase as compared to the modelled values. This motivates various types of analysis, including the analysis of continuity, robustness, monotonicity, stability, and sensitivity. In this regard, we mention the fundamental result obtained for telecommunication system models by B. A. Sevast’yanov [3], and the basic monographs [4,5].

The authors would like to use this opportunity to pay tribute to Professor Vladimir Zolotarev and to note his outstanding role as the founder of the International Seminar on Stability Problems for Stochastic Models. One of the authors had a great pleasure to communicate with Professor Zolotarev over many years, and all of the authors actively participated in the seminar he has founded. In the context of this paper, it is especially appropriate to emphasize an important role of Professor Zolotarev in the study of the stability and monotonicty of queuieng processes, see [6,7,8].

Information flows in modern telecommunications and computing systems have the form of a superposition of some sequential-parallel structures [9]. Ranging from small personal devices up to large scale high-performance computing systems, all of these may be modeled as multiserver queueing systems. Thus, it is highly important to study the performance of such systems and, in particular, the sensitivity of stationary performance indexes with respect to the variability of input parameters. However, direct output analysis of queueing systems is often tricky (see e.g., [10,11]), and explicit expressions for the distributions of steady-state performance indexes of a multiserver system are, in general, hardly available and, beyond classic models, known in some special cases only. In some cases, the analysis may be performed by obtaining asymptotic upper bounds, as in the paper [12], or studying the continuity of the process, as in [8], or stochastic stability of the queueing process, like in [4,13], or by means of simulation. In the present paper, we utilize the latter approach.

This paper is dedicated to sensitivity analysis of a steady-state performance index of a multiserver system with respect to service time distribution having the form of the so-called finite mixture [14]. However, instead of studying the direct parametric sensitivity, we focus on a more delicate analysis of the (combined) effect of the service time distribution on the steady-state performance estimate. That is, we compare the basic system to a disturbed one, using a sensitivity measure (Kolmogorov–Smirnov distance) both for the service time distributions, and for the steady-state performance estimate (queue size). The service time distribution perturbation is performed by changing the mixing coefficient and parameters defining the mixture components. We formalize this at the end of Section 3.

In general, the output distributions are hardly analytically available and, in this case, we must be able to obtain the steady-state performance indexes by simulation. As a basic model, we consider the classical $M/G/c$ model, where the steady-state distribution of the vector workload process is unknown as well; however, it can be estimated by means of the recently developed method of regenerative perfect simulation [15]. In more detail, as the target (perturbed) service time distribution we take the two-component mixture of Type-II Pareto distributions with support on the positive axis, which is known as Lomax distributions, as well as two-component exponential (hyperexponential) distribution. Such a choice also allows for obtaining some analytical expressions. Our interest to Pareto distribution is caused by the heavy tailed property of this distribution that is frequently observed in models of file size and flow duration [16].

This paper continues the study performed in [17] in the context of monotonicity. The key idea of the present paper is to study qualitatively the sensitivity of the steady-state distribution of the system performance index (steady-state queue size) to the variability of service time distribution by means of simulation. We also apply the auxiliary results on the failure rates comparison, which allows us to characterize the monotonicity of some stationary performance measures.

The structure of the paper is as follows. In Section 2, we introduce the two-component mixture of distributions and discuss some properties that are used in the subsequent analysis. Subsequently, we define the uniform distance between the mixture and the corresponding parent distribution. In Section 3, some known stochastic monotonicity properties of the multiserver system are collected, which further are specified for the considered mixture distributions. In Section 4, we describe the perfect sampling algorithm that is then used to sample from exact (but unknown) steady-state distribution of a multiserver queue $M/G/c$. The results of simulation are presented in Section 5. We study the sensitivity of the steady-state queue size distribution with respect to (w.r.t.) the shape parameters of mixture and the mixing coefficient and illustrate stochastic monotonocity of the system performance. The discussion of the simulation results finalizes the paper in the concluding Section 6.

2. Two-Component Mixture Distributions

The goal of this Section is to derive the uniform distance between the two-component mixture distribution and its parent distribution. First, we introduce the two-component mixture, and then give a few properties, including the stochastic monotonicity. This property is further used to obtain the monotonicity of the corresponding output queueing process.

Let $X_i$ be independent random variables having mean $\mathrm{E}{X}_i$, density $f_i$, tail distribution function (d.f.) ${\overline{F}}_i\left(x\right)=1-F_i\left(x\right)$, and failure rate

$$r_i\left(x\right)=\frac{f_i\left(x\right)}{{\overline{F}}_i\left(x\right)},i=1,2,$$

defined for such x that ${\overline{F}}_i\left(x\right)>0$. We assume that $F_1\ne F_2$ to avoid trivial case. Let I be a Bernoulli random variable independent of $X_i$, with success probability $P(I=1)=p$. Subsequently, it is called the random variable

$$X_M=I{X}_1+(1-I)X_2,$$

has the two-component mixture distribution [18] (we use the index M to denote the mixture). The mean $\mathrm{E}{X}_M$ and density, $f_M$ of $X_M$ equal, respectively,

$$\begin{array}{cc}\hfill \mathrm{E}{X}_M& =p\mathrm{E}{X}_1+(1-p)\mathrm{E}{X}_2,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill f_M\left(x\right)& =p{f}_1\left(x\right)+(1-p)f_2\left(x\right),\hfill \end{array}$$

(2)

and it is easy to see that the tail distribution is

$${\overline{F}}_M\left(x\right)=p{\overline{F}}_1\left(x\right)+(1-p){\overline{F}}_2\left(x\right).$$

(3)

Note that the d.f. $F_i$ may belong to the same family of distributions but have other parameters. In reliability analysis, such a mixture may be interpreted as a contaminated distribution [19], where $1-p$ is, as a rule, small enough. $F_1$ is called the parent distribution and $F_2$ is the contaminating distribution. In this Section, we focus on the distance between the mixture and its parent distribution.

A straightforward analysis shows that the failure rate of the mixture has the following form [20]:

$$r_M\left(x\right)=\frac{p{f}_1\left(x\right)+(1-p)f_2\left(x\right)}{p{\overline{F}}_1\left(x\right)+(1-p){\overline{F}}_2\left(x\right)}=a\left(x\right)r_1\left(x\right)+(1-a\left(x\right))r_2\left(x\right),$$

(4)

where

$$a\left(x\right)=\frac{p{\overline{F}}_1\left(x\right)}{p{\overline{F}}_1\left(x\right)+(1-p){\overline{F}}_2\left(x\right)},x\ge 0.$$

In particular, it follows from Equation (4) that

$$r_M\left(x\right)\ge min(r_1\left(x\right),r_2\left(x\right)),x\ge 0.$$

(5)

It is worth mentioning that the mixture preserves the monotonicity of failure rate in the following way: if both rates $r_i\left(x\right)$ are non-increasing, that is d.f.’s $F_i\left(x\right)$ are decreasing failure rate distributions (DFR), then the mixture $F_M\left(x\right)$ is DFR distribution as well [21]. Indeed, one can check that

$$r_M^{\prime }\left(x\right)=a\left(x\right)r_1^{\prime }\left(x\right)+(1-a\left(x\right))r_2^{\prime }\left(x\right)-a\left(x\right)(1-a\left(x\right)){(r_1\left(x\right)-r_2\left(x\right))}^2,$$

(6)

also see [1]. Subsequently, if $r_i^{\prime }\left(x\right)<0,i=1,2$, it follows from Equation (6) that $r_M^{\prime }\left(x\right)<0$, since $a\left(x\right)\in [0,1]$ for any $x\in (0,\infty )$. In particular, it follows from Equation (6) that the mixture of two exponential distributions is DFR (note that the exponential distribution has constant failure rate).

Another example of a DFR distribution is the Type-II Pareto distribution, denoted by $Pareto({\alpha }_i,x_0)$, having d.f. (see e.g., [22])

$$F_i\left(x\right)=1-{\left(\frac{x_0}{x_0+x}\right)}^{{\alpha }_i},x\ge 0,x_0>0,{\alpha }_i>0,i=1,2.$$

The failure rate of $Pareto({\alpha }_i,x_0)$ equals

$$r_i\left(x\right)=\frac{{\alpha }_i}{x_0+x},x\ge 0,i=1,2,$$

(7)

and it is monotonically decreases to 0 as $x\to \infty$. As has been noted above, the two-component mixture of Pareto distributions is DFR distribution. However, the failure rate of finite mixtures, in general, is a complicated function [20].

The uniform distance between distributions F and G, defined as [12]

$$\mathsf{\Delta }(F,G)=\underset{x}{sup}|F\left(x\right)-G\left(x\right)|,$$

(8)

is a recognised measure, which is actively used in the sensitivity analysis [12]. It is easy to see that the uniform distance $\mathsf{\Delta }(F_M,F_1)$ between the mixture distribution Equation (3) and its parent distribution is

$$\mathsf{\Delta }(F_M,F_1)=\underset{x\ge 0}{sup}|p{F}_1\left(x\right)+(1-p)F_2\left(x\right)-F_1\left(x\right)|=(1-p)\underset{x\ge 0}{sup}|F_1\left(x\right)-F_2\left(x\right)|.$$

(9)

Note that, if the densities $f_i$ exist, and there exists $x^{*}$ that delivers the supremum in Equation (9),

$$\mathsf{\Delta }(F_M,F_1)=|F_1\left(x^{*}\right)-F_2\left(x^{*}\right)|,$$

then $x^{*}$ satisfies the equality

$$f_1\left(x^{*}\right)=f_2\left(x^{*}\right).$$

(10)

By definition of the failure rates, $r_i$, it then follows that

$$r_1\left(x^{*}\right){\overline{F}}_1\left(x^{*}\right)=r_2\left(x^{*}\right){\overline{F}}_2\left(x^{*}\right).$$

Thus, expression Equation (9) can be written in the following convenient form

$$\mathsf{\Delta }(F_M,F_1)=(1-p)\frac{|r_2\left(x^{*}\right)-r_1\left(x^{*}\right)|}{r_2\left(x^{*}\right)}{\overline{F}}_1\left(x^{*}\right)=(1-p)\frac{|r_1\left(x^{*}\right)-r_2\left(x^{*}\right)|}{r_1\left(x^{*}\right)}{\overline{F}}_2\left(x^{*}\right).$$

(11)

Note that Equation (11) allows obtaining the following upper bound for the distance $\mathsf{\Delta }(F_M,F_1)$:

$$\mathsf{\Delta }(F_M,F_1)\le (1-p)\frac{|r_2\left(x^{*}\right)-r_1\left(x^{*}\right)|}{r_1\left(x^{*}\right)}=:\delta \left(x^{*}\right).$$

(12)

In particular, for the hyperexponential distribution, that is for two-component mixture of exponential distributions with densities $f_i\left(x\right)={\lambda }_i{e}^{-{\lambda }_{i}x},i=1,2$, it follows from Equation (10), that

$$x^{*}=\frac{log{\lambda }_1-log{\lambda }_2}{{\lambda }_1-{\lambda }_2},$$

and in this case expression Equation (11) becomes

$$\mathsf{\Delta }(F_M,F_1)=(1-p)\frac{|{\lambda }_2-{\lambda }_1|}{{\lambda }_2}{\left(\frac{{\lambda }_1}{{\lambda }_2}\right)}^{-\frac{{\lambda }_1}{{\lambda }_1-{\lambda }_2}}\le (1-p)\frac{|{\lambda }_2-{\lambda }_1|}{{\lambda }_2}.$$

(13)

Note that the last inequality in Equation (13) is a particular case of Equation (12). Expression Equation (13) is consistent with a more general result for the so-called univariate scale mixture $X_M$ having form [2]

$$X_M\stackrel{d}{=}\frac{X_1}{Y},$$

(14)

with d.f.

$${\widehat{F}}_M\left(x\right)={\int }_0^{\infty }{F}_1\left(\theta x\right)dG\left(\theta \right),$$

where $F_1$ is the parent distribution of the random variable $X_1$ and G is the distribution of a mixing random variable $Y\ge 0$. It is clear that the transformation Equation (14) is a scale change, and if $Y\in \{y_1,\dots ,y_m\}$ is a discrete random variable, then Equation (14) becomes

$$X_M=\sum _{i=1}^m\frac{1}{y_i}I(Y=y_i)X_1,$$

(15)

which is a finite mixture, where I is an indicator function. The aforementioned general result for the univariate scale mixture states that if $F_1$ is an exponential d.f., then an upper bound for the uniform distance may be obtained as follows [2,23]:

$$\mathsf{\Delta }({\widehat{F}}_M,F_1)\le \mathrm{E}|Y-1|.$$

(16)

To show that Equation (16) indeed coincides with Equation (13) for the two-component scale mixture case, let Y have point masses at 1 and ${\lambda }_2/{\lambda }_1$ with probabilities p and $1-p$, respectively. It immediately follows from Equations (15) and (16) that

$$\mathrm{E}(Y-1)=(1-p)\frac{|{\lambda }_1-{\lambda }_2|}{{\lambda }_2}.$$

Now, we return to the two-component $Pareto({\alpha }_i,x_0)$ mixture $F_M$. It follows from Equation (11) that in this case

$$\mathsf{\Delta }(F_M,F_1)=(1-p)\frac{|{\alpha }_1-{\alpha }_2|}{{\alpha }_2}{\left(\frac{{\alpha }_2}{{\alpha }_1}\right)}^{\frac{{\alpha }_1}{{\alpha }_1-{\alpha }_2}},$$

(17)

where the value $x^{*}$ satisfying equality Equation (10) equals

$$x^{*}=x_0\left({\left(\frac{{\alpha }_1}{{\alpha }_2}\right)}^{\frac{1}{{\alpha }_1-{\alpha }_2}}-1\right).$$

Note that the r.h.s. of Equation (17) is similar to the r.h.s. of Equation (13). This similarity is caused by the specific shape of the failure rate of the distribution $Pareto({\alpha }_i,x_0)$. Moreover, in such a case, the quantity $\delta \left(x^{*}\right)$ defined in Equation (12) does not depend on $x^{*}$ and, thus, for Pareto mixture, it readily follows from Equation (12) that

$$\mathsf{\Delta }(F_M,F_1)\le (1-p)\frac{|{\alpha }_1-{\alpha }_2|}{{\alpha }_2}=:\delta ({\alpha }_1,{\alpha }_2).$$

(18)

In Figure 1, to illustrate the dependence of the uniform distance on the parameter ${\alpha }_2$ of the contaminating distribution, we depict $\mathsf{\Delta }(F_M,F_1)$ jointly with $\delta ({\alpha }_1,{\alpha }_2)$ for fixed ${\alpha }_1=2$ and $p=0.9$ by varying ${\alpha }_2$ in the interval $(1,5)$.

3. Multiserver System Sensitivity

In this Section, we formalize our main goal for the numerical experiments conducted and discussed in Section 5. We then demonstrate how stochastic and failure rate ordering can be applied to multiserver systems with mixed service time distribution. The numerical experiments equipped with the stochastic comparison technique not only allow for obtaining the absolute value, but also characterizing the monotonicity of performance indexes.

Consider a classical First-Come-First-Served (FCFS) c-server $M/G/c$ queueing system that is fed by a Poisson input with rate $\lambda$, arrival instants $\{t_i,i\ge 1\}$ with $t_1=0$, independent and identically distributed (iid) interarrival times $T_i=t_{i+1}-t_i$ and iid service times $\{S_i,i\ge 1\}$. Note that $\lambda =1/\mathrm{E}T$, where T is generic interarrival time. Now, we consider the c-dimensional vector of the remaining workload process in such a system,

$$W_i=(W_{i,1},\dots ,W_{i,c}),$$

where $W_{i,k}$ is the kth smallest component of the vector which is observed by the ith arrival [24]. Thus, the vector components are kept in ascending order,

$$W_{i,1}\le \cdots \le W_{i,c},$$

and the quantity $W_{i,j}$, “observed” by the arriving customer i, equals the unfinished work which must be done by server j provided no new work arrives after arrival instant $t_i$ of customer $i;j=1,\dots ,c$. If there are no idle servers upon arrival of customer i, then s/he waits in a common infinite capacity queue until the server with minimal work, $W_{i,1}$, becomes free. It is easy to see that $W_{i,1}$ is the waiting time of customer i which starts being served at time $t_i+W_{i,1}$. It is well-known that the workload vector sequence follows the celebrated stochastic Kiefer–Wolfowitz recursion [25]:

$$W_{i+1}=R{(W_i+e_1{S}_i-\mathbf{1}{T}_i)}^{+},$$

(19)

where $e_1=(1,0,\dots ,0)$ and $\mathbf{1}$ is the vector of ones, operator R puts the components in an ascending order, and operation ${(·)}^{+}=max(0,·)$ is applied componentwise (we omit the sub-index for a generic element of a sequence). In what follows, we assume that the stability condition holds [25],

$$\begin{array}{c}\hfill \rho :=\lambda \mathrm{E}S<c.\end{array}$$

(20)

Define the departure instant of customer i by $d_i=t_i+W_{i,1}+S_i$. Now define the process

$$Q_n=\sum _{j\ge 1,j\ne n}I(t_j\le t_n<d_j),$$

(21)

counting the queue size (number of customers in the system) at the arrival instant $t_n$. Under condition Equation (20), $Q_n$ converges in distribution, as $n\to \infty$, to the steady-state queue size Q, with stationary distribution

$${\pi }_n=\mathrm{P}(Q=n),n\ge 0.$$

Note that when service times $S_i$ are exponential, the steady-state queue size distribution, ${\pi }_n,n\ge 0$, is well known [26]:

$${\pi }_n=\left\{\begin{array}{cc}{\left({\displaystyle \sum _{k=0}^{c-1}}\frac{{\displaystyle {\rho }^k}}{{\displaystyle k!}}+\frac{{\displaystyle {\rho }^c}}{{\displaystyle (c-1)!(c-\rho )}}\right)}^{-1},\hfill & n=0,\hfill \\ {\pi }_0\frac{{\displaystyle {\rho }^n}}{{\displaystyle n!}},\hfill & 1\le n\le c,\hfill \\ {\pi }_0\frac{{\displaystyle {\rho }^n}}{{\displaystyle c!c^{n-c}}},\hfill & n>c.\hfill \end{array}\right.$$

(22)

The operators $R(·)$ and ${(·)}^{+}$ in Equation (19) preserve ordering, and it allows for us to establish the monotonicity of the workload sequence in the multiserver system in the case when the driving sequences $\{T_n^{\left(i\right)},S_n^{\left(i\right)},n\ge 1\}$, $i=1,2$ satisfy stochastic order. We recall that the stochastic order $X_2{\le }_{st}{X}_1$ between two random variables $X_1,X_2$ means that the tail d.f.’s satisfy inequality

$$\mathrm{P}(X_2>x)\le \mathrm{P}(X_1>x),x\ge 0.$$

(23)

In is known [27] that, in two c-server systems with stochastically ordered input sequences, $T^{\left(2\right)}{\ge }_{st}{T}^{\left(1\right)}$ and $S^{\left(2\right)}{\le }_{st}{S}^{\left(1\right)}$, the workload sequences $\left\{W_n^{\left(i\right)}\right\},i=1,2$, are (componentwise) ordered in the following way

$$W_n^{\left(2\right)}{\le }_{st}{W}_n^{\left(1\right)}n\ge 1.$$

(24)

It also holds for the steady-state workloads:

$$W^{\left(2\right)}{\le }_{st}{W}^{\left(1\right)}.$$

If the input in both systems is the same, which is $T^{\left(2\right)}{=}_{st}{T}^{\left(1\right)}$, then the the queue length process at the arrival instants satisfy similar ordering both in path-wise sense and in steady-state [27]

$$Q^{\left(2\right)}{\le }_{st}{Q}^{\left(1\right)}.$$

(25)

The stochastic ordering ${\le }_{st}$ can be transformed into the ordering with probability 1 by the coupling technique [28]. In the context of this work, it is worth mentioning that the sufficient condition for the stochastic ordering $S^{\left(2\right)}{\le }_{st}{S}^{\left(1\right)}$ is the failure rate ordering [29]:

$$r_2\left(x\right)\ge r_1\left(x\right),x\ge 0,$$

(26)

where $r_i$ is the failure rate of r.v. $S^{\left(i\right)}$, $i=1,2$. We summarize the discussion in the following lemma which is a straightforward result of [27].

Lemma 1.

Consider two c-server systems with stochastically equivalent input, $T^{\left(2\right)}{=}_{st}{T}^{\left(1\right)}$, and failure rate ordered service time distributions, $r_2\left(x\right)\ge r_1\left(x\right),x\ge 0$. Subsequently, Equation (25) holds.

Now, we consider two $M/G/c$ queueing systems, denoted by ${\mathsf{\Sigma }}^{\left(1\right)}$ and ${\mathsf{\Sigma }}^{\left(M\right)}$, fed by (stochastically) identical Poisson process with rate $\lambda$. Let the first system ${\mathsf{\Sigma }}^{\left(1\right)}$ have the service time distribution $F_1$. We refer below to the first system as being basic. In the second (contaminated) system ${\mathsf{\Sigma }}^{\left(M\right)}$, we use service time distribution $F_M$ defined by Equation (3), with the same $F_1$ and some $F_2$ and $p\in (0,1)$. Let now $Q^{\left(1\right)}$ (with d.f. $F_{Q^{\left(1\right)}}$) be the steady-state queue size in the first system. Define similarly $Q^{\left(M\right)}$ and $F_{Q^{\left(M\right)}}$ for the system ${\mathsf{\Sigma }}^{\left(M\right)}$. We are interested in studying the sensitivity of the uniform distance

$$\mathsf{\Delta }(F_{Q^{\left(M\right)}},F_{Q^{\left(1\right)}})=\underset{x\ge 0}{sup}|F_{Q^{\left(1\right)}}\left(x\right)-F_{Q^{\left(M\right)}}\left(x\right)|.$$

(27)

More formally, we study the effect of $\mathsf{\Delta }(F_M,F_1)$ given by Equation (9), on the steady-state performance $\mathsf{\Delta }(F_{Q^{\left(M\right)}},F_{Q^{\left(1\right)}})$ defined in Equation (27), by varying the mixing coefficient p and parameters defining the mixture components $F_1$ and $F_2$. However, since the distributions $F_{Q^{\left(M\right)}}$ and $F_{Q^{\left(1\right)}}$ are not available explicitly in general, we use simulation to obtain the corresponding estimates. As such, we study a combined effect of the service time distribution on the steady-state performance estimate.

The generic service time $S^{\left(M\right)}$ in the contaminated system ${\mathsf{\Sigma }}^{\left(M\right)}$ has a two-component mixture d.f. $F_M$ and, thus, it follows from Equation (5) that the conditions of Lemma 1 are satisfied, since M, where $r_M$ is the failure rate of $S^{\left(M\right)}$. In particular, this means that the basic system ${\mathsf{\Sigma }}^{\left(1\right)}$ is heavier loaded than the contaminated system ${\mathsf{\Sigma }}^{\left(M\right)}$. It then follows from Lemma 1 and Equation (23) that the difference $F_{Q^{\left(1\right)}}\left(x\right)-F_{Q^{\left(M\right)}}\left(x\right)$ (see Equation (27)) is negative for all $x\ge 0$. In Section 5, we study the distance Equation (27) numerically.

4. Exact Steady-State Simulation by Regenerative Approach

In general, there are no closed form expressions for the steady-state distribution of the queue length and vector workload process in an $M/G/c$ system. Although a number of approximations exist [30,31,32,33], in general the accuracy of such methods is a point of discussion [34], especially when the service times distribution is heavy-tailed. Thus, to study the sensitivity we need to rely on simulation. A contribution of this work is that unlike classical discrete-event simulation (crude Monte-Carlo), which always has the so-called transient (warm-up) period during which an influence of initial conditions exists, we use the perfect simulation technique that allows exact sampling from the (unknown) steady-state distribution. In what follows, we rely on the regenerative approach designed for the $M/G/c$ system in the work [15] (although there are recently developed more sophisticated techniques based on backward coupling, for instant [35], which are valid for a more general $G/G/c$ system). Below, we outline the approach from [15].

This approach uses the so-called a Random Assignment (RA) system $M/G/c$ as a majorant for the original $M/G/c$ system. In the RA system, each new customer is assigned to arbitrary server randomly (that is with probability $1/c$). As a result, the remaining workload in server j that customer n meets, denoted by $V_{n,j}$, satisfies recursion

$$V_{n+1,j}={[V_{n,j}+I(U_n=j)S_n-T_n]}^{+},j=1,\dots ,c,$$

(28)

where iid random variables $\left\{U_n\right\}$ are uniformly distributed over $\{1,\dots ,c\}$, and $I(U_n=j)=1$ means that customer n is routed to server j. The RA system is indeed is a collection of $M/G/1$ systems, each with Poisson input with rate $\lambda /c$. As a result, in each, such a system the stationary workload, D, is distributed in accordance with the following version of the Pollaczek–Khintchine formula [15]

$$D=\sum _{i=1}^L{S}_i^{\left(e\right)},$$

(29)

where L has geometric distribution

$$P(L=k)={\left(\frac{\rho }{c}\right)}^k\left(1-\frac{\rho }{c}\right),$$

(30)

and $S^{\left(e\right)}$ has the so-called equilibrium (integrated tail) distribution,

$$\mathrm{P}\left(S^{\left(e\right)}>x\right)=\frac{1}{\mathrm{E}S}{\int }_x^{\infty }{\overline{F}}_S\left(t\right)dt,$$

(31)

where $F_S$ is the d.f. of original service time S. It is well-known that both workload process and queue size process in the RA system dominate the corresponding process in the original $M/G/c$ system [5,24,27]. Applying coupling, this dominance holds with probability (w.p.) 1. In particular, the regenerations of RA system (the instants when customers meet totally idle system) are also regeneration instants of the original system $M/G/c$. These results then are used to sample from the steady-state distribution of the RA system as follows:

• sample the values $L_i,i=1,\dots ,c$ according to geometric distribution Equation (30);
• sample $S_1^{\left(e\right)},\dots ,S_{L_i}^{\left(e\right)},i=1,\dots ,c$ according to integrated tail distribution Equation (31); and,
• construct the (stationary) components $D_i$ for $i=1,\dots ,c$, by formula Equation (29).

Subsequently, starting from the steady-state vector $V_1=(D_1,\dots ,D_c)$ containing iid components, the recursion Equation (28) is applied to each separate queue in the RA system until the event

$$V_{{\tau }_e}=(V_{{\tau }_e,1},\dots ,V_{{\tau }_e,c})=\mathbf{0}$$

happens at the (arrival) instant of some customer ${\tau }_e$. Thus, ${\tau }_e$ is the length of equilibrium (steady-state) remaining regeneration period. Note that by construction, at each step of this recursion, the workload vector has steady-state distribution in the RA system. Omitting unnecessary details, the remaining steps of algorithm are as follows [15]:

• sample stochastic copies $V^{\left(k\right)}=(V_1^{\left(k\right)},V_2^{\left(k\right)},\dots ),k=1,2,\dots$ of the sequence of workload vectors using recursion Equation (28); each sequence starts with $V_1^{\left(k\right)}=0$ and lasts until the event $V_{\tau \left(k\right)}^{\left(k\right)}=0$ happens at some instant $\tau \left(k\right)$; note that $\left\{\tau \right(k\left)\right\}$ are iid random variables distributed as a generic regeneration period $\tau$ of RA system;
• repeat previous step until the event $\tau \left(j\right)>{\tau }_e$ happens in some sample $V^{\left(j\right)}=(V_1^{\left(j\right)},V_2^{\left(j\right)},\dots )$; and,
• the value $V_{{\tau }_e}^{\left(j\right)}$ of the workload vector $V^{\left(j\right)}$ at instant ${\tau }_e$, has the target steady-state distribution of the workload in the original $M/G/c$ system.

We note that, although this approach allows to sample exactly from the steady-state distribution, the regeneration period in the dominated RA system can be very large in practice, and, thus, can lead to unacceptable long simulation. For further details on perfect sampling, see [15,35,36,37].

Now we explain how to sample from the equilibrium distribution of a two-component mixture. Let ${\overline{F}}_M$ be the tail of a two-component mixture Equation (3). Subsequently, it follows from Equation (3) that

$${\overline{F}}_M^{\left(e\right)}\left(x\right)=\frac{1}{\mathrm{E}{X}_M}{\int }_x^{\infty }{\overline{F}}_M\left(u\right)du=\frac{p\mathrm{E}{X}_1}{\mathrm{E}{X}_M}{\overline{F}}_1^{\left(e\right)}\left(x\right)+\frac{(1-p)\mathrm{E}{X}_2}{\mathrm{E}{X}_M}{\overline{F}}_2^{\left(e\right)}\left(x\right).$$

(32)

It is clearly seen from Equation (32) that the equilibrium distribution of a mixture is itself a two-component mixture of equilibrium distributions of the components. Thus, to sample from the equilibrium distribution (32), we sample from $F_1^{\left(e\right)}$ w.p. $q=p\mathrm{E}{X}_1/\mathrm{E}{X}_M$, and sample from $F_2^{\left(e\right)}$ w.p. $1-q$. Finally, note that, as easy to see, if the original distributions are $Pareto({\alpha }_i,x_0)$, then ${\overline{F}}_i^{\left(e\right)}$ are also $Pareto({\alpha }_i-1,x_0),i=1,2$ (also see [38]).

5. Simulation Results

As a sanity check of the perfect sampling $M/G/c$ model, we validate the algorithm via the $M/M/c$ system having input rate $\lambda =7.5$, service rate $\mu =1.5$, $c=10$ servers, and $\rho =\lambda /\mu =5$. We run $N=5000$ samples from steady-state process using perfect simulation and build the empirical queue size distribution vs. theoretical values that were obtained from Equation (22). We depict the results of validation on Figure 2. Note that the uniform distance between the empirical and theoretical distributions is $0.0091$.

5.1. Experiment 1: Hyperexponential Case

Now, we step away from the basic Markovian case, $M/M/c$, having service time distribution $F_1\left(x\right)=1-e^{-{\mu }_{1}x}$, by introducing a contaminated system $M/G/c$ with generic service time, $S^{\left(M\right)}$, having two-state hyperexponential distribution, $H_2$, with $F_i\left(x\right)=1-e^{-{\mu }_{i}x},$ and mixing coefficient $p\in (0,1)$. Note that such a case has computationally tractable solution, see [39]. However, we use the perfect sampling algorithm to check the accuracy of the sensitivity analysis. We fix

$${\mu }_1=2,c=\lambda =5,p=0.7,$$

and vary ${\mu }_2$ over range $(2,8]$ with step $0.4$. We obtain the empirical queue size d.f., ${\widehat{F}}_{Q^{\left(1\right)}}$, in the basic, and ${\widehat{F}}_{Q^{\left(M\right)}}$ in the contaminated system, and construct Equation (27) for each combination of the parameters while using N = 10,000 samples from the steady-state distribution. The linear dependence of $\mathsf{\Delta }({\widehat{F}}_{Q^{\left(M\right)}},{\widehat{F}}_{Q^{\left(1\right)}})$ on $\mathsf{\Delta }(F_M,F_1)$ is clearly seen in Figure 3.

5.2. Experiment 2a: Pareto Case, Sensitivity to Mixing Parameter

In the following experiments, we use an $M/G/c$ system with $c=4$, load $\rho =0.5$, and $Pareto({\alpha }_1,1)$ service time d.f., with ${\alpha }_1=2.1$ as the basic system for comparison. The input rate of the basic system is taken as $\lambda =\rho c({\alpha }_1-1)$ so as to guarantee the desired load $\rho =0.5$. Note that, to the best of our knowledge, there is no explicit expression for the steady-state queue size in such a system, and thus simulation is used to obtain the corresponding estimates of the steady-state queue size d.f. To obtain such an estimate, N = 10,000 samples from the corresponding steady-state distribution are obtained by the perfect sampling technique described in Section 4.

In the first experiment, we study the steady-state queue size distribution sensitivity to the mixing parameter, p. The mixing coefficient is iterated over the discrete values $p=0.95,0.9,\dots ,0.25$, and the empirical steady-state queue size d.f., $F_{Q_p^{\left(M\right)}}$, is constructed for the disturbed system with mixture service time d.f., $F_M$ given in Equation (3), consisting of $Pareto({\alpha }_1,1)$ and $Pareto({\alpha }_2,1)$ with mixing parameter p, where ${\alpha }_2=4.9$. The input rate $\lambda$ is fixed at the level $\lambda =2.2$, so as to guarantee the load $\rho =0.5$ in the basic system. Note that the parameter p is varied in such a way that the mixing proportion of $Pareto({\alpha }_2,1)$ distribution becomes larger with smaller p, and dominates the $Pareto({\alpha }_1,1)$, for $p<0.5$. Finally, we plot the values $\mathsf{\Delta }(F_M,F_1)$ vs. $\mathsf{\Delta }({\widehat{F}}_{Q_p^{\left(M\right)}},{\widehat{F}}_{Q^{\left(1\right)}})$ for the values p given. The results are depicted in Figure 4. Note that the dependence of the distance is approximately linear in mixing probability, p.

5.3. Experiment 2b: Pareto Case, Sensitivity to Contaminating Distribution

In the following experiment, we study the sensitivity of the steady-state queue size distribution on the parameter ${\alpha }_2$ of the mixture. Now $p=0.7$ is fixed, and ${\alpha }_2$ is iterated over the discrete set ${\alpha }_2\in \{2.1,2.3,\dots ,4.9\}$, ceteris paribus. As in the previous experiment, we build the empirical steady-state queue size distribution of the basic system, ${\widehat{F}}_{Q^{\left(1\right)}}$, by exact sampling from steady state using the method described in Section 4. We plot the values $\mathsf{\Delta }(F_M,F_1)$ vs. $\mathsf{\Delta }({\widehat{F}}_{Q_{{\alpha }_2}^{\left(M\right)}},{\widehat{F}}_{Q^{\left(1\right)}})$ for the given values of ${\alpha }_2$ as a parametric functions of ${\alpha }_2$. The results are depicted in Figure 5, where, unlike the previous scenarios, the nonlinear dependence on ${\alpha }_2$ is clear.

Note that the non-linear dependence of $\mathsf{\Delta }(F_{Q^{\left(M\right)}},F_{Q^{\left(1\right)}})$ on ${\alpha }_2$ may be caused by the non-linear dependence of the distance of service time distributions, $\mathsf{\Delta }(F_M,F_1)$, on ${\alpha }_2$, see Figure 1. Moreover, the mean service time, $S^{\left(M\right)}$, also differs from mean service time of the basic system, which causes appropriate changes in the load, $\rho$, in the disturbed system, see Figure 6.

Using the results of Experiment 2b, we illustrate the stochastic monotonicity property Equation (25) for selected values of parameter ${\alpha }_2$. Figure 7 depicts the results.

5.4. Experiment 2c: Pareto Case, Constant Load

In the final experiment, we study the joint effect of both the parent and the contaminating (Pareto) distributions. To do so, we change ${\alpha }_1=2.1,2.5,\dots ,4.9$, vary ${\alpha }_2={\alpha }_1,{\alpha }_1+0.4,\dots ,4.9$. To mitigate the effect of changing load illustrated by Figure 6, we simultaneously change the parameter $\lambda$, so as to guarantee constant load $\rho =0.5$ for all systems, keeping $p=0.7,c=4$ constant. The comparison is done to the system with the parent distribution of ${\alpha }_1=2.1$ of the service times. Each point is obtained then by N = 10,000 samples by the perfect sampling technique. Figure 8 depicts the results, where the color reflects the parent distribution parameter, ${\alpha }_1$, and size of a dot is proportional to ${\alpha }_2$. With increasing distance of the parent distribution from the contaminating distribution, the distance changes in a linear manner. Moreover, increased ${\alpha }_1$ changes the starting point (which in all lines corresponds to the parent distribution with parameter ${\alpha }_1$), and increasing ${\alpha }_2$ for fixed ${\alpha }_1$ increases the distance both in the input (for the mixture) and performance index (queue size distribution distance). Interestingly, for the lower line that corresponds to the fixed ${\alpha }_1=2.1$ and varying ${\alpha }_2$, there seems to be a slightly negative slope, which is likely to be the result of an increasing variance and, hence, decreasing accuracy. However, this effect might be interesting to study separately in the future.

Finally, we note that, to speedup the computation, we used parallel computation of the uniform distance for various system configurations using the resources of the High-Performance Datacenter of Karelian Research Centre of Russian Academy of Sciences.

6. Conclusions and Discussion

In this paper, the effect of the service time distribution perturbation on the steady-state performance measures of a multiserver queueing system is studied. The explicit form for the sensitivity measure (Kolmogorov-Smirnov distance) between the service time distribution functions was obtained, and the performance estimates were obtained by the regenerative perfect simulation technique. The simulation results outline the qualitative nature of the sensitivity, which is, in most cases, linear (possibly after appropriate scaling of the input rate to guarantee the constant load).

The approach to sensitivity analysis that is presented in this paper can be applied to more sophisticated, and more practically oriented systems, such as the simultaneous service multiserver system [40], which would result, though, in an increased dimension of the system state. However, we note that the steady-state exact sampling by regenerative simulation has several serious drawbacks. First, the average working time of the algorithm may be infinite [36], e.g., in a system with large number of servers (which indeed depends on the regenerative cycle length). This problem can be solved either by the coupling-from-the-past technique [35] (which, although, is rather technically tricky), or by non traditional regenerative techniques, such as the artificial regeneration [41] or regenerative envelopes [40]. Finally, the study may be extended to larger classes of service time distributions. At that, we leave these as opportunities for future research.