# Control of Traffic Intensity in Hyperexponential and Mixed Erlang Queueing Systems with a Method Based on SPRT

Authors: Müjgan Zobu; Vedat Sağlam
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The control of traffic intensity(ρ) is one of the important problems in the study of queueing systems. Rao et al. (1984) developed a method to detect changes in the traffic intensity in queueing systems of the M/G/1 and GI/M/c types based on the Sequential Probability Ratio Test (SPRT). In this paper, SPRT is theoretically investigated for two different phase-type queueing systems which consist of hyperexponential and mixed Erlang. Also, for testing H0:ρ=ρ0 against H1:ρ=ρ1, Operating Characteristic (OC) and Average Sample Number (ASN) functions are obtained with numerical methods using multipoint derivative equations according to different situations of α and β type errors. Afterward, numerical illustrations for each model are provided with Matlab programming.


---

## Body

## 1. Introduction

There are very few studies on the sequential detection of parameter changes of queueing systems in the literature. The theory of SPRT for a sequence of observations forming a finite Markov chain was given in [1]. After, based upon the theory of [1] was discussed statistical quality control and SPRT procedures for the control of traffic intensity in [2–6]. This method aims to detect changes in traffic intensity by observing only the number of customers in the system at successive departure epochs 𝒬n, which are embedded Markov points. Recently, a SPRT to regulate the traffic intensity based on the number of arrivals during the nth service periods for M/Ek/1 queue was proposed in [7] and an autoregressive process based on the number of customers at the departure point and its application to the queueing model were given in [8]. In this paper, the method based on the SPRT is first examined theoretically for the two different queueing phase types. Afterward, the observations obtained for different samples model are evaluated for each of the models by simulation through MATLAB 7.10.0 (R2010a) programming.


## 2. Queueing Process and SPRT

In many queueing applications, a performance characteristic of great importance and interest is the traffic intensityρ (the ratio of mean arrival rate λ to mean service rate μ). The purpose of this test is to determine as quickly as possible the changes in traffic intensity and to take the appropriate corrective actions. With this objective, a procedure has been developed for testing the hypothesis H0:ρ=ρ0 against H1:ρ=ρ1 using Wald’s SPRT for the systems M/G/1 and GI/M/c, in which queue length processes have imbedded Markov chains {𝒬n} in [6]. This procedure is explained as follows.Consider the single server queue where arrivals occur according to a Poisson process with rateλ per unit time. The service times of customers are independent identically distributed (i.i.d.) random variables with the distribution B(x). For this system, the queue lengths at service completion points form an imbedded Markov chain. 𝒬n will be the number of customers left behind by the nth departing customer. The capacity of the queueing system is restricted to N. Then the transition probability matrix of the imbedded Markov chain is given by
(1)P={Pij(ρ)}=012⋯N-101-∑0k-2kj11-∑0k-2kj21-∑0k-2kj·⋮··N-1[k01-∑0k-2kjk01-∑0k-2kj01-∑0k-2kj0⋮000k11-∑0k-2kjk11-∑0k-2kjk01-∑0k-2kj0⋮000k21-∑0k-2kjk21-∑0k-2kjk11-∑0k-2kjk0⋮000⋯1-∑0k-2kj⋯1-∑0k-2kj⋯1-∑0k-2kj⋯⋯⋯⋯1-∑0k-2kj1-∑0k-2kj1-∑0k-3kj···1-k0],
where kn=P{narrivals during a service period}=∫0∞e-λt((λt)j/j!)dB(t).Lett0<t1<t2<⋯<tn be the set of time points at which the 𝒬(t) process exhibits the Markov property; namely,
(2)P{𝒬(tn+1)=j∖𝒬(tn)=i,𝒬(tn-1)=in-1,…,𝒬(t0)=i0}=P{𝒬(tn+1)=j∖𝒬(tn)=i}.In theM/G/1 queue, tn is the time of service completion of the nth customer and, in the case of the GI/M/c queue, tn is the arrival epoch of the nth customer. In the case of the M/M/c/N queue, though the process is in Markovian continuous time, for purposes of control it is easier to examine the process at arrival points.For simplicity, we will use the notation𝒬n=𝒬(tn). Denoted by P={Pij(ρ)} the transition probability matrix of the Markov chain {𝒬n} defined over state space E={0,1,2,…,N}, where
(3)Pij(ρ)=P{𝒬n+1=j∖𝒬n=i;ρ}.It is assumed thatP is known except for the value of the parameter ρ. The problem is to test the hypothesis H0:ρ=ρ0 against H1:ρ=ρ1.Consider the sequence of observations𝒬0,𝒬1,𝒬2,…,𝒬n. The joint probability of observing this sequence under H0 and H1 is given by
(4)Pr{𝒬0,𝒬1,𝒬2,…,𝒬a;ρi}=P(𝒬0;ρi)∏j=1nP(𝒬j∖𝒬j-1;ρi),0000000000000000000000000000000000000000i=0,1.Then the likelihood ratio is(5)L=P(𝒬0;ρ1)∏j=1nP(𝒬j∖𝒬j-1;ρ1)P(𝒬0;ρ0)∏j=1nP(𝒬j∖𝒬j-1;ρ0).(6)Z0=lnP(𝒬0;ρ1)P(𝒬0;ρ0),(7)Zr=lnP(𝒬r∖𝒬r-1;ρ1)P(𝒬r∖𝒬r-1;ρ0),(r≥1).LetA=(1-β)/α and B=β/(1-α), where α and β are the probabilities of the errors of the first and second types. Then, Wald’s SPRT [9] to test H0:ρ=ρ0 against H1:ρ=ρ1 becomes as follows.Observe{𝒬i} (i=0,1,2,…) successively and at stage n≥1, (1)
acceptH0 if ∑0nZr≤lnB;(2)
acceptH1if ∑0nZr≥lnA;(3)
continue by observing𝒬n+1 if lnB<∑0nZr<lnA.If we assume𝒬0=i0 is specified and denote by nij the number of transitions i→j up to and including the nth transition, then the likelihood ratio given in (5) reduces to
(8)L=∏i,jPijnij(ρ1)∏i,jPijnij(ρ0),lnL=∑i,jnijlnPij(ρ1)Pij(ρ0).Then the SPRT for testing the hypothesisH0:ρ=ρ0 against H1:ρ=ρ1 will have its continuation region
(9)lnB<∑i,jnijlnPij(ρ1)Pij(ρ0)<lnA.
In the case of M/Ek/1, Ek/M/c, and M/M/c/N queues, we will show in the sequel that the logarithm of the likelihood ratio can be written in the form
(10)lnL=an+∑i,jnijcij,
where a, cij, and Pij are constants depending upon the parameters ρ1, ρ0, and the transition probabilities Pij. Thus (9) reduces to
(11)lnB-an<∑i,jnijcij<lnA-an.
From this it follows that the continuation region of the test is bounded by straight lines lnA-an ve lnB-an.


