key: cord-0952684-3ly0lznj
authors: Ameur, Lounes; Bachioua, Lahcene
title: Sensitivity analysis of queueing models based on polynomial chaos approach
date: 2021-09-27
journal: Arab
DOI: 10.1007/s40065-021-00344-y
sha: f85ca9af212afd02d26c1dbf39a9dea835b5077e
doc_id: 952684
cord_uid: 3ly0lznj

Queueing systems are modeled by equations which depend on a large number of input parameters. In practice, significant uncertainty is associated with estimates of these parameters, and this uncertainty must be considered in the analysis of the model. The objective of this paper is to propose a sensitivity analysis approach for a queueing model, presenting parameters that follow a Gaussian distribution. The approach consists in decomposing the output of the model (stationary distribution of the model) into a polynomial chaos. The sensitivity indices, allowing to quantify the contribution of each parameter to the variance of the output, are obtained directly from the coefficients of decomposition. The proposed approach is then applied to M/G/1/N queueing model. The most influential parameters are highlighted. Finally several numerical and data examples are sketched out to illustrate the accuracy of the proposed method and compare them with Monte Carlo simulation. The results of this work will be useful to practitioners in various fields of theoretical and applied sciences.

Queueing models have been greatly studied and successfully applied in fabrication and production systems, service systems, telecommunication, transport, and the service industry. Typically, a queueing model is a simplified representation of the real-life system under consideration. In practice, input parameters of the model are not known exactly and are subject to uncertainty, because they are determined by insufficient statistical data. For these reasons, sensitivity analysis of queueing systems (SAQS) based on polynomial chaos (PC) has been developed. SAQS allows to analyze a queueing model by studying the impact of the variability of the model parameters on the output variable. Determining the parameters responsible for this variability makes it possible to take the necessary measures to reduce the variance of the output when the uncertainty on the output is important. Determining the least influential parameters can simplify the model by fixing the parameters to a nominal value of their interval of variation without significant repercussions on the output. Different approaches to sensitivity analysis have been proposed in the literature of Saltelli [18] , Sobol [17] and Saltelli et al. [20] , [19] and [21] . Approaches can be local or global. In this paper, we are interested in global methods, very frequent, which are based on the analysis of the variance of the output (ANOVA, ANalysis Of VAriance), Saltelli et al. [20] and Sobol [17] , Jacques et al. [13] . The function of the model (here the stationary distribution π) is decomposed into a sum of functions of increasing dimensions Sobol [17] . This decomposition, called HDMR (High Dimensional Model Representation), allows to separate the effects of the different parameters, which are transmitted in the decomposition of the variance. To quantify the contribution of a parameter to the variance of the output, a sensitivity index is calculated. Sensitivity indices can sometimes be calculated formally, when the analytic form of the stationary distribution π of the model is known and relatively simple. in our queueing model the stationary distribution π is complex (nonlinear) and not analytically known. Unable to calculate these sensitivity indices, it is then necessary to estimate them Saltelli [18] and Gugole et al. [11] . Very often, the indices are estimated using Monte Carlo simulations. However, in some applications, this approach may require a high number of model evaluations, which makes the approach costly in computing time. Moreover, the measures are not always chosen freely because of the constraints of their implementation. In this case, the approach based on Monte Carlo simulations is no longer adapted. To overcome this problem, we will replace the queueing model by its analytical approximation based on the (PC) decomposition, which is less costly. The indices are then obtained directly from the expression of the coefficients of the polynomial decomposition. The (CP) decomposition of the output has recently been used in an original way for the sensitivity analysis of Blatman and Sudret [6] and Blatman [5] . Wiener [23] has shown that any random variable following a normal distribution, and which possesses a finite variance, may be represented by a decomposition into Hermite polynomials, thus forming an orthogonal basis. Xiu and Karniadakis [25] have extended this decomposition to other types of distributions, which are associated with specific families of polynomials.

