Communications in Statistics - Simulation and

Computation

ISSN: 0361-0918 (Print) 1532-4141 (Online) Journal homepage: https://www.tandfonline.com/loi/lssp20

Propagation of epistemic uncertainty in queueing

models with unreliable server using chaos
expansions

Katia Bachi, Cédric Chauvière, Hacène Djellout & Karim Abbas

To cite this article: Katia Bachi, Cédric Chauvière, Hacène Djellout & Karim Abbas (2019):
Propagation of epistemic uncertainty in queueing models with unreliable server using
chaos expansions, Communications in Statistics - Simulation and Computation, DOI:
10.1080/03610918.2019.1577966

To link to this article: https://doi.org/10.1080/03610918.2019.1577966

Published online: 11 Mar 2019.

Submit your article to this journal

Article views: 6

View Crossmark data

Full Terms & Conditions of access and use can be found at

https://www.tandfonline.com/action/journalInformation?journalCode=lssp20
COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R

https://doi.org/10.1080/03610918.2019.1577966

Propagation of epistemic uncertainty in queueing models

with unreliable server using chaos expansions
Katia Bachia , C
edric Chauviereb , Hacene Djelloutb , and Karim Abbasa
a
Operational Research Department, Research Unit LaMOS (Modeling and Optimization of Systems),
Faculty of Exact Sciences, Campus of Targua Ouzemour, Bejaia, Algeria; bLaboratoire de Math
ematiques
Blaise Pascal (LMBP), CNRS UMR 6620, Universit
e Clermont Auvergne, Campus Universitaire des
C
ezeaux, Cedex

ABSTRACT ARTICLE HISTORY

In this paper, we develop a numerical approach based on chaos Received 26 April 2018
expansions to analyze the sensitivity and the propagation of epi- Accepted 18 January 2019
stemic uncertainty through a queueing systems with breakdowns.
KEYWORDS
Here, the quantity of interest is the stationary distribution of the
Chaos expansions; Epistemic
model, which is a function of uncertain parameters. Polynomial uncertainty; Orthogonal
chaos provide an efficient alternative to more traditional Monte polynomial; Sobol’ indices;
Carlo simulations for modeling the propagation of uncertainty arising Unreliable queueing model
from those parameters. Furthermore, polynomial chaos expansion
affords a natural framework for computing Sobol’ indices. Such indi- AMS 2010 SUBJECT
ces give reliable information on the relative importance of each CLASSIFICATIONS
uncertain entry parameters. 90B22; 60K25; 65C05; 65C50

1. Introduction
Queueing systems with breakdowns are widely used to model problems occurring in
computer, manufacturing systems, and communication networks. Typically, in queueing
models the performances measurements are assessed for deterministic parameters val-
ues. However, in practice, the exact values of these parameters are not well known
(uncertain), and they are generally estimated empirically from few experimental obser-
vations. This lack of precise information propagates uncertainties in the output measure
and such study is the purpose of this work.
Since the pioneering work of Thiruvengadam (1963), Avi-Itzhak and Naor (1963),
many authors have been investigating the queueing systems with server breakdowns, see
for example Cao and Cheng (1982), Li, Shi, and Chao (1997) and Wang, Cao, and Li
(2001) and references therein. In this paper, we choose a functional approach for the
analysis of the dependence of the performance measure of certain queues, such as the
M/M/1/N and the M/G/1/N queue with breakdowns, with respect to some input param-
eters. More specifically, denoting the probability of a server breakdown by h, we seek to
compute the stationary distribution of the queue-length process, denoted by ph . For a
fixed value of h and a finite waiting capacity, ph can be found numerically by solving

CONTACT Hacene Djellout [email protected] Laboratoire de Math
ematiques Blaise Pascal (LMBP), CNRS
UMR 6620, Universit
e Clermont Auvergne, Campus Universitaire des C
ezeaux, France 3 place Vasarely, TSA 60026, CS
60026, 63 178 Aubiere Cedex.
Color versions of one or more of the figures in the article can be found online at www.tandfonline.com/lssp.
ß 2019 Taylor & Francis Group, LLC
2 K. BACHI ET AL.

P
ph Nh ¼ ph and ph ¼ 1, where Nh denotes the transition probability matrix of an
embedded jump chain. The definition of the embedded jump chain will depend on the
type of the queue: for the M/M/1/N queue, we use the sample-chain, embedded at
appropriate Poisson-times, and for the M/G/1/N queue, we will embed the chain at
departure and repair moments; details will be provided later in the paper. In case of a
large or infinite waiting capacity, ph can be obtained via a Laplace-Stieltjes transform,
see Abate, Choudhury, and Whitt (2000). Solving for ph involves numerical inversion,
which can be computationally demanding, see Shortle, Fischer, and Brill (2007). The
computation of ph is a challenging problem and a variety of approaches have been pro-
posed in the literature for approximately or indirectly solving the stationary distribution.
The predominant approach is to obtain either the generating function of ph or an ana-
lytical expression for ph containing a Laplace-Stieltjes transform, see, for example,
Abate, Choudhury, and Whitt (2000) and Baccelli and Znati (1981). Also numerical sol-
utions by means of the matrix geometric method (Neuts 1981) are available, see
Mitrany and Avi-Itzhak (1968) and Neuts and Lucantoni (1979) for details.
In performance analysis, one is not only interested in evaluating the system for spe-
cific set of parameters but also in the sensitivity of the performance with respect to
these parameters. For example, in a queueing model with breakdowns, the breakdown
probability is a key parameter of interest and, in this paper, we will analyze the depend-
ence of ph on h, which is significantly more challenging than evaluating ph for a fixed
h. Most of the time, it is assumed that these stochastic models are solved for fixed
parameters values. However, the parameters of the model are determined through insuf-
ficient statistical data (a limited number of observations), leading to uncertainty in the
assessment of their values. This parametric uncertainty, induced from the incomplete
information of the parameter, is called epistemic uncertainty (Limbourg 2008; Mishra
and Trivedi 2011; Winkler 1996).
In order to estimate the uncertainty of the parameters in the performance measure of
the model, two complementary approaches may be used: uncertainty analysis and sensi-
tivity analysis. Uncertainty analysis consists in modeling input parameters as random
variables, see, for example, Cacuci (2003) and Helton (1994, 1997). Then, the sensitivity
analysis aims at determining the relative contribution of individual parameters. Several
approaches for the propagation of uncertainty have been developed, including interval
arithmetic (Moore 1979; Rocco and Klindt 1998), Taylor series expansion (Dhople and
Dom
ınguez-Garc
ıa 2012; Granger and Henrion 1993; Ouazine and Abbas 2016;
Shooman 1990; Takhedmit and Abbas 2017), moments (Granger and Henrion 1993;
Mattia, Sergio, and Dimitrov 2007; Shooman 1990) and Monte Carlo analysis (Granger
and Henrion 1993; Mattia, Sergio, and Dimitrov 2007; Shooman 1990). Overviews of
these approaches are available in several reviews (Winkler 2004; Helton 1993; Helton
et al. 2006; Iman and Helton 1988; Ionescu-Bujor and Cacuci 2004; Ronen 1988).
Since the seminal work of Ghanem and Spanos (1991) and later generalized by Xiu
and Karniadakis (2002), polynomial chaos has spread over a broad scientific commu-
nity, see Sepahvand, Marburg, and Hardtke (2010) and Chauviere, Hesthaven, and
Lurati (2006). The interested reader may also refer to Shuxing, Xiong, and Wang (2017)
for improved versions of the polynomial chaos technique to reduce computational cost.
Then, due to the specific form of these polynomials, it quickly became apparent that
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
3