## 3. Operating Characteristic and Average Sample Number

Approximate formulas for the OC and ASN functions are given in [1]. In order to evaluate the OC and ASN functions, λ0(t,ρ) is determined as the largest real positive latent root of the matrix
(12)P(t)={pij(ρ)[pij(ρ1)pij(ρ0)]t},
its derivate at t=0 and the nonzero root t(ρ) of the equation λ0(t,ρ)=1.A numerical search technique was employed to find the root in [10]. For evaluating the derivate of λ0(t,ρ) at t=0, the following five-point numerical differentiation formula has been used:
(13)λ0′(t,ρ)≅112h{λ-2-8λ-1+8λ1-λ2},(14)λ0′′(t,ρ)≅112h2{-λ-2+16λ-1-30λ0+16h1-λ2}.When the stated space of{𝒬n} is finite, the OC function for the SPRT can be obtained as
(15)L(ρ)≅At0(ρ)-1At0(ρ)-Bt0(ρ)ift0(ρ)≠0≅lnAlnA-lnBift0(ρ)=0,
where t0(ρ) is the nonzero real root of the equation λ0(t,ρ)=1. The ASN can then be obtained as
(16)E(n;ρ)≅L(ρ)lnB+{1-L(ρ)}lnAλ0′(0)ifλ0′≠0(17)E(n;ρ)AaMm≅L(ρ){lnB}2+{1-L(ρ)}{lnA}2λ0′′ifλ0′=0.


## 4. SPRT for Phase-Type Distribution

The exponential distribution is very widely used in performance modelling. The reason, of course, is that mathematical tractability flows from the memoryless property of this distribution. But sometimes mathematical tractability is not sufficient to overcome the need for a model process in which the exponential distribution is simply not adequate. This leads us to explore ways in which we can develop more general distributions while maintaining some of the tractability of the exponential. This is precisely what phase-type distributions permit us to do [11].


### 4.1. Hyperexponential-2 Service Model: TheM/H2/1 Queue