In this paper, we propose a numerical approach under uncertainty for M/G/1/N queueing model based on polynomials chaos, whose input parameters follow normal distributions; otherwise, in the case where the input parameters follow an arbitrary distribution, the transformation of these parameters to normal distributions is possible; see (Hamza et al. [12] , Zhao and Ono [26] ).

In this regard, we approach the expression of the stationary distribution π, the expected, and variance of different performance measures. To illustrate the accuracy of the proposed methods, we compare the obtained results with those obtained from repeated Monte Carlo simulations.

The rest of this paper is organized as follows: The embedded Markov chain of the M/G/1/N model is presented in Sect. 2. In Sect. 3, the expression of sensitivity indices is given. Their estimation based on decomposition in CP is also detailed. In Sect. 4, the proposed approach is applied to M/G/1/N queue, and numerical results of the computation of the stationary distribution and some performance measures under epistemic uncertainty are discussed in Sect. 5. Finally, Sect. 6 concludes the paper.

Consider a computational model M whose input parameters are represented by a random vector X = (X 1 , X 2 , . . . , X d ), (the parameters X i , i = 1, 2, . . . , d are assumed to be independent real random variables) and the associated quantity of interest

the response Y is also uncertain, and its statistics are unknown and have to be estimated. Denote by ( , A, P) the probability space, where is the set of all possible outcomes, A is a σ -algebra over , and P is a function A → [0, 1] that gives a probability measure on A Consider an independent random vector X = (X 1 , X 2 , . . . , X d ) that describes input uncertainties. The probability law of X may be defined by the probability density function

where f i (x i ) is the marginal probability density of X i . Assuming that Y defined in (1) has a finite variance (which is a physically meaningful assumption when dealing with geotechnical systems), it belongs to the so-called Hilbert space of second-order random variables, which allows for the following representation [24] :

where ψ p (.) is the polynomial chaos of order p.

After some rearranging (see [10] ), this equation can be rewritten in a more convenient way as

with deterministic coefficients b j which are unknown, to be estimated, ψ j are multivariate polynomials that depend of X 1 , . . . , X d . It is given by

where the multi-indices (also called tuples) α ∈ N d which are ordered lists of integers

and φ α i (X i ) is a univariate polynomials depend on X i . The two-dimensional counterpart of Eq. (3) is rewritten, in a fully expanded form, as

The Eq. (7) can be recast in terms of ψ j (.) as follows:

Now, for example the term a 211 ψ 3 (X 2 , X 1 , X 1 ) of Eq. (7) is identified with the term b 7 ψ 7 of Eq. (8) . For a given degree p, the standard truncation scheme consists in selecting all polynomials, such that

The number of elements of a complete multivariate polynomial basis of degree p is

The set of complete multivariate orthogonal basis {ψ k } k≥0 is consistent with the density of X . This family can be obtained by applying the Gram-Schmidt orthogonalization procedure to the canonical family of monomials {1, x 2 , x 2 , . . .}. For standard distributions, the associated family of orthogonal polynomials is well known. If X i ∼ U (−1, 1) has a uniform distribution over [−1, 1], the resulting family is that of so-called Legendre polynomials. If X i ∼ N (0, 1) has a standard normal distribution with zero mean value and unit standard deviation, the resulting family is that of Hermite polynomials. The families associated with standard distributions are summarized in table (Xiu and Karniadakis [24] ). Remark An input variable X i with a different distribution F i (continuous, strictly increasing) may be transformed to a random variable X i (see [12] ) with distribution F belonging to one of the classical families given in Fig. 1 .

The series (4) is truncated by keeping terms up to a degree p

where P p d is given in (10) . Note that P p d increases polynomially with d and p. Thus, the number of terms in the Chao polynomial series, i.e., the number of coefficients to be computed, increases dramatically when d > 10, say. This complexity is referred to as curse of dimensionality, and there are advanced truncation schemes that allow one to bypass this problem (see [3, 14, 15] ). Now, let us give the Hermite polynomial chaos which is the very common example studied in the literature (the input random variables are N (0, 1)).