this representation was suitable to compute Sobol’ indices (Sudret 2008). To the best of
our knowledge, this is the first attempt to use chaos expansion to investigate the propa-
gation of the epistemic uncertainty through the performance measure of the queueing
systems using Sobol’ indices. Others studies have been carried out to model and propa-
gate uncertainty in software systems, see Aleti et al. (2018), for example.
The main objective of this work is to develop a numerical procedure to investigate
the sensitivity and the propagation of epistemic uncertainty in the input parameters of a
queueing systems with unreliable server. As a typical example, a customer arrives at the
queueing system according to a Poisson process with rate k, which we consider as an
input parameter of the model. We assume that this parameter is not precisely known
and therefore it can be modeled as a random variable of known probability density
function. Clearly, the output measure is a function of this input random variable. The
epistemic uncertainty is propagated through a functional relationship of the type
YðxÞ ¼ gðkðxÞÞ which links the input parameter to the output measures. In this work,
special attention is given to the output stationary distribution. From a statistical point
of view the quantities of interest may be the moments of the functional stationary dis-
tribution. Such quantities can easily be computed from the coefficient of the polynomial
chaos representation of gðkðxÞÞ.
Over the past few years, several techniques have emerged to compute the coefficients of
polynomial chaos. They can be classified in two main families: intrusive methods and non-
intrusive methods (Pettersson, Iaccarino, and Nordstr€ om 2015). In our context, YðxÞ can
be explicitly expressed as a functional relationship and therefore the two approaches are of
equal complexity. In this paper, projection method based on the Gauss quadrature rules
are used to obtain the chaos coefficients. Then statistical moments and Sobol’ sensitivity
indices can easily be computed from polynomial chaos expansion.
The paper is organized as follows. The next section, is devoted to restate some basic
notions about orthogonal polynomial and chaos expansion. Then, the Sobol’ indices
which are an important tool in sensitivity analysis are introduced. Epistemic uncertainty
in input parameters of two queueing systems with unreliable server is considered in
Sec. 4. In the first queue, a single input parameter is considered as random. On the
other hand, the second queueing system considers four input random variables and a
sensitivity analysis is performed. For the two queues, numerical results for the computa-
tion of statistical moments are obtained using both chaos expansion and Monte-Carlo
simulation. Finally, Sec. 5 draws some conclusions.

2. Polynomial chaos expansion

Polynomial chaos expansion shares many elements with power series involving orthogonal
polynomials, which are commonly used in the spectral community. This is why the follow-
ing subsection is devoted to restate some results concerning orthogonal polynomials.

2.1. Orthogonal polynomials: univariate case

Let V be the real vector space of all polynomials in one variable with real coefficients
and let <; > be the positive inner product on V defined as
4 K. BACHI ET AL.

ð
<u; v> ¼ uðxÞvðxÞf ðxÞdx 8u; v 2 V;
I

where f : I  R ! Rþ is a nonnegative integrable function of x hereinafter referred to

as a ‘weight function’. The set of polynomials fWn gn0 are said to be orthogonal with
respect to the weight function f(x), if Wn is a polynomials of degree n and
ð
<Wn ; Wm > ¼ Wn ðxÞWm ðxÞf ðxÞdx ¼ h2n dn;m ; n; m 2 N; (2.1)
I