Consider the configuration presented in Figure1, in which α1 is the probability that the upper phase is taken and α2=1-α1 the probability that the lower phase is taken. If such a distribution is used to model a service facility, then a customer entering a service will, with probability α1, receive the service that is exponentially distributed with parameter μ1 and then exit the server or else with probability α2 receive the service that is exponentially distributed with parameter μ2 and then exit the server. Once again, only one customer can be in the process of receiving service at any one time; that is, both phases cannot be active at the same time [11].Figure 1
Two exponential phases in parallel.Theorem 1.
The density function of the service time is given by(18)dB(x)dx=(α1μ1e-μ1x+α2μ2e-μ2x),(x≥0).
In this case, one has
(19)kn=α1(ρ1ρ1+1)n(1ρ1+1)+α2(ρ2ρ2+1)n(1ρ2+1),
where kn=P{narrivalsduringaserviceperiod}.Proof.
Consider(20)kn=∫0∞e-λt(λt)nn!(α1μ1e-μ1t+α2μ2e-μ2t)d(t)=∫0∞e-λt(λt)nn!α1μ1e-μ1td(t)+∫0∞e-λt(λt)nn!α2μ2e-μ2td(t)=α1μ1λnn!∫0∞e-λttne-μ1td(t)+α2μ2λnn!∫0∞e-λttne-μ2td(t)=α1μ1λnn!∫0∞e-t(λ+μ1)tnd(t)+α2μ2λnn!∫0∞e-t(λ+μ2)tnd(t)=α1μ1λnn!(1λ+μ1)n+1Γ(n+1)+α2μ2λnn!(1λ+μ2)n+1Γ(n+1)=α1μ1λnn!n!(1λ+μ1)n+1+α2μ2λnn!n!(1λ+μ2)n+1=α1(λλ+μ1)n(μ1λ+μ1)+α2(λλ+μ2)n(μ2λ+μ2)=α1(λ/μ1λ/μ1+μ1/μ1)n(μ1/μ1λ/μ1+μ1/μ1)+α2(λ/μ2λ/μ2+μ2/μ2)n(μ2/μ2λ/μ2+μ2/μ2)=α1(ρ1ρ1+1)n(1ρ1+1)+α2(ρ2ρ2+1)n(1ρ2+1).
Here
(21)Pij(ρ)=kj-i+1,i=1,2,…,N,000j=0,1,…,N-1,j≥i-1,P0j(ρ)=kj,j=0,1,…,N-1,P0N(ρ)=1-∑0N-1kn,PiN(ρ)=1-∑0N-ikn,i=1,2,…,N,PiN(ρ)=1-∑r=0N-iα1(ρ1ρ1+1)r(1ρ1+1)+α2(ρ2ρ2+1)r(1ρ2+1).
With these values for Pij(ρ), the logarithm of the likelihood ratio will be
(22)lnL=∑i,jnijlnPij(ρ1)Pij(ρ0),lnL=∑i,jnijln(+α2(ρ2(0)ρ2(0)+1)n(1ρ2(0)+1))-1(α1(ρ1(1)ρ1(1)+1)n(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)n(1ρ2(1)+1))×(α1(ρ1(0)ρ1(0)+1)n(1ρ1(0)+1)+α2(ρ2(0)ρ2(0)+1)n(1ρ2(0)+1))-1),
where
(23)cij=ln(+α2(ρ2(0)ρ2(0)+1)j-i+1(1ρ2(0)+1))-1(α1(ρ1(1)ρ1(1)+1)j-i+1(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)j-i+1(1ρ2(1)+1))×(α1(ρ1(0)ρ1(0)+1)j-i+1(1ρ1(0)+1)+α2(ρ2(0)ρ2(0)+1)j-i+1(1ρ2(0)+1))-1),000000i=1,2,…,N,j=0,1,…,N-1,j≥i-1,c0j=ln(α2(ρ2(0)ρ2(0)+1)j(1ρ2(0)+1))-1(α1(ρ1(1)ρ1(1)+1)j(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)j(1ρ2(1)+1))×(α1(ρ1(0)ρ1(0)+1)j(1ρ1(0)+1)+α2(ρ2(0)ρ2(0)+1)j(1ρ2(0)+1))-1),000000000000000000000j=0,1,…,N-1,ciN=ln{+α2(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1))∑r=0N-i(α1(ρ1(1)ρ1(1)+1)r(1ρ1(1)+1))-1(1-∑r=0N-i(α1(ρ1(1)ρ1(1)+1)r(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)r(1ρ2(1)+1)))×(1-∑r=0N-i(α1(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1)+α2(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1))∑r=0N-i(α1(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1))-1},i=1,2,…,N,c0N=c1N.


### 4.2. Mixed Erlang-k Service Model: The M/mixEk/1 Queue