The original polynomial chaos, also termed as the homogeneous chaos, was first introduced by Wiener [22] . It employs the Hermite polynomials in terms of Gaussian random variables.

According to (3), the random response Y can be represented in the form

where H d (X i 1 , . . . , X i d ) denote the Hermite polynomials of order d in terms of the multi-dimensional independent standard Gaussian random variables X = (X i 1 , . . . , X i d ) with zero mean and unit variance. The above 

equation is the discrete version of the original Wiener polynomial chaos expansion, where the continuous integrals are replaced by summations. The general expression of the Hermite polynomials is given by

For example, the one-dimensional Hermite polynomials are

For notational convenience, Eq. (12) can be rewritten as

where there is a one-to-one correspondence between the functions H d (X i 1 , . . . , X i d ) and ψ j (X ), and also between the coefficients b j and a i 1 ...i r . In Eq. (12), the summation is carried out according to the order of the Hermite polynomials, while in Eq. (15) , it is simply a re-numbering with the polynomials of lower order counted first. For clarity, the two-dimensional expansion is shown here, both in the fully expanded form [see Eq. (12)]

and the simplified form [see Eq. (15)]

Example 2.1 The two-dimensional hermite polynomial chaos for different degree p is given by (Table 1) . [10] give the three-dimensional and four-dimensional hermite polynomial chaos for different degree p.

Several approaches have been developed to estimate the coefficients of the CP that can be classified as intrusive (Xiu and Karniadakis [25] ) or non-intrusive (Berveiller et al. [4] ). Among the non-intrusive approaches, we quote: the method of projection, the method of regression, and the method of the least squares which will be retained in the sequel, for its simplicity of implementation. We denote by B the vector of the coefficients b j , j = 0, . . . , d, to be determined. Let U be the vector containing the N samples of the output Y and Z the matrix containing the polynomials ψ j

k represents the ith sample of X k . The output vector U is given by U = Z B. The estimator of the vector B, denoted B, is given by the least-squares method

The estimated output U = ( Y (1) , . . . , Y (n) ) T is thus given by the following expression:

Consider a mathematical model, consisting of a set of random input variables, a deterministic function, and a set of random output variables(or responses). We write this model in the following form:

where

The function Y is the response of the model, and the set of input variables X = (X 1 , X 2 , . . . , X d ) groups all the entities considered random in the model. Quantities X i , i = 1, . . . , p are assumed to be independent random variables. The purpose is to study the sensitivity analysis of each component Y . An indicator of the influence of the parameter X i on the output Y denoted S i (see Sobol [17] ) is defined by

This index is called first-order sensitivity index of Sobol [17] . It quantifies the sensitivity of the output Y to the input variable X i , or the part of the variance of Y due to the variable X i . The sensitivity index S i is between 0 and 1; more the index is elevated (close to 1), more the variable X i is influential (will have importance).

The variance of the model with independent inputs can be decomposed as follows:

with

based on this decomposition (see Sobol [17] ), the sensitivity indices of second order are given by

These indices express the sensitivity of the variance of Y to the interaction of the variables X i and X j , that is, the sensitivity of Y to the variables X i and X j which is not taken into account in the effect of the variables alone. The higher order sensitivity indices are defined on the same principle. The total sensitivity index (see Gugole et al. [11] ), Cacuci et al. [7] which expresses the total sensitivity of the variance Y to a variable X i , that is, sensitivity to this variable in all its forms (sensitivity to the variable alone and sensitivity to the interactions of this variable with other variables), is defined as the sum of all the sensitivity indices relative to the variable X i

For example, for a model with three input variables, we have

Since the model function is not known analytically, it is necessary to estimate these indices. In the literature, different estimation methods have been proposed by Cukier et al. (1978) . The approach frequently used for its simplicity of implementation is based on simulations of Monte Carlo. However, the number of evaluations of the model can become high and costly in calculation time. To avoid this disadvantage, the model is replaced by its decomposition on (PC), which is an analytic representation in an orthogonal polynomial basis. The sensitivity indices are obtained directly from the PC decomposition coefficients.