where dn;m is the Kronecker’s delta function, and hn are non-zero constants. We recall
that, for orthogonal polynomials of degree d ¼ 0, W0 is always equal to one (W0 ¼ 1Þ.
Furthermore, the system (2.1) is called orthonormal if hn ¼ 1.
The most general way to build such polynomials rely on the three-term recurrence
relation
Wnþ1 ðxÞ ¼ ðxan ÞWn ðxÞbn Wn1 ðxÞ;
(2.2)
W0 ðxÞ ¼ 1; W1 ðxÞ ¼ 0;
with
<xWn ; Wn >
an ¼ ; n2N
< Wn ; Wn >
<Wn ; Wn > : (2.3)
bn ¼ ; n 2 N
< Wn1 ; Wn1 >

(see Gautschi 2004 for more details). Furthermore, according to Favard’s theorem, for a
given weight function f, corresponds a unique set of coefficients ðan ; bn Þn2N . From a
numerical point of view, the integrals appearing into (2.3) can be evaluated using Fej
er
quadrature rule, see Rahman (2009), for example. Then, taking an orthogonal polyno-
mial basis of degree n, fW0 ; W1 ; :::; Wn g; any sufficiently regular function u : I  R ! R
may be represented by its projection Pn u on such basis, i.e.
Xn
Pn uðxÞ ¼ ^ i Wi ðxÞ;
u (2.4)
i¼0
<u;Wi >
where the coefficients of the projection can be computed by evaluating u ^ i ¼ <Wi ;Wi >
.
Depending on the regularity of the function u and the choice of the polynomial basis,
upper bounds of the truncation error ku  Pn uk can be derived (Gottlieb and Orszag
1977). In particular, when u is C1 , one can observe the so-called spectral convergence
(the truncation error decays exponentially fast with respect to n).
The construction of multivariate orthogonal polynomial basis rely on univariate poly-
nomial basis as will be described in the following section.

2.2. Orthogonal polynomials: multivariate case

The setting is now the real vector space Vd of all polynomials in d variables with real
coefficients, together with the positive inner product on Vd defined as
ð
<u; v> ¼ uðxÞvðxÞf ðxÞdx 8u; v 2 V d ; (2.5)
I
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
5

where, now, f : I  Rd ! Rþ is a nonnegative integrable function of x ¼ ðx1; x2; :::; xd Þ.

We further assume that the weight function f can be written as the product of univari-
ate weight functions, i.e. f ðx1; x2; :; xd Þ ¼ f1 ðx1 Þf2 ðx2 Þ:::fd ðxd Þ, where fi : Ii  R ! Rþ are
nonnegative integrable function. It would be tempting to construct the elements of the
multivariate polynomial basis in a similar fashion i.e. as a tensor product of monovari-
ate elements. However, such construction leads to uncomplete basis (all monomials up
to a certain degrees are not represented). To ensure that fW0 ðxÞ; W1 ðxÞ; :; WPp ðxÞg is a
p d
multivariate polynomial basis of degree p, with ðPd þ 1Þ elements, the procedure
described in the sequel is usually adopted. We first define Wa1 :::ad ðxÞ as the tensor prod-
uct of the elements of a univariate polynomial basis, of degree n, i.e.
Y
d
Wa1 :::ad ðxÞ ¼ Wai ðxi Þ; (2.6)
i¼1

where ai 2 f0; 1; :::; ng. However, not all the elements of the form (2.6) are retained
when constructing the multivariate polynomial basis (that would lead to an uncomplete
d
Pd ðn þ 1Þ elements). Instead, for a given degree p, only the ele-
basis of degree nd with
ments that satisfy i¼0 ai  p in (2.6) are kept, and a one to one correspondence
between the multi-index ða0 ; :::; ad Þ and the ith element Wi ðxÞ of the multivariate basis
is set. Proceeding that way, it can be shown that the number of elements of a complete
multivariate polynomial basis of degree p is
 
p pþd ðp þ dÞ!
Pd þ 1 ¼ ¼ : (2.7)
d d!p!
Similarly to (2.4), any function u : Rd ! R may be represented by its projection Pp u
on such basis, i.e.
p
X
Pd
Pp uðxÞ ¼ ^ i Wi ðxÞ:
u (2.8)
i¼0

Having set the basic framework of orthogonal polynomials, we now look into the
way they can be efficiently used in the field of probability, where they are usually
referred as Polynomial Chaos (PC).

2.3. Polynomial chaos expansion

Denote by ðX; A; PÞ the probability space, where as usual X is the set of all possible outcomes,
A is a r-algebra over X; and P is a function A ! ½0; 1 that gives a probability measure on
A. Consider an Rd -valued independent random vector X ¼ ðX1 ; :::; Xd Þ that describes input
uncertainties. The probability law of X may be defined by the probability density function
Y
d
fX ðxÞ ¼ fi ðxi Þ; (2.9)
i¼1

where fi ðxi Þ is the marginal probability density of Xi defined on ðXi ; Ai ; P i Þ. Let us now
denote L2 ðXi ; Ai ; P i Þ the space of real random variables with finite second order
moments, i.e. such that
6 K. BACHI ET AL.

ð
 
E Xi2 ¼ xi2 fi ðxi Þdxi <1; (2.10)

where E stands for the mathematical expectation. L2 ðXi ; Ai ; P i Þ is a Hilbert space which
can be provided with a set of complete orthogonal basis fWij ðxÞgj0 that are consistent
with the density of Xi. For example, Legendre polynomials are associated with uniform
distribution whereas Hermite polynomials are associated with Gaussian distribution.
Similarly, L2 ðX; A; PÞ is provided with a set of complete multivariate orthogonal basis
fWj ðxÞgj0 which, in turn, is consistent with the density of X.
Let Y ¼ ðY1 ; :::; Y‘ Þ : X ! R‘ such that Yi 2 L2 ðX; A; PÞ, for i ¼ 1; :::; ‘, be a math-
ematical model. For the sake of simplicity, we will consider only one component of this
model, denoted by Y, since the same procedure apply identically to all the other compo-
nents. Since Y is assumed to belong to L2 ðX; A; PÞ, it can be represented as (Ghanem
and Spanos 1991, Xiu and Karniadakis 2002)
X
1 X
1 X
i1
Y ðX1 ; :; Xd Þ ¼z0 W0 þ zi1 W1 ðXi1 Þ þ zi1 i2 W2 ðXi1 ; Xi2 Þ
i1 ¼1 i1 ¼1 i2 ¼1
(2.11)
1 X
X i1 X
i2
þ zi1 i2 i3 W3 ðXi1 ; Xi2 ; Xi3 Þ þ    ;
i1 ¼1 i2 ¼1 i3 ¼1

which, after some rearranging, can be rewritten in a more convenient way as

X
1
Y ðXÞ ¼ yj Wj ðXÞ: (2.12)
j¼0

Then, similarly to Eq. (2.8), this series is truncated by keeping terms up to a degree p
p
X
Pd
Y ðXÞ Pp Y ðXÞ ¼ yj Wj ðXÞ: (2.13)
j¼0

Example 2.1. Let X ¼ ðXi1 ; Xi2 Þ be a random normal vector. The functional approxima-
tion of the response Y ¼ YðXÞ may be approximated with Hermite-Chaos expansions
according to (2.13). If we fix the degree of the chaos expansions to be d ¼ 2, we obtain the
series expansion with P22 þ 1 ¼ 6 terms as follows:
X
5
Y yi Wi ðX1 ; X2 Þ
i¼0 (2.14)
¼ y0 W00 þ y1 W10 þ y2 W01 þ y3 W11 þ y4 W20 þ y5 W02 ;
where Wi1 i2 ðX1 ; X2 Þ ¼ Wi1 ðX1 ÞWi2 ðX2 Þ is the product of univariate Hermite polynomials of
degree i1 and i2 satisfying i1 þ i2  2. Such polynomials can easily be computed according
to the recurrence formula (2.2) together with (2.3) with an ¼ 0 and bn ¼ n. Should we
have used Legendre polynomials instead of Hermite’s ones, they would be built by taking
2
an ¼ 0 and bn ¼ ð4nn2 1Þ. Table 1 shows the bivariate Hermite polynomial basis con-
structed from univariate basis f1; X1 ; X12 1g and f1; X2 ; X22 1g and the link between the
main index j and the multi-index ði1 ; i2 Þ. For completeness, the same table also provides
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
7

Table 1. Bivariate Hermite and Legendre polynomial basis.

j ði1 ; i2 Þ Wi1 i2 (Hermite) Wi1 i2 (Legendre)
0 (0, 0) 1 1
1 (1, 0) X1 X1
2 (0, 1) X2 X2
3 (1) X1 X2 X1 X2
4 (2, 0) X12 1 X12 1=3
5 (0, 2) X22 1 X22 1=3

bivariate Legendre polynomial basis constructed from univariate basis f1; X1 ; X12 1=3g
and f1; X2 ; X22 1=3g.

2.4. Computing PC coefficients

The way to compute the PC coefficients fyj g0jPp appearing into Eq. (2.13) can be
d
casted into two different families: projection methods and regression methods
Pettersson, Iaccarino, and Nordstr€ om (2015). Here we have used the first one which
consists in premultiplying (2.13) by Wj ðXÞ and by taking the expectation of the resulting
product. Using the orthogonality of the PC basis, most of the terms cancel and we end
up with
ð
 
ð Þ ð
yj ¼ E Y X Wj X ¼ Þ p
Y ðxÞWj ðxÞf ðxÞdx for j ¼ 0; 1:; Pd : (2.15)
IRd

The above integral can be evaluated through different techniques: going from rough
Monte-Carlo sampling simulation to Gaussian quadrature rule or sparse quadrature
rules when the dimension d of the input random is high. Here, we evaluate such inte-
grals using Gaussian quadrature rules which take the form
ð X X X
Ng1 Ng2 Ngd
Y ðxÞf ðxÞdx :: xi1 xi2 :xid Y ð~x i1 ; ~x i2 :; ~x id Þ; (2.16)
i1 ¼1 i2 ¼1 id ¼1
IRd

where fxik g1ik Ng are the quadrature weight and f~x ik g1ik Ng are the quadrature
k k
points. The Gaussian quadrature rules are such that the evaluation of the integral is
exact if Y is a multivariate polynomial containing monomials xik of maximum degree
2Ngk 1. From Eqs. (2.15) and (2.16), we see that it suffices to evaluate the response Y
at Ng1 Ng2 ::: Ng‘ deterministic quadrature points in order to compute the PC coeffi-
cients. Furthermore since W0 ¼ 1, setting j ¼ 0 in Eq. (2.15) leads to
y0 ¼ EðY ðXÞÞ; (2.17)
i.e. the first coefficient of the PC expansion is the expectation of the random response
of the system. Similarly, by considering the approximation of Y2
p p

2
X X
Pd Pd
Y ðXÞ yi yj Wi ðXÞWj ðXÞ; (2.18)
i¼0 j¼0

and taking the expectation on each side, the orthogonality of the PC basis leads to
8 K. BACHI ET AL.

p
 2
 X
Pd
ð Þ
E Y X ¼ yi2 ; (2.19)
i¼0

from which the variance of the random response of the system can easily be deduced.
Furthermore, the PC decomposition provides a convenient way of computing Sobol’
indices, as explained in the next section.

3. Sensitivity analysis
The purpose of sensitivity analysis is to investigate the influence of each input parameter
and their possible interactions onto the output measures. They can be casted into two
main families: local analysis based on a local perturbation around an average value and
global analysis that consider input parameters as random variables and decompose the
output variance into several components. The Sobol’ indices belong this last type of family.

3.1. Sobol’ indices

As in the previous section, we consider the mathematical model Y ¼ YðXÞ, where the
input parameter X ¼ ðX1 ; :::; Xd Þ are d independent random variables belonging to
L2 ðXi ; Ai ; P i Þ for i ¼ 1; :; d, and similarly, Y is assumed to belong to L2 ðX; A; PÞ. Sobol
(1993) proposed an indicator of the influence of the input parameter Xi defined by
V ðEðY=Xi ÞÞ Vi
Si ¼ ¼ ; (3.1)
V ðY Þ V
commonly termed ‘first order Sobol’ indices’. Vi ¼ VðEðY=Xi ÞÞ is the conditional vari-
ance of Y with respect to Xi and V ¼ VðYÞ is the total variance of Y. Similarly, sensitiv-
ity indices of higher order can be defined by first introducing the following
decomposition of the total variance, which is valid for any output Y belonging to
L2 ðX; A; PÞ
X
d X X
V¼ Vi þ Vi1 i2 þ Vi1 i2 i3 þ    þ Vi1 :::id ; (3.2)
i1 ¼1 1i1 <i2 d 1i1 <i2 <i3 d

where
  
Vi1 ¼ V E Y=Xi1 ;
  
Vi1 i2 ¼ V E Y=Xi1 ; Xi2 Vi1 Vi2 ;
  
Vi1 i2 i3 ¼ V E Y=Xi1 ; Xi2 ; Xi3 Vi1 i2 Vi1 i3 Vi2 i3 Vi1 Vi2 Vi3 ;
(3.3)
;
X
d X X
Vi1 :::id ¼ V Vi  Vi1 i2      Vi1 :::id1 :
i¼1 1i1 <i2 d 1i1 <i2 :::<id1 d

Then, Sobol’ indices of order k are given by

Vi1 :::ik
Si1 :::ik ¼ : (3.4)
V
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
9

Although Sobol’ indices could be computed by estimating integrals appearing into

Eq. (3.3), such a procedure would be both computationally expensive and hardly tract-
able. Instead, computing those indices from the PC representation of the random out-
put Y turns out to be an efficient alternative and has been the method of choice for
many years since the seminal work of Sudret (2008). For that, we define by Ii1 ;i2 ;:::;is
(s  d) the set of d-dimensional vectors a ¼ ða1 ; :::; ad Þ with a 6¼ 0 and 0  a1  ::: 
ad  d that selects elements of the PC basis Wa1 :::ad defined by (2.6) containing solely
the variables Xi1 ; Xi2 ; :::; Xis .
This way, the multi-indices defined by Ii will select elements of the PC basis depend-
ing only on the variable Xi . Similarly, Ii;j will select elements of the PC basis that
depend on Xi and Xj ; to the exclusion of any other variable. With this notation, Sobol’s
indices of order k are simply given as a function of the PC coefficients as follows:
1 X
Si1 :::ik ¼ y2 : (3.5)
V ða1 ;:::;ad Þ2Ii ;:::;i a1 ;:::;ad
1 k

For a problem with d input random parameters, it can be shown that ð2d 1Þ Sobol’
indices can be computed for each output random quantity of interest. When the num-
ber of input r.v. is large, the number of Sobol’ indices grows exponentially and it
becomes difficult to draw information from these statistics. This is why, Homma and
Saltelli (1996) introduced the total sensitivity indices STi ði ¼ 1; :::; dÞ which measures
the total effect of the ith random input parameter. It is defined as the sum of all sensi-
tivity indices Si1 :::ik ðk ¼ 1; :::; dÞ for which, one of the indices i1 ; i2 ; :::; ik is equal to i,
i.e.
X
d
STi ¼ Si1 ;:::;ik ;
k¼1
ði1 ;:::;ik Þ2Jik

where Jik is the set of k-dimensional vectors ði1 ; :::; ik Þ with 1  i1 <:::<ik  d, such that
one of its components is equal to i. We can also compute these indices from the PC
coefficients as follows:
1 X
STi ¼ 1 y2 ;
V a ;:::;a 2I a1 ;:::;ad
ð 1 dÞ i

where I i is the complementary set of Ii .

Example 3.1. The computation of the sensitivity indices for the example 2.1 where,
p ¼ 2 and d ¼ 2 gives
1 2 
S1 ¼ y1;0 þ y2;0
2
;
V

1 2 
S2 ¼ y þ y0;22
; (3.6)
V 0;1
1 2
S1;2 ¼ y1;1 ;
V
with
V ¼ y1;0
2
þ y0;1
2
þ y1;1
2
þ y2;0
2
þ y0;2
2
: (3.7)
10 K. BACHI ET AL.

In accordance with (3.2), we see that the variance V of the output can be divided into
three parts: V ¼ S1 þ S2 þ S1;2 : The first part S1 represents the influence of the first input
random variable alone; the second S2 represents the effect of the second one alone whereas
S1;2 accounts for the combined effect of the two input random variables.
Similarly, according to the definition of the total sensitivity indices, we have
ST1 ¼ S1 þ S1;2 ;
and
ST2 ¼ S2 þ S1;2 :

4. Uncertainty analysis in unreliable queueing models

In this section, we consider the PC expansion for propagating the uncertainty in performance
measure of queueing models with breakdowns and repairs, due to epistemic uncertainties in
the model input parameters. Section 4.1 studies the M/G/1/N queue with breakdowns, when
the perturbation of a single parameter is introduced. Then, Sec. 4.2 addresses the case of the
M/M/1/N queue with breakdowns, for multiple random input parameters.

4.1. The M/G/1/N queue with breakdowns and repairs

Consider a finite capacity M/G/1/N queue with server subject to breakdowns and repairs.
Assume that the customers arrive at the system according to a Poisson stream with rate k,
the service times are general and independent identically distributed with mean 1=l and
we denote the service time distribution by F(x). There can be at most N customers in the
queue (including the one being served), and customers attempting to enter the queue
when there are already N customers present are lost. The service discipline is assumed to
be FCFS. At the beginning of each service, there is a probability h that the server breaks
down (and the customer is sent back to the queue) and enters into a repair state. The time
of repair is exponentially distributed with rate r. The only instant of time when a server
can breaks down is right at the beginning of a service. This system is modeled by a Markov
chain X ¼ fXn : n 2 Ng, embedded at service completions and completion of a repair and
has state space S ¼ f0; 1; :::; N1g. Let N ¼ ðai;j Þi;j2S be the transition probability matrix
of the Markov chain X. The matrix N has the following form:
0 1
X
N 2
B a0 a1 a2 a3    aN2 1 ak C
B C
B k¼0 C
B N 2 C
X
B C
B a0 a1 a2 a3    aN2 1 ak C
B C
B k¼0 C
B N 3 C
X
B C
B 0 a0 a1 a2    aN3 1 ak C
N¼B B k¼0
C;
C (4.1)
B X
N 4 C C
B
B 0 0 a0 a1    aN4 1 ak C
B C
B C
B k¼0
C
B .. .. .. .. C
B . . . . C
B C
@ 0 0 0 0  a0 1a0 A
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
11

where
ð1  k
kx
ðkxÞk r k
ak ¼ h e dF ðxÞ þ ð1hÞ ; k ¼ 0; :::; N2: (4.2)
0 k! rþk kþr
Note that the Markov chain X is unichain and therefore the stationary distribution
exists. Let p denote the stationary distribution of the queue-length process embedded at
service completions and completion of a repair in the M/G/1/N queue with breakdowns
and repairs (see Abbas, Berkhout, and Heidergott (2016) for details). In the sequel, we
consider p as a function of the breakdown probability h, denoted by pðhÞ. Assume that
the probability h is estimated from insufficient statistical data, and hence has uncer-
tainty associated with it. In the next section, we will discuss a functional approach based
on PC expansion for computing the epistemic uncertainty in stationary distribution
pðhÞ, due to epistemic uncertainties in h.

4.1.1. Epistemic uncertainty in queueing system

In this model, we will assume that the uncertainty affects only one parameter. To ana-
lyze the propagation of epistemic uncertainty, the mean and the variance of the
responses will be computed using the PC expansion.
Let us consider pðhÞ ¼ ðp0 ðhÞ; :::; pN1 ðhÞÞ be the random response or the model out-
put of the queueing system. We consider the following perturbation model on the
uncertain parameter:
hðxÞ ¼ ^h þ rh eðxÞ; (4.3)
where hðxÞ is a random variable with uniform density on the interval ½0; 1. The param-
eters ^
h and rh represent the mean and the standard deviation of h and they can be esti-
mated by statistical methods. eðxÞ can be considered as a random noise imposed to h
and it modelizes epistemic uncertainty. We assume that eðxÞ is a random variable with
uniform distribution on ½1; 1. Then, the approximation of any component p‘ ðhÞ of
the stationary distribution by the PC expansions writes
p
X
Pd
p‘ ðhÞ yi Wi ðhÞ;
i¼0

where fWi ðhÞg0iPp form an orthonormal Legendre polynomial basis, and yi are the
d
coefficients of the approximating series expansion, see (2.13).

4.1.2. Numerical results

In this section, the parameters N, k and r are deterministic and their values are set to
N ¼ 7, k ¼ 1, r ¼ 0.4, respectively. The only input random parameter is h and it is given by
h ¼ 0:5 þ 0:28 e; e U ð1; 1Þ:
For numerical computations, we propose two types of distribution for the service
time, namely, generalized Erlang ðE2 Þ and Hyperexponential of second order ðH2 Þ.
Those two distributions are described in the following subsections.
12 K. BACHI ET AL.

Table 2. Expected value of the steady-state vector in M/E2 /1/7 queue, together with its 95% confi-
dence interval.
Eðpi Þ PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 0.0237 0.0230 0.0237 [0.0229, 0.0245] [0.0236, 0.0237]
p1 0.0413 0.0422 0.0413 [0.0402, 0.0425] [0.0412, 0.0415]
p2 0.0660 0.0667 0.0660 [0.0646, 0.0674] [0.0658, 0.0661]
p3 0.1018 0.1026 0.1018 [0.1004, 0.1032] [0.1017, 0.1019]
p4 0.1557 0.1563 0.1558 [0.1549, 0.1566] [0.1556, 0.1558]
p5 0.2393 0.2401 0.2392 [0.2384, 0.2402] [0.2392, 0.2394]
p6 0.3722 0.3719 0.3722 [0.3673, 0.3771] [0.3717, 0.3727]

4.1.3. Epistemic uncertainty in the M/E2 /1/N queue with server breakdowns
and repairs
In (4.2), the density F 0 ðxÞ ¼ f ðxÞ of the service time is assumed to be generalized
Erlang of second order, i.e.
ll
f ðxÞ ¼ 1 2 ðel1 x el2 x Þ1½0;þ1Þ ðxÞ;
l2  l1
where the rates l1 ¼ 4 and l2 ¼ 2.
We determine the stationary distribution P vector with respect to the random variable
hðxÞ by solving the system ph Nh ¼ ph and ph ¼ 1, where Nh denotes the transition
probability given by (4.1). Therefore, the random outputs of interest are the stationary
distributions pi ; i 2 f0; :::; 6g: Their projection is performed on a monovariate PC basis
of degree n ¼ 4 and the coefficients are computed using a 6-points Gaussian quadrature
rule as explained in Sec. 2.4. Numerical results for the expectations and the variances
using a Monte-Carlo simulations of a sample size NMC ¼ 1000 and NMC ¼ 100000 are
given in Tables 2 and 3, respectively. Results are compared with those of the PC expan-
sion. We note that the Monte Carlo simulations (MC in short) converge to the PC
results; however the later are obtained at a fraction of the cost of the MC simulations.
Similarly moments of higher order (skewness and kurtosis) are given in Tables 4 and 5,
and the same conclusions can be drawn.
Figure 1 shows the density of the output stationary distributions for a random uni-
form input of the form (4.3). The shape of the densities show that they are far from
being uniform. Furthermore, the plots are coherent with Tables 4 and 5. Indeed, a posi-
tive skewness reflects a random variable with higher probabilities in the left part of its
support, and vice versa.
Regarding the computational times, numerical experiments are performed on a
3.30 GHZ IntelR coreTM i3 CPU processor with 4 GB memory in the MATLAB environ-
ment. Tables 6 and 7 show the execution times of the PC and the MC approaches. For
a similar accuracy, we see that the PC simulation are about fifty time faster than the
MC simulations.

4.1.4. Epistemic uncertainty in the M/H2 /1/N queue with server breakdowns
and repairs
In this section, the density F 0 ðxÞ ¼ f ðxÞ of the service time is now assumed to be
Hyperexponential of second order, i.e.
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
13

Table 3. Variance of the steady-state vector in M/E2 /1/7 queue, together with its 95% confi-
dence interval.
Vðpi Þ 104 PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 1.7753 1.6686 1.7809 [1.6322, 1.9458] [1.7598, 1.7909]
p1 3.6413 3.6789 3.6423 [3.3478, 3.9911] [3.6097, 3.6734]
p2 5.2349 5.1333 5.2175 [4.8130, 5.2676] [5.1894, 5.2812]
p3 5.0800 4.9922 5.0508 [4.6706, 5.5680] [5.0358, 5.1249]
p4 1.8737 1.8585 1.8732 [1.7227, 2.0537] [1.8574, 1.8902]
p5 2.1440 2.0101 2.1534 [1.9712, 2.3500] [2.1253, 2.1630]
p6 58.7699 61.0693 58.7199 [54.0336, 64.416] [58.2592, 59.2896]

Table 4. Skewness coefficient of the steady-state vector in M/E2 /1/7 queue, together with its 95%
confidence interval.
Skewðpi Þ PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 0.5053 0.5429 0.5064 [0.3534, 0.6571] [0.4901, 0.5204]
p1 0.3565 0.2911 0.3550 [0.2046, 0.5083] [0.3413, 0.3717]
p2 0.1730 0.1435 0.1701 [0.0211, 0.3248] [0.1578, 0.1881]
p3 0.0722 0.0927 0.0730 ½0:2241; 0:0796 ½0:0874;0:0571
p4 0.5848 0.6616 0.5922 ½0:7366;0:4330 ½0:6000;0:5696
p5 0.7140 0.7968 0.7036 ½0:8658;0:5622 ½0:7292;0:6988
p6 0.0543 0.0692 0.0540 ½0:0975; 0:2061 [0.0391, 0.0695]

Table 5. Kurtosis coefficient of the steady-state vector in M/E2 /1/7 queue, together with its 95%
confidence interval.
Kurtðpi Þ PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 2.0384 2.1425 2.0341 [1.7347, 2.3420] [2.0080, 2.0688]
p1 1.8968 1.9122 1.8953 [1.5931, 2.2004] [1.8664, 1.9271]
p2 1.7864 1.7682 1.7900 [1.4828, 2.0901] [1.7561, 1.8168]
p3 1.7463 1.7513 1.7495 [1.4427, 2.0499] [1.7159, 1.7767]
p4 2.0301 2.0290 2.0385 [1.7265, 2.3338] [1.9997, 2.0605]
p5 2.1921 2.3603 2.1972 [1.8884, 2.4957] [2.1617, 2.2224]
p6 1.7759 1.8007 1.7724 [1.4722, 2.0795] [1.7455, 1.8062]

 
f ðxÞ ¼ cl1 el1 x þ ð1cÞl2 el2 x 1½0;1Þ ðxÞ;

where the rates are l1 ¼ 3=2; l2 ¼ 3; and c ¼ 0:3. All the other parameters are set as in
the previous section, including the random input. Similarly, Tables 8–11 show the
moments of the stationary distributions when the Hyperexponential law is used instead
of the generalized Erlang one for the service time. Here also, we can observe the good
convergence properties of the PC, showing the robustness of the method independently
of the law of the service time.
Similarly, Figure 2 shows the density of the output stationary distributions for a ran-
dom uniform input of the form (4.3).

4.2. The M/M/1/N queue with server breakdowns and threshold-based

recovery policy
Consider a finite capacity M/M/1/N queue with server breakdowns and threshold-based
recovery policy. If a server provides service for a customer, the server may experience a
breakdown, and the repair begins when the number of customers present in the system
exceeds some prespecified threshold level q (1  q  N). The stream customers arrive
14 K. BACHI ET AL.

Figure 1. Marginal probability density fpi of the steady-state vector pðhÞ in M/E2 /1/7 queue.

Table 6. Execution time for the computation of pi using the PC approach and MC simulation of the
M/E2 /1/7.
pi n Runtime approach PC MC (NMC ¼ 105 )
p0 0.058433 1.184682
p1 0.022896 1.100644
p2 0.023910 1.089840
p3 0.024508 1.130768
p4 0.022965 1.101523
p5 0.023862 1.139988
p6 0.024123 1.135433

Table 7. Execution time for the computation of statistical moments of the model response M/E2 /1/7
using the two approaches PC and MC with NMC ¼ 105 .
Eðpi Þ Vðpi Þ Skewðpi Þ Kurtðpi Þ
Moments
Approach PC MC PC MC PC MC PC MC
p0 0.0218 1.0678 0.0236 1.0681 0.0237 1.1037 0.0237 1.1039
p1 0.0216 1.0744 0.0241 1.0747 0.0244 1.1102 0.0245 1.1104
p2 0.0228 1.0596 0.0251 1.0599 0.0253 1.0956 0.0254 1.0956
p3 0.0234 1.0874 0.0254 1.0876 0.0257 1.1232 0.0257 1.1234
p4 0.0231 1.0816 0.0250 1.0818 0.0250 1.1178 0.0251 1.1180
p5 0.0237 1.1186 0.0258 1.1189 0.0261 1.1544 0.0261 1.1547
p6 0.0266 1.0955 0.0279 1.0958 0.0279 1.1312 0.0280 1.1314

at the queue according to a Poisson law of parameter k. Arriving customers form a sin-
gle waiting line based on the order of their arrivals and the server can serve only one
customer at a time. There can be at most N customers present at the queue (including
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
15

Table 8. Expected value of the steady-state vector in M/H2 /1/7 queue, together with its 95% confi-
dence interval.
Eðpi Þ PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 0.0397 0.0407 0.0396 [0.0383, 0.0411] [0.0396, 0.0398]
p1 0.0602 0.0588 0.0602 [0.0586, 0.061] [0.0600, 0.0603]
p2 0.0853 0.0847 0.0853 [0.0836, 0.0870] [0.0851, 0.0854]
p3 0.1179 0.1179 0.1179 [0.1165, 0.1192] [0.1177, 0.1180]
p4 0.1619 0.1615 0.1619 [0.1613, 0.1624] [0.1618, 0.1619]
p5 0.2233 0.2238 0.2233 [0.2219, 0.2247] [0.2232, 0.2235]
p6 0.3118 0.3103 0.3117 [0.3069, 0.3167] [0.3113, 0.3123]

Table 9. Variance of the steady-state vector in M/H2 /1/7 queue, together with its 95% confi-
dence interval.
Vðpi Þ 104 PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 4.7367 4.5451 4.7370 [4.43549, 5.1917] [4.6955, 4.7785]
p1 6.7973 6.5517 6.8120 [6.2495, 7.4503] [6.7382, 7.8574]
p2 7.1124 7.1342 7.1237 [6.5392, 7.7958] [7.0506, 7.1753]
p3 4.7642 4.5830 4.7726 [4.3802, 5.2218] [4.7228, 4.8063]
p4 0.7665 0.8177 0.7653 [0.7047, 0.8401] [0.7598, 0.7733]
p5 5.4708 5.2979 5.4699 [5.0299, 5.9964] [5.4232, 5.5191]
p6 64.0300 62.3300 64.0396 [58.8790, 70.1926] [63.4835, 64.6063]

Table 10. Skewness coefficient of the steady-state vector in M/H2 /1/7 queue, together with its 95%
confidence interval.
Skewðpi Þ PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 0.3481 0.2807 0.3531 [0.1962, 0.4999] [0.3329, 0.3632]
p1 0.1478 0.2379 0.1474 ½0:0040; 0:2996 [0.1326, 0.1630]
p2 0.0881 0.0925 0.0896 ½0:2399; 0:0637 ½0:1033;0:0729
p3 0.4210 0.4300 0.4230 ½0:5728;0:2691 ½0:4361;0:4058
p4 1.2389 1.1806 1.2381 ½1:3907;1:0871 ½1:2541;1:2237
p5 0.3119 0.3561 0.3101 ½0:4637;0:1601 ½0:3271;0:2967
p6 0.2751 0.3242 0.2747 [0.1233, 0.4270] [0.2600, 0.2903]

Table 11. Kurtosis coefficient of the steady-state vector in M/H2 /1/7 queue, together with its 95%
confidence interval.
Kurtðpi Þ PC MC (NMC ¼ 103 ) MC (NMC ¼ 105 ) 95% CI (NMC ¼ 103 ) 95% CI (NMC ¼ 105 )
p0 1.8643 1.9138 1.8624 [1.5606, 2.1679] [1.8339, 1.8947]
p1 1.7524 1.7200 1.7568 [1.4487, 2.0560] [1.7220, 1.7827]
p2 1.7297 1.6881 1.7266 [1.4261, 2.0333] [1.6993, 1.7601]
p3 1.8885 1.8168 1.8935 [1.5849, 2.1922] [1.8581, 1.9189]
p4 3.3357 3.7223 3.3412 [3.0320, 3.6393] [3.3053, 3.3660]
p5 1.7598 1.7821 1.7586 [1.4562, 2.0635] [1.7295, 1.7902]
p6 1.8392 1.8446 1.8372 [1.5356, 2.1429] [1.8089, 1.8696]

the one in service), and customers attempting to enter the queue when N customers are
already present are lost. A single server with exponential distributed service times with
rate l is considered. We assume that the server can break down only if the system is
not empty. The lengths of breakdowns are identically distributed and follow an expo-
nential distribution with rate a. The duration of reparation of the server is assumed to
follow an exponential law with parameter b. After the server being repaired, it switches
to working state and continues to provide service until the system becomes empty. The
arrival flow customers, service times, breakdown times and repair times are assumed to
be mutually independent input random variables. The service discipline is FCFS.
16 K. BACHI ET AL.

Figure 2. Marginal probability density fpi of the steady-state vector in M/H2 /1/7 queue.

A typical state of this system may be denoted by XðtÞ ¼ fðNðtÞ; YðtÞÞ; t  0g, where
N(t) is the number of customers in the queue at time t, and YðtÞ is a random variable
representing the server state at time t. If at time t, the server is experiencing a break-
down, then Y(t) ¼ 1; otherwise, the server is in the working state and Y(t) ¼ 0. The
stochastic process X(t) is a continuous time Markov chain whose state space
S ¼ fðn; 0Þ : n ¼ 0; :::; Ng [ fðn; 1Þ : n ¼ 1; :::; Ng. The infinitesimal generator matrix Q
of the continuous time Markov chain X(t) has the following block-tridiagonal form:
0 1
A0 B0 0 0 0 0 0 0 0 :::
B D1 B1 F1 0 0 0 0 0 0 ::: C
B C
B 0 D2 B2 F2 0 0 0 0 0 ::: C
B C
B 0 0 D3 B3 F3 0 0 0 0 ::: C
B C
B .. .. .. .. .. .. .. .. .. .. C
Q¼B B . . . . . . . . . . C;
C (4.4)
B 0 0 0 ::: Ds Ms Fs 0 0 C
::: C
B
B 0 0 0 0 ::: D M F 0 ::: C
B sþ1 sþ1 sþ1 C
B .. .. .. .. .. .. .. .. .. .. C
@ . . . . . . . . . .A
0 0 0 0 0 0 0 ::: DN EN
where
   
l ðk þ l þ aÞ a
B0 ¼ k; A0 ¼ k; D1 ¼ ; Bi ¼ ; i ¼ 1; :::q1;
0 0 k
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
17

Figure 3. Markov transition diagram of the M/M/1/N queue with server breakdowns and threshold-
based recovery policy.
! !
k 0 l 0
Fj ¼ ; j ¼ 1; :::; N1; Dk ¼ ; k ¼ 2; :::N
0 k 0 0
! !
ðk þ a þ lÞ a ða þ lÞ a
Ms ¼ ; s ¼ q; :::N1; EN ¼ :
b ðb þ kÞ b b

The Markov transition diagram of the process is given in Figure 3.

Denote by pm;n ; m ¼ 0; 1; n  0 the stationary distribution of the Markov chain
X(t). We have p0 ¼ p0;0 , and pn ¼ ðp0;n ; p1;n Þ; n ¼ 1; :::; N. The steady-state analysis of
the Markov chain X(t) can be performed via a recursive scheme (Yang, Chiang, and
Tsou 2013), or by using matrix analytic methods (Neuts (1981)). In the sequel, we fol-
low a different method, and we will provide a functional approach based on the PC
expansion for obtaining the stationary distribution pm;n in terms of some parameters,
which are obtained under the epistemic uncertainties.

4.2.1. Epistemic uncertainty in queueing system

In this model, we consider the epistemic uncertainty in all the input parameters of the
queueing system, k, l, a, b. However, now pm;n ðhÞ is assumed to be the random output
response of the model of the queueing system, where h ¼ ½h1 ¼ k; h2 ¼ l; h3 ¼ a; h4 ¼ b
is the vector that describes the model parameters. We consider an epistemic uncertainty in
the input parameter h of known probability density function. Due to the uncertainty in
input parameters, the responses pm;n ðhðxÞÞ are considered as random variables. The
uncertainty analysis of the functional pm;n may be quantified by the computation of its
mean and variance and moments of higher order.
Furthermore, we proceed to a sensitivity analysis in order to determine the most sen-
sitive parameters with respect to the performance measures. A convenient way to com-
pute those quantities is to approximate the functional pm;n ðhðxÞÞ under the form of a
PC expansion. The sensitivity analysis introduced in Sec. 3 help us to point out which
parameters generate a higher sensitivity with respect to the stationary distribution. We
consider the following perturbation models for all the uncertain parameters:
hj ðxÞ ¼ h^j þ rj ej ðxÞ; j ¼ 1; :::; 4;
where hj ðxÞ is random variable, h^j ; rj , are respectively, estimated mean and the stand-
ard deviation of the random variable hj obtained by statistical methods. ej ðxÞ is a ran-
dom noise imposed to hi which models epistemic uncertainty. Here, ej ðxÞ is assumed to
18 K. BACHI ET AL.

Table 12. First order Sobol’ indices for stationary distribution vector in M/M/1/5 queue with server
breakdowns and threshold-based recovery policy.
pm;n nShi Sa Sl Sb Sk
p0;0 0.3799 0.0366 0.0046 0.5789
p0;1 0.7942 0.0047 0.0096 0.1913
p1;1 0.0144 0.0038 0.0078 0.9740
p0;2 0.0632 0.0160 0.0080 0.9126
p1;2 0.1650 0.0003 0.0109 0.8236
p0;3 0.0445 0.0473 0.0073 0.9008
p1;3 0.7812 0.0151 0.0418 0.1619
p0;4 0.0502 0.0264 0.0006 0.9229
p1;4 0.3120 0.0124 0.0411 0.6345
p0;5 0.0528 0.0198 0.0044 0.9229
p1;5 0.1512 0.0085 0.0396 0.8005

follow a standard normal distribution. Let h ¼ ½h1 ðxÞ ¼ k; h2 ðxÞ ¼ l; h3 ðxÞ ¼

a; h4 ðxÞ ¼ b; pm;n ðhðxÞÞ and e ¼ ½1 ; :::; 4  be random fields.
The approximation of the measure pm;n ðhðxÞÞ in the form of multivariate PC expan-
sions is given by:
p
X
P4
pm;n ðhÞ yi Wi ðhÞ;
i¼0

where Wi ðhÞ is the multivariate orthonormal Hermite polynomial and yi are the coeffi-
cients of the approximating series expansion, (2.13).

4.2.2. Numerical application

As an illustrative example, the queueing system described in Sec. 4.2 is studied. The
capacity is set to N ¼ 5 and the threshold is equal to q ¼ 3. In order to study the sensi-
tivity of the queue, we assume that the exact value of all the input parameters, a, b, k,
l, are not precisely known and we propose the following perturbation on these parame-
ters:
a ¼ ^a þ ra e1 ; e1 N ð0; 1Þ ^ þ rb e2 ;
b¼b e2 N ð0; 1Þ;
k ¼ ^k þ rk e3 ; e3 N ð0; 1Þ ^ þ rl e4 ;
l¼l e4 N ð0; 1Þ:
so a, b, k, l, are Gaussian random variables. The mean and the standard deviation of
these random variables are set respectively to l ^ ¼ 7:3; ^k ¼ 2; ^a ¼ 3; b
^ ¼ 4; ra ¼
0:04; rb ¼ 0:02; rk ¼ 0:04 and rl ¼ 0:02.
Table 12 gives the first order Sobol’ indices for all the components of the stationary
distribution vector. From this table, we see that for all but p0;1 and p1;3 , the most influ-
ential input random parameter is k. Table 13, which gives the second order Sobol’ indi-
ces, shows that the combined effect of two input random parameter is negligible in
comparison with the effect of a single parameter. The same is also true for higher order
indices since the total Sobol’ indices are almost equal to first order ones, as shown by
Table 14. Note that for the components p0;1 and p1;3 ; it is the parameter a which is
most influential.
When considering p0;1 and p1;3 , we can set a as being the only random parameter
and when considering all the other components of the stationary distribution, we can
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
19

Table 13. Second order Sobol’ indices for stationary distribution vector in M/M/1/5 queue with ser-
ver breakdowns and threshold-based recovery policy.
pm;n nShi ;hj Sa;l 108 Sa;b 108 Sa;k 106 Sl;b 109 Sl;k 108 Sb;k 107
p0;0 0.6972 1.3636 0.2051 1.8411 4.9130 1.1427
p0;1 0.0144 2.8434 0.0127 0.4145 0.0016 0.6002
p1;1 0.0361 0.0175 0.5309 0.3338 0.0569 1.9217
p0;2 0.0516 0.7513 0.0183 0.0181 0.0249 0.8644
p1;2 0.0263 0.0365 0.0143 0.0109 0.0013 7.3456
p0;3 0.8840 0.0198 0.3924 0.0283 0.4031 0.8201
p1;3 9.8593 1.6277 0.3311 0.1568 0.0997 0.0147
p0;4 0.1210 7.1050 4.0803 4.9641 0.0012 8.0038
p1;4 0.0137 0.0368 0.4395 0.9306 0.6662 8.6512
p0;5 0.0113 0.2362 0.5077 0.6212 0.0029 0.1801
p1;5 0.0250 0.0959 0.0187 0.0433 0.0010 0.0427

Table 14. Total order Sobol’ indices for steady-state vector in M/M/1/5 queue with server break-
downs and threshold-based recovery policy.
pm;n nSThi STa STl STb STk
p0;0 0.3799 0.0366 0.0046 0.5789
p0;1 0.7943 0.0047 0.0096 0.1914
p1;1 0.0144 0.0038 0.0078 0.9741
p0;2 0.0634 0.0160 0.0080 0.9128
p1;2 0.1652 0.0002 0.0109 0.8237
p0;3 0.0445 0.0473 0.0074 0.9009
p1;3 0.7812 0.0151 0.0419 0.1619
p0;4 0.0502 0.0264 0.0005 0.9229
p1;4 0.3121 0.0124 0.0411 0.6345
p0;5 0.0528 0.0198 0.0044 0.9230
p1;5 0.1514 0.0085 0.0396 0.8007

Table 15. Expected value of the stationary distribution in M/M/1/5 queue with server breakdowns
and threshold-based recovery policy.
Eðpm;n Þ PC (1 r.v.) PC (4 r.v.)
p0;0 0.2639 0.2639
p0;1 0.0723 0.0723
p1;1 0.1084 0.1084
p0;2 0.0495 0.0495
p1;2 0.1827 0.1827
p0;3 0.0636 0.0636
p1;3 0.0927 0.0927
p0;4 0.0428 0.0428
p1;4 0.0523 0.0523
p0;5 0.0261 0.0261
p1;5 0.0457 0.0457

Table 16. Variance of the stationary distribution in M/M/1/5 queue with server breakdowns and
threshold-based recovery policy.
Vðpm;n Þ 106 PC (1 r.v.) PC (4 r.v.)
p0;0 14.711 25.338
p0;1 0.7178 0.8879
p1;1 2.4845 2.5455
p0;2 0.4445 0.4905
p1;2 4.2497 5.0989
p0;3 0.7925 0.8793
p1;3 0.6677 0.8450
p0;4 1.4357 1.5550
p1;4 0.9757 1.5432
p0;5 1.2546 1.3592
p1;5 3.4476 4.3111
20 K. BACHI ET AL.

Table 17. Skewness coefficient of the steady-state vector in M/M/1/5 queue with server breakdowns
and threshold-based recovery policy.
Skewðpm;n Þ PC (1 r.v.) PC (4 r.v.)
p0;0 0.0350 0.0395
p0;1 0.0443 0.0415
p1;1 0.0350 0.0371
p0;2 0.0722 0.0411
p1;2 0.0200 0.0362
p0;3 0.0536 0.0523
p1;3 0.0316 0.0556
p0;4 0.0011 0.0042
p1;4 0.0622 0.0116
p0;5 0.0488 0.0576
p1;5 0.0374 0.0668

Table 18. Kurtosis coefficient of the steady-state vector in M/M/1/5 queue with server breakdowns
and threshold-based recovery policy.
Kurtðpm;n Þ PC (1 r.v.) PC (4 r.v.)
p0;0 3.0016 3.0020
p0;1 3.0037 3.0046
p1;1 3.0016 3.0024
p0;2 3.0071 3.0064
p1;2 2.9999 3.0023
p0;3 3.0041 3.0035
p1;3 3.0017 3.0060
p0;4 2.9983 2.9979
p1;4 3.0044 2.9988
p0;5 3.0013 3.0028
p1;5 2.9995 3.0052

consider the parameter k as sole input random variable. This is what we do as a second
numerical experiment. Then, Tables 15–18 compare the statistical moments of the sta-
tionary distribution vector when the four input parameters are random (first column)
to the case where only the most influential parameter is random (second column). As
expected, the two columns of the tables match rather well.

5. Conclusion
In this paper, we have developed a numerical approach based on polynomial chaos
expansion, to study the sensitivity and the propagation of the epistemic uncertainty in
queueing models that occurs with unreliable servers. To illustrate the applicability of the
proposed approach, two models of queueing systems have been investigated. In the first
model (M/G/1/N queue with breakdowns and repairs), the epistemic uncertainty only
affects one input parameter, whereas in the second model (M/M/1/N queue with server
breakdowns and threshold-based recovery policy), it affects four parameters. In the lat-
ter case, a sensitivity analysis using Sobol’ indices is performed. When considering the
stationary distribution as output quantity of interest, it was shown that the parameter k
of the Poisson law modeling the customers arrival at the queue, was the most influential
factor. This finding was then confirmed by considering k as the only random input par-
ameter and by setting the three remaining ones to their average values. In that case, it
was shown that the statistical moments of the output measure were relatively insensitive
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
21

to the other parameters. Finally, comparisons with Monte-Carlo simulations showed the
good convergence properties of the chaos expansion. Using dependent input r.v. would
lead to more realistic model, however this would require a radically different approach
to propagate randomness into the system. This is the topic of ongoing research.

Acknowledgments
The authors would like to thank the anonymous reviewers for their valuable remarks and con-
structive comments that contributed to improve the final version of the paper.

ORCID
Katia Bachi http://orcid.org/0000-0002-0344-3718
C
edric Chauviere http://orcid.org/0000-0002-2870-5083
Hacene Djellout http://orcid.org/0000-0002-5304-4926
Karim Abbas http://orcid.org/0000-0001-6118-2889

References
Abate, J., G. L. Choudhury, and W. Whitt. 2000. An introduction to numerical transform inver-
sion and its application to probability models. In: Computational probability, ed. W. K.
Grassmann, Kluwer Academic Publishers.
Abbas, K., J. Berkhout, and B. Heidergott. 2016. A critical account of perturbation analysis of
markov chains. Markov Processes And Related Fields. 522 (2):227–65.
Aleti, A., C. Trubiani, A. van Hoorn, and P. Jamshidi. 2018. An efficient method for uncertainty
propagation in robust software performance estimation. Journal of System Software 138:
222–35. doi:10.1016/j.jss.2018.01.010.
Avi-Itzhak, B., and P. Naor. 1963. Some queuing problems with the service station subject to
breakdown. Operations Research 11 (3):303–20. doi:10.1287/opre.11.3.303.
Baccelli, F., and T. Znati. 1981. Queueing systems with breakdowns in data base modeling. In
Proceedings of performance 81 (8 th IFIP international symposium on comp. Perf. Model. North
Holland: Amsterdam.
Cacuci, D. G. 2003. Sensitivity and uncertainty analysis. Vol. I. Boca Raton, FL: Chapman & Hall/
CRC.
Cao, J. H., and K. Cheng. 1982. Analysis of an M=G=1 queueing system with repairable service
station. Acta Mathematicae Applicatae Sinica 5 (2):113–27.
Chauviere, C., J. S. Hesthaven, and L. Lurati. 2006. Computational modeling of uncertainty in
time-domain electromagnetics. SIAM Journal on Scientific Computing 28 (2):751–75. doi:
10.1137/040621673.
Dhople, S. V., Y. C. Chen, and A. D. Dom
ınguez-Garc
ıa. 2013. A set-theoretic method for para-
metric uncertainty analysis in Markov reliability and reward models. IEEE Transactions on
Reliability 62 (3):658–69. doi:10.1109/TR.2013.2270421.
Dhople, S. V., and A. D. Dom
ınguez-Garc
ıa. 2012. A parametric uncertainty analysis method for
Markov reliability and reward models. IEEE Transactions on Reliability 61 (3):634–48. doi:
10.1109/TR.2012.2208299.
Gautschi, W. 2004. Orthogonal polynomials: computation and approximation. Numerical
Mathematics and Scientific Computation. Oxford University Press, New York. Oxford Science
Publications.
Ghanem, R. G., and P. D. Spanos. 1991. Stochastic finite elements: a spectral approach. New York:
Springer-Verlag.
22 K. BACHI ET AL.

Gottlieb, D., and S. A. Orszag. 1977. Numerical analysis of spectral methods: theory and applica-
tions. CBMS-NSF Regional Conference Series in Applied Mathematics, No. 26. Philadelphia,
PA Society for Industrial and Applied Mathematics.
Granger, M., and M. Henrion. 1993. Uncertainty: A guide to dealing with uncertainty in quantita-
tive risk and policy analysis. pp. 332. Cambridge, MA: Cambridge University Press.
Helton, J. C. 1993. Uncertainty and sensitivity analysis techniques for use in performance assess-
ment for radioactive waste disposal. Reliability Engineering and System Safety 42 (2–3):327–67.
doi:10.1016/0951-8320(93)90097-I.
Helton, J. C. 1994. Treatment of uncertainty in performance assessments for complex systems.
Risk Analysis 14 (4):483–511. doi:10.1111/j.1539-6924.1994.tb00266.x.
Helton, J. C. 1997. Uncertainty and sensitivity analysis in the presence of stochastic and subject-
ive uncertainty. Journal of Statistical Computation and Simulation. 57 (1-4):3–76. doi:10.1080/
00949659708811803.
Helton, J. C., J. D. Johnson, C. J. Sallaberry, and C. B. Storlie. 2006. Survey of sampling-based
methods for uncertainty and sensitivity analysis. Reliability Engineering and System Safety 91
(10-11):1175–209. doi:10.1016/j.ress.2005.11.017.
Homma, T., and A. Saltelli. 1996. Importance measures in global sensitivity analysis of nonlinear
models. Reliability Engineering and System Safety 52 (1):1–17. doi:10.1016/0951-8320(96)00002-6.
Iman, R. L., and J. C. Helton. 1988. An investigation of uncertainty and sensitivity analysis tech-
niques for computer models. Risk Analysis 8 (1):71–90. doi:10.1111/j.1539-6924.1988.tb01155.x.
Ionescu-Bujor, M., and D. G. Cacuci. 2004. A comparative review of sensitivity and uncertainty
analysis of large-scale system, I: Deterministic methods. Nuclear Science and Engineering 147
(3):189–203. doi:10.13182/NSE03-105CR.
Li, W., D. Shi, and X. Chao. 1997. Reliability analysis of M=G=1 queueing systems with server
breakdowns and vacations. Journal of Applied Probability 34 (2):546–55. doi:10.2307/3215393.
Limbourg, P. 2008. Dependability modelling under uncertainty: An imprecise probabilistic
approach. Studies in Computational Intelligence, 148, Berlin, Heidelberg: Springer.
Mattia, P., C. M. Sergio, and G. M. Dimitrov. 2007. Comparative analysis of uncertainty propaga-
tion methods for robust engineering design. In: Proceedings of the international conference on
engineering design, ICED’07. 28–31 August, Cit
e Des Sciences et de l’industrie, Paris.
Mishra, K., and K. S. Trivedi. 2011. An unobtrusive method for uncertainty propagation in sto-
chastic dependability models. IJRQP 3 (1):49–65.
Mitrany, I., and B. Avi-Itzhak. 1968. A Many-Server queue with service interruptions. Operation
Research 16 (3):628–38. doi:10.1287/opre.16.3.628.
Moore, R. 1979. Methods and applications of interval analysis. Philadelphia, PA: SIAM Studies in
Applied Mathematics.
Neuts, M. F. 1981. Matrix-geometric solutions in stochastic models, Volume 2 of johns hopkins ser-
ies in the mathematical sciences. An algorithmic approach. Baltimore, MD: Johns Hopkins
University Press.
Neuts, M. F., and D. M. Lucantoni. 1979. A markovian queue with N servers subject to break-
downs and repairs. Management Science. 25 (9):849–61. 1980. doi:10.1287/mnsc.25.9.849.
Ouazine, S., and K. Abbas. 2016. A functional approximation for retrial queues with two way com-
munication. Annals of Operations Research 247 (1):211–227. doi:10.1007/s10479-015-2083-2.
Pettersson, M. P., G. Iaccarino, and J. Nordstr€ om. 2015. Polynomial chaos methods for hyperbolic
partial differential equations. Mathematical engineering. Numerical techniques for fluid dynam-
ics problems in the presence of uncertainties. Cham: Springer.
Rahman, S. 2009. Extended polynomial dimensional decomposition for arbitrary probability dis-
tributions. Journal of Engineering Mechanics 135 (12):1439–1451. doi:10.1061/(ASCE)EM.1943-
7889.0000047.
Rocco, C. M., and W. Klindt. 1998. Distribution systems reliability uncertainty evaluation using
an interval arithmetic approach. In: Proceedings of the second IEEE international CARACAS
conference on devices, circuits and systems., Isla de Margarita, Venezuela.
Ronen, Y. 1988. Uncertainty analysis. Boca Raton, FL: CRC Press.
 COMMUNICATIONS IN STATISTICS - SIMULATION AND COMPUTATIONV
R
23

Sepahvand, K., S. Marburg, and H. J. Hardtke. 2010. Uncertainty quantification in stochastic sys-
tems using polynomial chaos expansion. International Journal of Applied Mechanics 02 (02):
305–53. doi:10.1142/S1758825110000524.
Shooman, M. 1990. Probabilistic reliability: An engineering approach. Second edition. Malabar,
FL: R. Krieger Pub. Co.
Shortle, J. F., M. J. Fischer, and P. H. Brill. 2007. Waiting-time distribution of M/DN/1 queues
through numerical laplace inversion. INFORMS Journal of Computing 19 (1):112–20. doi:
10.1287/ijoc.1050.0148.
Shuxing, Y., F. Xiong, and F. Wang. 2017. Polynomial chaos expansion for probabilistic uncer-
tainty propagation. Uncertainty Quantification and Model Calibration, Jan Peter Hessling,
IntechOpen. doi:10.5772/intechopen.68484. Available from: https://www.intechopen.com/books/
uncertainty-quantification-and-model-calibration/polynomial-chaos-expansion-for-probabilistic-
uncertainty-propagation.
Sobol, I. M. 1993. Sensitivity estimates for nonlinear mathematical models. Mathematcis
Modeling Computing Experiment 1 (4):407–14.
Sudret, M. 2008. Global sensitivity analysis using polynomial chaos expansion. Reliability
Engineering and System Safety 93 (7):964–79. doi:10.1016/j.ress.2007.04.002.
Takhedmit, B., and K. Abbas. 2017. A parametric uncertainty analysis method for queues with
vacations. The Journal of Computational and Applied Mathematics. 312:143–55. doi:10.1016/
j.cam.2016.02.031.
Thiruvengadam, K. 1963. Queuing with breakdowns. Operation Research 11 (1):62–71. doi:
10.1287/opre.11.1.62.
Wang, J., J. Cao, and Q. Li. 2001. Reliability analysis of the retrial queue with server breakdowns
and repairs. Queueing Systems 38 (4):363–80. doi:10.1023/A:1010918926884.
Winkler, R. L. 1996. Uncertainty in probabilistic risk assessment. Reliability Engineering System
Safety 54 (2/3):127–32. doi:10.1016/S0951-8320(96)00070-1.
Winkler, R. L. 2004. A comparative review of sensitivity and uncertainty analysis of Large-Scale
Systems - II: Statistical methods. Nuclear Science Engineering 147 (3):204–17.
Xiu, D., and G. E. Karniadakis. 2002. The Wiener-Askey polynomial chaos for stochastic differen-
tial equations. SIAM Journal on Scientific Computing 24 (2):619–44. doi:10.1137/
S1064827501387826.
Yang, D. H., Y. C. Chiang, and C. Tsou. 2013. Cost analysis of a finite capacity queue with server
breakdowns and threshold-based recovery policy. Journal of Manufacturing System 32 (1):
174–9. doi:10.1016/j.jmsy.2012.06.002.