This is illustrated in Figure2, where, with probability α, the top series of k-1 exponential phases is taken and, with probability 1-α, the bottom series of k phases is taken.Figure 2
Mixed Erlang representation.Theorem 2.
The density function of the service time is given by(24)dB(x)dx=αμ(μx)k-2e-μx(k-2)!+(1-α)μ(μx)k-1e-μx(k-1)!,(x≥0),
where k≥2 is an integer.
In this case, one has(25)kn=(ρρ+k)n(kρ+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ+k)),
where kn=P{narrivalsduringaserviceperiod}=∫0∞e-λt((λt)j/j!)dB(t).Proof.
Consider(26)kn=∫0∞e-λt(λt)nn!×(α(kμt)k-2(k-2)!e-ktxkμ+(1-α)(kμt)k-1(k-1)!e-kμtkμ)d(t)=∫0∞e-λt(λt)nn!α(kμt)k-2(k-2)!e-ktxkμd(t)+∫0∞e-λt(λt)nn!(1-α)(kμt)k-1(k-1)!e-ktxkμd(t)=αkμλn(kμ)k-2n!(k-2)!∫0∞e-λttne-kμttk-2d(t)+(1-α)kμλn(kμ)k-1n!(k-1)!∫0∞e-λttne-kμttk-1d(t)=αkμλnkμ(k-2)n!(k-2)!∫0∞e-t(λ+kμ)tn+k-2d(t)+αkμλnkμ(k-1)n!(k-1)!∫0∞e-t(λ+kμ)tn+k-1d(t)=αλn(kμ)k-1n!(k-2)!∫0∞e-t(λ+kμ)tn+k-2d(t)+αλn(kμ)kn!(k-1)!∫0∞e-t(λ+kμ)tn+k-1d(t)=kμλnkμ(k-1)n!(k-1)!∫0∞e-t(λ+kμ)tn+k-1d(t)⟶GamadistributionΓ(n+k)=αλn(kμ)k-1n!(k-2)!(1λ+kμ)n+k-1Γ(n+k-1)+(1-α)λn(kμ)kn!(k-1)!(1λ+kμ)n+kΓ(n+k)=(n+k-2)!αλn(kμ)k-1n!(k-2)!(1λ+kμ)n+k-1+(n+k-1)!(1-α)λn(kμ)kn!(k-1)!(1λ+kμ)n+k=(n+k-2n)(1λ+kμ)n(1λ+kμ)k-1α(kμ)k-1λn+(n+k-1n)(1λ+kμ)n(1λ+kμ)k(1-α)(kμ)kλn=α(n+k-2n)(λλ+kμ)n(kμλ+kμ)k-1+(1-α)(n+k-1n)(λλ+kμ)n(kμλ+kμ)k.
If the numerator and denominator are divided by μ,
(27)kn=α(n+k-2n)(λ/μ(λ+kμ)/μ)n((kμ)/μ(λ+kμ)/μ)k-1+(1-α)(n+k-1n)(λ/μ(λ+kμ)/μ)n((kμ)/μ(λ+kμ)/μ)k=α(n+k-2n)(ρρ+k)n(kρ+k)k-1+(1-α)(n+k-1n)(ρρ+k)n(kρ+k)k
and cobrackets are taken as follows:
(28)kn=(ρρ+k)n(kρ+k)k-1×(α(n+k-2)!n!(k-2)!+(1-α)(n+k-1)(n+k-2)!n!(k-1)(k-2)!(kρ+k))=(ρρ+k)n(kρ+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ+k)).
Here
(29)Pij(ρ)=kj-i+1,i=1,2,…,N,000j=0,1,…,N-1,j≥i-1,P0j(ρ)=kj,j=0,1,…,N-1,P0N(ρ)=1-∑0N-1kn,PiN(ρ)=1-∑0N-ikn,i=1,2,…,N,=1-∑r=0N-i(ρρ+k)r(kρ+k)k-1(n+k-2n)×(α+(1-α)(r+k-1)(k-1)(kρ+k)).
With these values for Pij(ρ), the logarithm of the likelihood ratio will be
(30)lnL=∑i,jnijlnPij(ρ1)Pij(ρ0)lnL=an+∑i,jnijcij,lnL=∑i,jnijln((α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1(n+k-2)!n!(k-2)!)-1((α+(1-α)(n+k-1)(k-1)(kρ1+k))(ρ1ρ1+k)n(kρ1+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ1+k)))×((ρ0ρ0+k)n(kρ0+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1(n+k-2)!n!(k-2)!)-1),lnL=∑i,jnijln((α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1)-1((α+(1-α)(n+k-1)(k-1)(kρ1+k))(ρ1ρ1+k)n(kρ1+k)k-1×(α+(1-α)(n+k-1)(k-1)(kρ1+k)))×((ρ0ρ0+k)n(kρ0+k)k-1×(α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1)-1),
where
(31)kn=n(k-1)ln(ρ0+kρ1+k)+∑i,jnijln(ρ1ρ0·ρ0+kρ1+k)n+∑i,jnij×ln(α+(1-α)((n+k-1)/(k-1))(k/(ρ1+k))α+(1-α)((n+k-1)/(k-1))(k/(ρ0+k))).
Here
(32)a=(k-1)ln(ρ0+kρ1+k),cij+bij=ln(ρ1ρ0·ρ0+kρ1+k)j-i+1+ln(α+(1-α)((j-i+k)/(k-1))(k/(ρ1+k))α+(1-α)((j-i+k)/(k-1))(k/(ρ0+k))),i=1,2,…,N,j=0,1,…,N-1,j≥i-1,c0j+b0j=ln(ρ1ρ0·ρ0+kρ1+k)j+ln(α+(1-α)((j+k-1)/(k-1))(k/(ρ1+k))α+(1-α)((j+k-1)/(k-1))(k/(ρ0+k))),j=0,1,…,N-1,ciN=ln{×(α+(1-α)(r+k-1)(k-1)(kρ0+k))(1-∑r=0N-i(ρ0ρ0+k)r(kρ0+k)k-1(r+k-2r))-1(1-∑r=0N-i(ρ1ρ1+k)r(kρ1+k)k-1(r+k-2r)×(α+(1-α)(r+k-1)(k-1)(kρ1+k))(1-∑r=0N-i(ρ1ρ1+k)r(kρ1+k)k-1(r+k-2r))×(1-∑r=0N-i(ρ0ρ0+k)r(kρ0+k)k-1(r+k-2r)×(α+(1-α)(r+k-1)(k-1)(kρ0+k))(1-∑r=0N-i(ρ0ρ0+k)r(kρ0+k)k-1(r+k-2r))-1}i=1,2,…,N,c0N=c1N.


## 4.1. Hyperexponential-2 Service Model: TheM/H2/1 Queue

Consider the configuration presented in Figure1, in which α1 is the probability that the upper phase is taken and α2=1-α1 the probability that the lower phase is taken. If such a distribution is used to model a service facility, then a customer entering a service will, with probability α1, receive the service that is exponentially distributed with parameter μ1 and then exit the server or else with probability α2 receive the service that is exponentially distributed with parameter μ2 and then exit the server. Once again, only one customer can be in the process of receiving service at any one time; that is, both phases cannot be active at the same time [11].Figure 1
Two exponential phases in parallel.Theorem 1.
The density function of the service time is given by(18)dB(x)dx=(α1μ1e-μ1x+α2μ2e-μ2x),(x≥0).
In this case, one has
(19)kn=α1(ρ1ρ1+1)n(1ρ1+1)+α2(ρ2ρ2+1)n(1ρ2+1),
where kn=P{narrivalsduringaserviceperiod}.Proof.
Consider(20)kn=∫0∞e-λt(λt)nn!(α1μ1e-μ1t+α2μ2e-μ2t)d(t)=∫0∞e-λt(λt)nn!α1μ1e-μ1td(t)+∫0∞e-λt(λt)nn!α2μ2e-μ2td(t)=α1μ1λnn!∫0∞e-λttne-μ1td(t)+α2μ2λnn!∫0∞e-λttne-μ2td(t)=α1μ1λnn!∫0∞e-t(λ+μ1)tnd(t)+α2μ2λnn!∫0∞e-t(λ+μ2)tnd(t)=α1μ1λnn!(1λ+μ1)n+1Γ(n+1)+α2μ2λnn!(1λ+μ2)n+1Γ(n+1)=α1μ1λnn!n!(1λ+μ1)n+1+α2μ2λnn!n!(1λ+μ2)n+1=α1(λλ+μ1)n(μ1λ+μ1)+α2(λλ+μ2)n(μ2λ+μ2)=α1(λ/μ1λ/μ1+μ1/μ1)n(μ1/μ1λ/μ1+μ1/μ1)+α2(λ/μ2λ/μ2+μ2/μ2)n(μ2/μ2λ/μ2+μ2/μ2)=α1(ρ1ρ1+1)n(1ρ1+1)+α2(ρ2ρ2+1)n(1ρ2+1).
Here
(21)Pij(ρ)=kj-i+1,i=1,2,…,N,000j=0,1,…,N-1,j≥i-1,P0j(ρ)=kj,j=0,1,…,N-1,P0N(ρ)=1-∑0N-1kn,PiN(ρ)=1-∑0N-ikn,i=1,2,…,N,PiN(ρ)=1-∑r=0N-iα1(ρ1ρ1+1)r(1ρ1+1)+α2(ρ2ρ2+1)r(1ρ2+1).
With these values for Pij(ρ), the logarithm of the likelihood ratio will be
(22)lnL=∑i,jnijlnPij(ρ1)Pij(ρ0),lnL=∑i,jnijln(+α2(ρ2(0)ρ2(0)+1)n(1ρ2(0)+1))-1(α1(ρ1(1)ρ1(1)+1)n(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)n(1ρ2(1)+1))×(α1(ρ1(0)ρ1(0)+1)n(1ρ1(0)+1)+α2(ρ2(0)ρ2(0)+1)n(1ρ2(0)+1))-1),
where
(23)cij=ln(+α2(ρ2(0)ρ2(0)+1)j-i+1(1ρ2(0)+1))-1(α1(ρ1(1)ρ1(1)+1)j-i+1(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)j-i+1(1ρ2(1)+1))×(α1(ρ1(0)ρ1(0)+1)j-i+1(1ρ1(0)+1)+α2(ρ2(0)ρ2(0)+1)j-i+1(1ρ2(0)+1))-1),000000i=1,2,…,N,j=0,1,…,N-1,j≥i-1,c0j=ln(α2(ρ2(0)ρ2(0)+1)j(1ρ2(0)+1))-1(α1(ρ1(1)ρ1(1)+1)j(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)j(1ρ2(1)+1))×(α1(ρ1(0)ρ1(0)+1)j(1ρ1(0)+1)+α2(ρ2(0)ρ2(0)+1)j(1ρ2(0)+1))-1),000000000000000000000j=0,1,…,N-1,ciN=ln{+α2(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1))∑r=0N-i(α1(ρ1(1)ρ1(1)+1)r(1ρ1(1)+1))-1(1-∑r=0N-i(α1(ρ1(1)ρ1(1)+1)r(1ρ1(1)+1)+α2(ρ2(1)ρ2(1)+1)r(1ρ2(1)+1)))×(1-∑r=0N-i(α1(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1)+α2(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1))∑r=0N-i(α1(ρ0(1)ρ0(1)+1)r(1ρ0(1)+1))-1},i=1,2,…,N,c0N=c1N.