According to (11) , each component Y of Y can be developed as follows:

After determining the coefficients of the CP and by reordering the terms of (27) according to the parameters on which they depend, we obtain the following decomposition:

with Γ k 1 ,...,k s the set of multi-indices j which corresponds to the polynomials dependent only on variables X k 1 ,...,k s . According to expression (28), the estimate of the first-order sensitivity index S i is given by

where b j are the coefficients estimated by (18) and the set Γ i corresponds to the polynomials ψ j depending only on X i .

In the same way, the estimate of the sensitivity index S ir due to the interaction between the variables X i and X r is given by

Moreover, the estimated total sensitivity index is written in this form (see [11] )

The CP forms an orthogonal basis, which guarantees the uniqueness of the decomposition and hence the uniqueness of the sensitivity indices (29)-(31).

In this section, we consider the PC expansion for propagating the uncertainty in performance measures of the M/GI/1 queue, and study the influence of each input parameter and their possible interactions onto the output measures.

We consider an M/G/1 queue with finite capacity N, where clients are served according to the discipline FIFO. The customers arrive according to a Poisson process with rate λ. The service times is a sequence of random variables distributed according to a general common distribution G(t) of mean 1/μ. In the sequel, we discuss an example motivating the practical features of this queueing model [2] . Storage Area Networks (SAN) have become one of the most widely used solutions for high-performance enterprise storage applications. However, the infrastructures of SANs are very complex and varied among different vendors. It is important to understand the behavior of these SANs in order of enhancing performance of the applications that run on them. The main idea of this paper is to understand the behavior of a server-SAN system that manages multiple applications stored in the same volume. The aim is modeling a typical workload of a small company which consists of one license of database system software and shares that same license for multiple databases. These multiple databases may be used potentially for diverse applications. The objective is then to see how the system would behave when managing queries from different types of databases at the same time. An M/G/1 Queueing Model has been used by [Emmanuel Arzuaga and David R. Kaeli] to analyze the system and they give a report set of experiments which are composed of two different database workloads based on the TPC-C and TPC-H benchmarks; they also show that the model can be used to anticipate the impact of growth in the arrival rate on the behavior of the system. This model has the potential of predicting the system behavior when adding more workloads to the mix.

Let Q n = Q(t + n ) be the number of customers in the system just after the departure of the nth customers; it is clear that Q n is an induced Markov chain with space state 0, 1, . . . , N − 1. The one-step transition matrix associated of the Markov chain Q n (see [9] ) is given by

where

A service-time distribution (G), having a finite ith moment m (i) , such that its coefficient of variation CV, can be approximated by a phase type distribution PH. Therefore, we use M/PH/1/N queue as approximation to M/GI/1/N queue 1. if the coefficient of variation of the service-time distribution CV = 0, the distribution is deterministic; 2. if the coefficient of variation of the service-time distribution CV = 1, the PH distribution is exponential; 3. if the coefficient of variation of the service-time distribution CV < 1, the PH distribution is hypoexponential; 4. if the coefficient of variation of the service-time distribution CV > 1, the PH distribution is hyperexponential.

To test the efficiency of the polynomial chaos decomposition and explore particularly the sensitivity of the input parameters on the considered performance measures. We limit our study with particular case when the service time is exponential and we provide a variety of numerical results for various performance measures and the results are exhibited via tables and graphs. All the calculations are done on Matlab Software Package.

In the next section, we will discuss a functional approach based on PC expansion for computing the epistemic uncertainty in stationary distribution π(λ, μ) = (π 1 (λ, μ) , . . . , π N (λ, μ)), due to epistemic uncertainties in λ and μ. In our numerical computation, we set X 1 = λ the arrival rate, and X 2 = μ the service rate, which are supposed random variables of normal distribution, with mean λ, μ, respectively. The output of the model represents the stationary distribution π(λ, μ). In the M/M/1 system, we consider the following situations: λ = μ = 2/3, σ 1 = σ 2 = 0.1 and N = 4, The aim is to study the influence of the λ and μ parameters on the variance of the response of the model π, and we apply the approach proposed in the Sect. 2.

Each component of the output of the model π is decomposed into PC (27)

The number of parameters d = 2, and the desired degree for the polynomials is fixed to p = 3. According to (10) , the number of coefficients b j , = 1, . . . , 4 is equal to 10. {ψ j } 0≤ j≤9 form an orthonormal Hermite polynomial basis:

The evolutions of the real output π and the estimated output U (Sect. 2.2) are shown in Fig. 2 , and the point cloud of the real output versus the estimated output is given in Fig. 3 .

The sensitivity indices of first order, S 1 and S 2 , corresponding, respectively, to the arrival rate, X 1 = λ, and the service rate, X 2 = μ, are estimated by (29). Similarly, the index of second order, S 12 , due to the interaction between λ and μ, is estimated by (30). Finally, the total indices, S T 1 and S T 2 , are estimated by (31). These indices are presented in Table 5 .

We observe that the parameter λ has the highest sensitivity index for each component of the stationary distribution. Then comes the parameter μ. This explains the importance of the contribution of the arrival rate on the performance of the system. Moreover, the total sensitivity indices are close to the first-order indices, since the contribution due to the interaction between λ and μ is low. This means that the contributions of λ and μ are additive.

In this section, we simplify the model by fixing the least influential parameters (μ = μ 0 ).

We will approximate the expectation and the variance of the stationary distribution as a function of the parameter λ which is obtained under epistemic uncertainty. λ can be written as follows:

where λ is the mean of arrival rate λ, σ its standard deviation and ε is the random noise inflicted on λ. This noise is called exogenous noise and follows the standard normal distribution. Figure 4 represents histogram and plot for the parameter λ given in (35) corresponding to the arrival rate. The function π , = 1, . . . , N where λ is given in (35), is approximated by a polynomial chaos [(see (27)]

The expectation and the variance of π for each = 1, . . . , N are given by

and

Since

(38) becomes as

In particular, we take N = 4, the desired degree for the polynomials is fixed to p = 3. According to (10) , the number of coefficients b j , = 1, . . . , N is equal to 4, and then, the polynomials ψ j (λ) are: ψ 0 (λ) = 1,

The formula (36) becomes

the expectation and the variance of each component of the stationary distribution are given by

and

where

The coefficients b j , j = 0, 1, 2, 3; = 1, 2, 3, 4 are estimated by (18) .

For our numerical experiment, we set λ = 2/3, μ 0 = 1, σ = 0.1, Furthermore, we validate the numerical results by Monte Carlo simulation. 

Let S be the number of customers in the system that is given as a function of the random variable λ [see (35)], and it is defined by

where π (λ) is given in (41). We define the average number of customers and the average waiting time of an arbitrary customers in the system respectively as follows (see [1] ):

and

where f λ is the probability density function of the random variable λ and λ is the effective arrival rate which is given by

Now, we calculate both performance measures L and W , and we will compare them with Monte Carlo simulation. The results obtained are given in Table 5 .

The comparison of the obtained results shows that the performance of the system calculated by polynomial chaos approximation are very close to that obtained by simulation. This explains that a third-order polynomial chaos yields a good approximation of the stationary distribution. 

Many countries are currently dealing with the COVID-19 epidemic and are searching for an exit strategy, such that life in society can return to normal. The load has increased tremendously in the hospitals where the resources are minimal. The queuing model is applied to identify the queue time of the patients in hospitals for the identification and confirmation of disease. M/M/1/N model is adopted for the study, to model the queue time for the patients in a single-server system. This model is defined as the M/M/1 model with the finite capacity, N . In this model, λ is defined as the arrival rate of patients (COVID patients only); μ is the mean service time of patients.

The patient approach to the hospitals for the treatment follows the Poisson probability distribution function at time t. The Poisson probability distribution function is defined for the λ (patients) per unit t (time), and all the patients will be served at a first-come, first-served (FCFS) basis. The time taken to treat the patients would be applied exponentially. The server time is assessed as μ (patients) per unit t (time). The service rate is not depending upon the arrival rate of patients, where the service rate is the time required for the treatment.

Data of COVID-19 patients were taken for the study in the period of 16th June 2021-8th August 2021 from the website (https://www.coronatracker.com).

The model has been designed specifically for the estimation of processing time of patients diagnosed with the novel Coronavirus. The number of cases of patients suffering from the Corona has increased exponentially, The single-server finite queuing model is designed to compute the time required for each patient. The MATLAB software is used to design the model. The patient data are used for testing the model using the arrival rate of patients (λ) as 3, and the service rate of patients (μ) per hour as 7. Figure 5 shows the increase in the daily cases of coronavirus in Algeria, and the number of recovered patients and casualties occurred due to COVID-19. It is visible from the graph that the rate of infections is similar to the rate of recovery of patients. However, the number of Corona positive cases is increasing at a rapid rate. Whereas, the fatalities were found low throughout the period.

The single-server queuing model is designed and results obtained based on the stationary distribution of the system are shown in Table 6 . The average number of patients in the system is found extensively high as 17.1131786315087, with the average time devoted by the patient in the system which is found as 6.348107106018424 min.

For a large and more general study, one can consider the multi-server queuing model (M/M/s/N) to identify the optimal number of server requirement with the less waiting time and minimum time spent by the patient in the system (see [16] ).

In this paper, a sensitivity analysis approach for an queueing model with parameters that follow a normal distribution has been proposed. The sensitivity indices are estimated by the method based on the PC decomposition of the output of the model. The advantage of this approach is that the number of evaluations of the model is few elevated and the indices are directly obtained from the coefficients of the PC. In addition, the proposed approach was applied to M/GI/1 queuing system, where we have demonstrated that the arrival rate is the most influential parameter on the system and we are fixed the service rate. We were able to estimate the values of each component of the stationary distribution while characterizing its expectation and its variance. Some important performance measures of the queueing model such as the mean system length and the mean waiting time of an arbitrary customer are discussed and evaluated under propagation of epistemic uncertainty in the model input parameter. Moreover, numerical results were presented and compared to the corresponding Monte Carlo simulations ones. The model adopted depends on two parameters; it would be interesting to study more complex models with similar approach.

Sensitivity analysis of the M/M/1 retrial queue with working vacations and vacation interruption

An M/G/1 queue model for multiple applications on storage area networks

Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis

Stochastic finite element: a non intrusive approach by regression. Revue Europeenne de

Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis. These de doctorat

Sparse polynomial chaos expansion and adaptive stochastic finite element using a regression approach

Sensitivity and Uncertainty Analysis

Nonlinear sensitivity analysis of multiparameter model systems

The basic M/G/1 queueing system

Stochastic Finite Elements: A Spectral Approach

Uncertainty quantification and sensitivity analysis of COVID-19 exit strategies in an individual-based transmission model

Analyse de sensibilité basée sur les polynomes du chaos pour des modéles de type boite noire

Sensitivity analysis in presence of model uncertainty and correlated inputs

The Art of Computer Programming: Generating All Combinations and Partitions

Fascicle 2: Generating All Tuples and Permutations (Art of Computer Programming)

A probabilistic model for the assessment of queuing time of coronavirus disease (COVID-19) patients using queuing model

Sensitivity estimates for nonlinear mathematical models

Making best use of model evaluations to compute sensitivity indices

Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models

Global Sensitivity Analysis: The Primer

The homogeneous chaos

The homogeneous chaos

The Wiener-Askey polynomial chaos for stochastic differential equations

Modeling uncertainty in flow simulations via generalized polynomial chaos

Third-moment standardization for structural reliability analysis

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

The authors are thankful to the anonymous referees and the editor for their valuable comments and suggestions which helped improve the quality of this paper.Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.