## 4.2. Mixed Erlang-k Service Model: The M/mixEk/1 Queue

This is illustrated in Figure2, where, with probability α, the top series of k-1 exponential phases is taken and, with probability 1-α, the bottom series of k phases is taken.Figure 2
Mixed Erlang representation.Theorem 2.
The density function of the service time is given by(24)dB(x)dx=αμ(μx)k-2e-μx(k-2)!+(1-α)μ(μx)k-1e-μx(k-1)!,(x≥0),
where k≥2 is an integer.
In this case, one has(25)kn=(ρρ+k)n(kρ+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ+k)),
where kn=P{narrivalsduringaserviceperiod}=∫0∞e-λt((λt)j/j!)dB(t).Proof.
Consider(26)kn=∫0∞e-λt(λt)nn!×(α(kμt)k-2(k-2)!e-ktxkμ+(1-α)(kμt)k-1(k-1)!e-kμtkμ)d(t)=∫0∞e-λt(λt)nn!α(kμt)k-2(k-2)!e-ktxkμd(t)+∫0∞e-λt(λt)nn!(1-α)(kμt)k-1(k-1)!e-ktxkμd(t)=αkμλn(kμ)k-2n!(k-2)!∫0∞e-λttne-kμttk-2d(t)+(1-α)kμλn(kμ)k-1n!(k-1)!∫0∞e-λttne-kμttk-1d(t)=αkμλnkμ(k-2)n!(k-2)!∫0∞e-t(λ+kμ)tn+k-2d(t)+αkμλnkμ(k-1)n!(k-1)!∫0∞e-t(λ+kμ)tn+k-1d(t)=αλn(kμ)k-1n!(k-2)!∫0∞e-t(λ+kμ)tn+k-2d(t)+αλn(kμ)kn!(k-1)!∫0∞e-t(λ+kμ)tn+k-1d(t)=kμλnkμ(k-1)n!(k-1)!∫0∞e-t(λ+kμ)tn+k-1d(t)⟶GamadistributionΓ(n+k)=αλn(kμ)k-1n!(k-2)!(1λ+kμ)n+k-1Γ(n+k-1)+(1-α)λn(kμ)kn!(k-1)!(1λ+kμ)n+kΓ(n+k)=(n+k-2)!αλn(kμ)k-1n!(k-2)!(1λ+kμ)n+k-1+(n+k-1)!(1-α)λn(kμ)kn!(k-1)!(1λ+kμ)n+k=(n+k-2n)(1λ+kμ)n(1λ+kμ)k-1α(kμ)k-1λn+(n+k-1n)(1λ+kμ)n(1λ+kμ)k(1-α)(kμ)kλn=α(n+k-2n)(λλ+kμ)n(kμλ+kμ)k-1+(1-α)(n+k-1n)(λλ+kμ)n(kμλ+kμ)k.
If the numerator and denominator are divided by μ,
(27)kn=α(n+k-2n)(λ/μ(λ+kμ)/μ)n((kμ)/μ(λ+kμ)/μ)k-1+(1-α)(n+k-1n)(λ/μ(λ+kμ)/μ)n((kμ)/μ(λ+kμ)/μ)k=α(n+k-2n)(ρρ+k)n(kρ+k)k-1+(1-α)(n+k-1n)(ρρ+k)n(kρ+k)k
and cobrackets are taken as follows:
(28)kn=(ρρ+k)n(kρ+k)k-1×(α(n+k-2)!n!(k-2)!+(1-α)(n+k-1)(n+k-2)!n!(k-1)(k-2)!(kρ+k))=(ρρ+k)n(kρ+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ+k)).
Here
(29)Pij(ρ)=kj-i+1,i=1,2,…,N,000j=0,1,…,N-1,j≥i-1,P0j(ρ)=kj,j=0,1,…,N-1,P0N(ρ)=1-∑0N-1kn,PiN(ρ)=1-∑0N-ikn,i=1,2,…,N,=1-∑r=0N-i(ρρ+k)r(kρ+k)k-1(n+k-2n)×(α+(1-α)(r+k-1)(k-1)(kρ+k)).
With these values for Pij(ρ), the logarithm of the likelihood ratio will be
(30)lnL=∑i,jnijlnPij(ρ1)Pij(ρ0)lnL=an+∑i,jnijcij,lnL=∑i,jnijln((α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1(n+k-2)!n!(k-2)!)-1((α+(1-α)(n+k-1)(k-1)(kρ1+k))(ρ1ρ1+k)n(kρ1+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ1+k)))×((ρ0ρ0+k)n(kρ0+k)k-1(n+k-2)!n!(k-2)!×(α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1(n+k-2)!n!(k-2)!)-1),lnL=∑i,jnijln((α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1)-1((α+(1-α)(n+k-1)(k-1)(kρ1+k))(ρ1ρ1+k)n(kρ1+k)k-1×(α+(1-α)(n+k-1)(k-1)(kρ1+k)))×((ρ0ρ0+k)n(kρ0+k)k-1×(α+(1-α)(n+k-1)(k-1)(kρ0+k))(kρ0+k)k-1)-1),
where
(31)kn=n(k-1)ln(ρ0+kρ1+k)+∑i,jnijln(ρ1ρ0·ρ0+kρ1+k)n+∑i,jnij×ln(α+(1-α)((n+k-1)/(k-1))(k/(ρ1+k))α+(1-α)((n+k-1)/(k-1))(k/(ρ0+k))).
Here
(32)a=(k-1)ln(ρ0+kρ1+k),cij+bij=ln(ρ1ρ0·ρ0+kρ1+k)j-i+1+ln(α+(1-α)((j-i+k)/(k-1))(k/(ρ1+k))α+(1-α)((j-i+k)/(k-1))(k/(ρ0+k))),i=1,2,…,N,j=0,1,…,N-1,j≥i-1,c0j+b0j=ln(ρ1ρ0·ρ0+kρ1+k)j+ln(α+(1-α)((j+k-1)/(k-1))(k/(ρ1+k))α+(1-α)((j+k-1)/(k-1))(k/(ρ0+k))),j=0,1,…,N-1,ciN=ln{×(α+(1-α)(r+k-1)(k-1)(kρ0+k))(1-∑r=0N-i(ρ0ρ0+k)r(kρ0+k)k-1(r+k-2r))-1(1-∑r=0N-i(ρ1ρ1+k)r(kρ1+k)k-1(r+k-2r)×(α+(1-α)(r+k-1)(k-1)(kρ1+k))(1-∑r=0N-i(ρ1ρ1+k)r(kρ1+k)k-1(r+k-2r))×(1-∑r=0N-i(ρ0ρ0+k)r(kρ0+k)k-1(r+k-2r)×(α+(1-α)(r+k-1)(k-1)(kρ0+k))(1-∑r=0N-i(ρ0ρ0+k)r(kρ0+k)k-1(r+k-2r))-1}i=1,2,…,N,c0N=c1N.


## 5. Numerical Results

Consider aM/H2/1 queue (with poisson arrivals, two-phase-type exponential service, fixed N).Letα1=0.4 and α2=0.6 be the probability for the upper phase and the probability for the lower phase and α=0.005 and β=0.05 the first and second types of errors, respectively. Let N=4 be the capacity of the queuing system.The mean value ofρ is calculated using the following formula:
(33)ρ=λ((α1/μ1)+(α2/μ2)).Suppose that we wish to maintainρ at the level 0.5 and we wish to detect whether its value has increased. Then, the hypothesis test is H0:ρ0=0.5 against alternative H1:ρ1=0.8.Lett0=0,t1,t2,t3,… be a discrete set of the number of customers remaining at points of departure in the M/H2/1 queue. The number of customers remaining at the 48-point departure is given as follows:
(34)000321010332210013211010000000001000001211002100.State space isE={0,1,2,3}. Table 1 shows nij that the number of transitions i→j, and Table 2 shows the value of cij calculating for SPRT.Table 1
i
j
0
1
2
3
0
17
5
1
2
n
0
=
25
1
8
2
1
1
n
1
=
12
2
0
5
1
0
n
2
=
6
3
0
0
3
1
n
3
=
4
n
=
47
ln
L
=
∑
i
,
j
n
i
j
c
i
j
=
-
1.1062.Table 2
The value  ofcij  calculating  for  SPRT.
i
j
0
1
2
3
0
−0.15871
0.114953
0.348346
0.744993
1
−0.15871
0.114953
0.348346
0.744993
2
0
−0.15871
0.114953
0.523806
3
0
0
−0.15871
0.280086
ln
L
=
∑
i
,
j
n
i
j
c
i
j
=
-
1.1062.FromlnB<∑i,jnijcij<lnA decision region, ∑i,jnijcij is between lnA and lnB so that -2.9907<-1.1062<5.2470. Therefore, the system is observed again and the number of customers remaining at the 49 points of departure is given as follows:
(35)00011000000111000000101101021000001000000000210022.Table3 shows nij that the number of transitions i→j.Table 3
i
j
0
1
2
3
0
42
11
4
2
n
0
=
59
1
16
6
1
1
n
1
=
24
2
0
7
2
0
n
2
=
9
3
0
0
3
1
n
3
=
4
n = 96
ln
L
=
∑
i
,
j
n
i
j
c
i
j
=
-
4.35.FromlnB<∑i,jnijcij<lnA, the decision region, ∑i,jnijcij, is small lnB so that -4.3517<-2.9907. Therefore, the hypothesis H0 is accepted.For the computation of the OC and ASN of the SPRT, the largest latent rootλ0(t,ρ) of the P(t) matrix has been computed by fixing the values for t and ρ. Figure 3 shows the graphs of λ0(t,ρ) which have been plotted against the values of t for different values of ρ.Figure 3
Graph ofλ(t,ρ) against t for testing H0:ρ0=0.5, H1:ρ1=0.8.As can be seen from Figure3, there exists one and only one real t0≠0 such that λ0(t0)=1. The derivative of λ(t,ρ) at t=0 has been computed. The OC and ASN functions are then evaluated using the expressions between (13) and (17). The results of the OC and ASN functions for testing H0:ρ0=0.5 against H1:ρ1=0.8 are given in Tables 4 and 5 for α=0.005, β=0.05 and ve α=0.05, β=0.1, respectively.Table 4
The OC and ASN values forα=0.005  and  β=0.05.
ρ
t
(
ρ
)
λ
0
′
(
0
,
ρ
)
L
(
ρ
)
E
(
n
,
ρ
)
0.5
1
−0.0359
0.9951
82.1395
0.6
0.2431
−0.0087
0.8332
184.7906
0.7
−0.4308
0.0165
0.2543
191.3422
0.8
−1
0.0395
0.0496
122.4006Table 5
The OC and ASN values forα=0.05  and β=0.1.
ρ
t
(
ρ
)
λ
0
′
(
0
,
ρ
)
L
(
ρ
)
E
(
n
,
ρ
)
0.50.60.70.8
10.2431−0.4308−1
−0.0359−0.00870.01650.0395
0.95030.70740.30310.0994
55.572985.382680.853760.1863Figures4 and 5 give the OC curve and the ASN curve for two different error probabilities in the same figure.Figure 4
The OC curve for testingH0:ρ0=0.5 against H1:ρ1=0.8.Figure 5
The ASN curve for testingH0:ρ0=0.5 against H1:ρ1=0.8.In this example, we wish to maintainρ at the level 0.5 and we wish to detect whether its value has increased. Then, the hypothesis test is H0:ρ0=0.5 against alternative H1:ρ1=0.8. We accepted H0 hypothesis. There is no need for any change in the system. So traffic intensity (ρ) is the desired level. If any change in ρ, the queueing system would need to be organized. Namely, if ρ has shifted to ρ1 from ρ0 (ρ1>ρ0), the control action could be to increase the mean μ or add an extra server. If ρ has shifted to ρ1 from ρ0 (ρ1<ρ0), the control action could be to decrease the mean μ or reduce the number of servers.We can understand from Table4 that the number of samples required for this hypothesis to be accepted is 82.1395 and the number of samples required for it to be rejected is 122.4006.As can be understood from the tables and Figures4 and 5, the smaller the error probability needed, the more the sampling number and power of a test are decreased.


## 6. Conclusions

The objective of the control technique is to detect changes in the traffic intensityρ as quickly as possible, then take appropriate corrective action, and determine how much of a sample size is needed in the applications. Thus, the sequential probability ratio test provides a saving of up to fifty percent in the sample size according to traditional methods. At the same time, the use of SPRT is easy for observing only the number of customers in the system at successive departure epochs 𝒬n, which are embedded Markov points.In this paper, SPRT is theoretically investigated for two different phase-type queueing systems which consist of Hyperexponential and Mixed Erlang. Also the necessary coefficients have been obtained for SPRT.In the application part, the queue lengths have been generated randomly for different arrival rates(λ), different service rates (μ), and selected system capacities (N) by using MATLAB 7.10.0 (R2010a) programming. Traffic intensities have also been calculated. With the obtained data and predetermined α and β, a simple hypothesis has been established. Accepting or rejecting the hypothesis has been examined by SPRT. Afterward, the largest latent root λ(t,ρ) of the P(t) matrix has been computed by fixing the values for t and ρ. Their graphs have been drawn. It has been found that there exists one and only one real t0≠0 such that λ0(t0)=1. The OC and ASN have been calculated with obtained values, and their graphs have been drawn by Microsoft Office Excel 2007 programming.As a result, the length queue observed at the departure points can be maximumN, where ρ>1, and the length queue observed at the departure points can be minimum N, where ρ<1.According to the givenα and β, differences in sample size have been observed. It has been observed that when the error probability decreased the sample size grows and the power of test (1-β) increases. Here, the value of μ needed to calculate the traffic intensity was taken from the mean μ (33). If you recall, the purpose of a control technique is to detect changes in the traffic intensity ρ. When ρ has shifted to ρ1 from the design level ρ0 (ρ1<ρ0), an appropriate action is taken to bring ρ1 back to the design level ρ0. In the same way, when ρ has shifted to ρ1 from the design level ρ0 (ρ1>ρ0), an appropriate action is taken to bring ρ1 back to the design level ρ0. Consequently, the control action could be to increase mean μ or decrease mean μ for the phase-type queueing systems.



---

*Source: 241241-2013-07-21.xml*