Numerical behavior - E F F I C I E N T M AT R I X-A N A LY T I C S O L U T I O N O F M U L T I

The steps of the presented analysis procedure are significantly less computationally demand-ing than the past methods published in the literature considerdemand-ing the same queuedemand-ing system.

In this section we compare our procedure with three prior methods: the method of [3]

(transformed to continuous time), its improved version published in [45], and the procedure of [48]. Note that the latter two procedures are far less general than [3] or the proposed one.

They can handle only preemptive resume service, they do not analyze the sojourn time at all, and [48] is only able to provide the moments of the number of customers.

Since all involved procedures are exact, only the scalability is investigated, that is, the analysis time as the function of the number of phases.

For this purpose let us define theMMAPmatrices as

D₀⁽^K⁾ =

and matrix D_H⁽^K⁾ is defined similarly. The diagonal entries denoted by• are determined uniquely such that the row sums ofD₀⁽^K⁾+DL(K)+DH(K)are zeroes.

The service times are characterized by order-2PHdistributions with parameters σ_H =^h0.16667 0.83333

having service ratesµL = 2.8andµH = 2. The utilization depends onK, it varies between 0.6and0.75.

Figure34depicts the analysis time required to obtain the first 10 moments of the number of low priority customers in the system in the preemptive case as the function ofK². (This is the only performance measure that is supported by all the procedures). It is clearly visible that the presented method is at least an order of magnitude faster than the prior ones, and is able to solve systems with a large number of phases. No numerical problems were encountered even with the largest model. Additionally, as opposed to [45] and [48], the presented procedure can provide sojourn time related performance measures, and is able to handle the case of non-preemptive service as well.

7.4 departure process analysis of priority qeues

This section describes a method to obtain the multi-class lag-1joint moments of the inter-departure times of priority queues. From these moments it is possible to create aMMAP (based on the results of Section 3.2.3) in order to approximate the departure process in a Markovian way.

The approach to derive the multi-class joint moments of the inter-departure times is similar to the one presented in5.3.3 for the single-class case. For simplicity, we discuss only the

2 In our MATLAB implementation theNAREproblems are solved by the ADDA procedure [86] and the Sylvester equations are solved by thelyapfunction of MATLAB, which is based on the Hessenberg-Schur algorithm [33].

7.4 departure process analysis of priority qeues 113

0 100 200 300 400

0 100 200

Size of the MMAP (K)

Executiontime[s]

Presented Method [48]

Method [45]

Method [3]

Figure 34.: Comparison of the execution times of various procedures

two-class case with preemptive resume priority, but the procedure itself can be extended to handle more general systems as well.

The main idea is that the stochastic behavior of two consecutive departure intervals is inde-pendent on the number of customers in the queue when there are at least two customers. The reason is that the system cannot become idle during the two consecutive departure intervals in this case. As a consequence, we have to distinguish just six cases, as follows:

• 0, 0: the last departure left the system empty,

• 1, 0: at the last departure one high and zero low priority customers are left in the system,

• 1, 1+: at the last departure one high and at least one low priority customers are left in the system,

• 2+, 0+: at the last departure at least two high priority customers are left in the system,

• 0, 1: at the last departure zero high and one low priority customers are left in the system,

• 0, 2+: at the last departure zero high and at least two low priority customers are left in the system.

(Note that there were only three cases to distinguish in the single-class MAP/MAP/1 system.) The analysis method presented in Sections7.1and7.2provides only per-class performance measures and does not allow the analysis of the joint behavior of the priority classes. Hence, a different approach is needed to obtain the probabilities of the above listed these six cases.

The approach we are going to use is based on [3], where the analysis of the discrete time DMAP/PH/1 priority queue is presented. In contrast to [3], here we consider the continuous time model and extend the results of that paper in several ways. We provide new closed formulas and more efficient algorithms than the existing ones.

7.4.1 The MMAP[K]/PH[K]/1 preemptive priority queue as a QBD process

It is possible to define a three dimensionalCTMCto model the queue length behavior. One dimension keeps track of the length of the high priority queue, the second one the length of the low priority queue, and the third dimension describes the phase of the arrivalMMAP together with the phases of the low and high priority servicePHdistributions.

114 analysis of the mmap[k]/ph[k]/1 priority qeue

With proper numbering of the states the structure of the generator of this Markov chain is

where the blocks of the generator are infinite matrices corresponding to the same number of high priority customers but different number of low priority customers and different phases of the arrival and service processes. The blocks ofQare defined as

F =diaghE_Hi (280)

With these definitions matricesE_L(andE_H) contain the transition rates accompanied by low (and high) priority arrivals, whileJ₁⁽^L⁾(andJ₁⁽^H⁾) are the ones accompanied by low (and high) priority services, respectively.

Since the generator is a QBD (with infinite number of phases), the solution is matrix-geometric (see Section4.1.3), thus we have

y_k =y₀R^k , k≥0, (285)

wherey_k is the vector of the steady state probability of the states withkhigh priority cus-tomers. This vector can be partitioned according the number of low priority customers, y_k = {y_k,j, j ≥ 0}, wherey_k,j denotes the vector of steady state probabilities for the states withkhigh andjlow priority customers. Furthermore, we denote the marginal steady state probability vectors of the classes as

y⁽_i^H⁾=

∑

∞ ^y^i,j^, ^y⁽ⁱ^L⁾ ⁼

∑

^∞ ^y^j,i^.

7.4 departure process analysis of priority qeues 115 Due to the definition of the blocks of the generator, both matricesRandGexhibit an upper-block-Toeplitz structure, since the number of low priority arrivals and the number of cus-tomers in the system are independent given the phase of the arrival and service processes (shown in [3]). Thus we have:

As in every homogeneousQBD, these matrices are the minimal non-negative solutions to the matrix quadratic equations (see (123) and (128))

0=B+LG+FG². (287)

0= F+RL+R²B, (288)

Applying the definitions of matrices B,LandF, and exploiting the upper-block-Toeplitz structure of matricesRandGwe can derive relationships for matricesR⁽_i^L⁾andG⁽_i^L⁾,i≥0 (see [3]).

The equations for matricesG_i⁽^L⁾are as follows:

fori=0 : 0= J₁⁽^H⁾+ (E₀+J₀⁽^H⁾)G₀⁽^L⁾+E_HG₀⁽^L⁾², (289) fori>_{0 :} ₀=E_LG_i⁽₋^L⁾₁+ (E₀+J₀⁽^H⁾)G⁽_i^L⁾+E_H

∑

i k=0

G⁽_k^L⁾G⁽_i₋^L⁾_k. (290) MatricesG_i⁽^L⁾have important probabilistic interpretations. Entry(a,b)of matrixG_i⁽^L⁾is the conditional probability that starting from levelnwith the background process being in phase a, (1) the first visit to level n−1 occurs in phase b, (2)ilow probability customers arrive during the first passage time.

The following expressions can be obtained for matricesR⁽_i^L⁾:

fori=0 : 0=E_H+R⁽₀^L⁾(E₀+J₀⁽^H⁾) +R⁽₀^L⁾²J₁⁽^H⁾, (291) fori>0 : 0=R⁽_i₋^L⁾₁E_L+R⁽_i^L⁾(E₀+J₀⁽^H⁾) +

∑

i k=0

R⁽_k^L⁾R⁽_i₋^L⁾_kJ₁⁽^H⁾. (292) Interestingly, summing up equations (290) from i = 1to ∞ and adding (289) leads to a matrix quadratic equation forGˆ⁽^L⁾=_∑^∞_i₌₀G⁽_i^L⁾, since

0= J₁⁽^H⁾+ (E₀+EL+J₀⁽^H⁾)G^ˆ⁽^L⁾+EHGˆ⁽^L⁾². (293) Similarly, the sum of matricesR⁽_i^L⁾, denoted byRˆ⁽^L⁾=_∑^∞_i₌₀R⁽_i^L⁾, can be obtained as the minimal non-negative solution of a matrix quadratic equation as well, as

0=E_H+R^ˆ⁽^L⁾(E₀+E_L+J₀⁽^H⁾) +R^ˆ⁽^L⁾²J₁⁽^H⁾. (294) The generating function of matrix seriesG⁽_i^L⁾defined byG⁽^L⁾(z) =_∑^∞_k₌₀z^kG_k⁽^L⁾will be used several times in the sequel. From (290) and (289) we have thatG⁽^L⁾(z)is the solution of the following matrix-quadratic equation:

0= J₁⁽^H⁾+ (E₀+zE_L+J₀⁽^H⁾)G⁽^L⁾(z) +E_HG⁽^L⁾(z)². (295)

116 analysis of the mmap[k]/ph[k]/1 priority qeue

Finally, based on (285) and the block structure of matrix R the steady state probability vector of the number of high and low priority customers in the system can be expressed as

y_i,j =

∑

j k=0

y_i₋_1,kR⁽_j₋^L⁾_k, i≥1,j≥0, (296)

To completely characterize the joint distribution of the number of high and low priority customers, it remains to derive the stationary probability vector at level0,y₀ ={y_0,j, j≥0}_. 7.4.2 Analysis of level zero

Relations for vectory₀can be derived from the boundary equationsyQ=0as

y₀L₀+y₀RB =0, y_0,01=1−λ⁽^H⁾/µ⁽^H⁾−λ⁽^L⁾/µ⁽^L⁾. (297) Due to the structure of matrices L₀, B and R, (297) is equivalent to the solution of an M/G/1 typeCTMC(see [90]). However, by using the relationRB = FG(see e.g. [57], page 144) it is possible to re-formulate (297) and the corresponding M/G/1 typeCTMCto a more appropriate form. The equations fory₀by usingGinstead ofRare then

y₀L₀+y₀FG =0, y_0,01=1−λ⁽^H⁾/µ⁽^H⁾−λ⁽^L⁾/µ⁽^L⁾, (298) and the generator of the related M/G/1 typeCTMCbecomes

Q0=L0+FG= The solution of M/G/1-type Markov chains is based on their invariant matrixGH₀whose entry in position(i,j)is the probability that starting in phaseiat levelnthe first phase visited at leveln−1is statej[57]. (Hence, this matrix plays the same role as the matrixGofQBDs).

In this particular M/G/1-type system matrixGH₀ is the minimal non-negative solution to the matrix polynomial equation

0= J₁⁽^L⁾+ (E₀+J₀⁽^L⁾)GH₀ +ELGH₀2+EH

∑

∞ i=₀

G_i⁽^L⁾GH₀i·GH₀. (300) The stationary probability vectorsy_0,ican be calculated recursively using the Ramaswami formula [70]. Tailoring it to this particular system gives

y_0,i =

where matricesTiare defined by T_i=

7.4 departure process analysis of priority qeues 117

and, vectory_0,0is the solution of the linear system y_0,0

(Observe that the matrix in the parenthesis is the generator of the background process re-stricted to level(0, 0)).

7.4.3 The joint moments of the departure process

As shown in Section5.3.3, the departure process of a MAP/MAP/1 queue is an infinite state MAP. Similarly, the departure process of the MMAP[K]/PH[K]/1 priority queue is an infinite MMAPas well. ThisMMAPgenerates a high priority (low priority) arrival when a high pri-ority (low pripri-ority) departure occurs in the system. The computation of the joint moments of the departure process,η_i,j⁽^L⁾,η_i,j⁽^H⁾, is rather difficult based on this infiniteMMAP representa-tion. Instead of using this representation directly, we construct a finiteMMAP(considering the 6 cases listed above in the beginning of Section7.4), such that the joint distribution of the first two arrivals of thisMMAP– starting from the appropriate initial distribution – is identi-cal with the one of two consecutive stationary departures of the MMAP[K]/PH[K]/1 priority queue.

The blocks of the matrices of this finiteMMAPare constructed from the block matrices of theQBDdescribing the queue length process ((280)-(283)), such that the transitions between the 6 listed cases are taken into consideration. The initial probability distribution is computed according to the stationary distribution of the queue length process just after a departure considering the 6 listed cases.

The matrices of the resultingMMAPrepresentation are as follows:

H₀ =

118 analysis of the mmap[k]/ph[k]/1 priority qeue

where the diagonal blocks ofH₀are:

M₁= E₀,

MatricesH₀,HL andHH can be interpreted as follows. The first set of states correspond to case0, 0. From this state, a high priority arrival moves the Markov chain to state1, 0, and a low priority arrival to state0, 1. No service events can occur in0, 0. From state1, 0(second set of states) there are transitions to1, 1+due to a low priority arrival (EL), to2+, 0due to a high priority arrival (EH), and to0, 0due to a high priority service (J₁⁽^H⁾). The latter one is accompanied by a high priority departure event, hence the corresponding matrix block is located inH_H. The rest of the blocks can be interpreted similarly.

The steady-state distribution of the MMAP[K]/PH[K]/1 queue just after a departure can be calculated from the stationary distribution as

x_i,j = are computed based on (307) and (296) as

x_0,0= ^y^1,0^J

7.4 departure process analysis of priority qeues 119

x_0,1= ^y^1,1^J

(H)

1 +y_0,2J⁽₁^L⁾ λ⁽^L⁾+λ⁽^H⁾

= (y_0,0R⁽₁^L⁾+y_0,1R⁽₀^L⁾)J₁⁽^H⁾+y_0,2J₁⁽^L⁾

λ⁽^L⁾+λ⁽^H⁾ , (312)

x_0,2₊= ^∑

∞j=₃y_0,jJ⁽₁^L⁾+_∑^∞_j₌₂y_1,jJ⁽₁^H⁾ λ⁽^L⁾+λ⁽^H⁾

= (y⁽₀^H⁾−y_0,0−y_0,1−y_0,2)J⁽₁^L⁾+ (_∑^∞_j₌₀y_1,j−y_0,0R⁽₀^L⁾−y_0,0R⁽₁^L⁾−y_0,1R⁽₀^L⁾)J⁽₁^H⁾ λ⁽^L⁾+λ⁽^H⁾

= (y⁽₀^H⁾−y_0,0−y_0,1−y_0,2)J⁽₁^L⁾ λ⁽^L⁾+λ⁽^H⁾

+ (_∑^∞_j₌₀_∑_k^j₌₀y_0,kR⁽_j₋^L⁾_k−y_0,0R₀⁽^L⁾−y_0,0R⁽₁^L⁾−y_0,1R⁽₀^L⁾)J⁽₁^H⁾ λ⁽^L⁾+λ⁽^H⁾

= (y⁽₀^H⁾−y_0,0−y_0,1−y_0,2)J⁽₁^L⁾+(y⁽₀^H⁾Rˆ⁽^L⁾−y_0,0R⁽₀^L⁾−y_0,0R⁽₁^L⁾−y_0,1R⁽₀^L⁾)J⁽₁^H⁾ λ⁽^L⁾+λ⁽^H⁾

, (313) where vectory⁽₀^H⁾denotes the probability vector that there are no high priority customers in the system, thusy⁽₀^H⁾= _∑^∞_i₌₀y_0,i.

Having computed the vectorvand the matricesH₀,H_H andH_L, the joint moments of the departure process are obtained according to (64), hence

η⁽_i,j^c⁾=i!j!v(−H₀)⁻ⁱ⁻¹H_c(−H₀)⁻^j₁, i,j≥0,c={H,L}. (314) 7.4.4 An efficient, truncation-free procedure

Section7.4.2describes how to calculate the steady state joint distribution of the number of customers, based on which Section7.4.3 derives the joint moments of the inter-departure times. These results, however, are not straight forward to implement in an efficient way. The potential numerical pitfalls are:

• To obtain matrixGH₀ the matrix polynomial equation (300) has to be solved, but this matrix equation has infinitely many terms, and relies on infinitely many elements of matrix seriesG_i⁽^L⁾.

• To obtain vectorsy_0,j,j ≥ 0, matricesTi are required (see (301)), that are defined by infinite summations (302).

• To obtain vectorsy_i,j,i > 0,j ≥ 0, a convolution like formula needs to be evaluated that gets slower and slower with increasingi(see (296)).

• To obtain vectorv, vectors likex₂₊_,0₊are necessary to compute, which involves infinite summations of vectorsy_0,j.

The naive numerical implementation of such a procedure evaluates the infinite summa-tions by truncation. In order to avoid the loss of accuracy, the truncation threshold must be high, depending on the parameters (in particular, the load) of the system, which makes the procedure slow.

120 analysis of the mmap[k]/ph[k]/1 priority qeue

Fortunately, it is possible to develop a numerically efficient procedure that addresses all the above listed critical steps without using any truncation. The efficient solution of (300) is addressed in Section7.4.4.1, while the accurate calculation of the vectorsy_0,i and further related quantities is discussed in Section7.4.4.2.

7.4.4.1 Two fundamental matrices and their relations

There are two matrices that play key roles in the efficient analysis of the system. One of them isG_H₀, that is the solution of (300), and the other one is matrixZ, defined by

∑

∞ i=0

G_i⁽^L⁾G_H₀ⁱ. (315)

Theorem 33. If the algebraic multiplicities of the eigenvalues of matrices G_H₀ andZare one, G_H₀Z=ZG_H₀ holds.

To prove this theorem we first need the following Lemma:

Lemma 7. If the algebraic multiplicities of the eigenvalues of matrixG⁽^L⁾(z)_{and matrix}E₀+ zEL+EHG⁽^L⁾(z)are one, the eigenvectors of G⁽^L⁾(z)and E0+zEL+EHG⁽^L⁾(z)are the same.

Proof of the Lemma. The proof uses the same techniques as in [47].

Let ν_i andu_i be the eigenvalue and the corresponding right eigenvector of G⁽^L⁾(z)(for simplicity we assume thatG⁽^L⁾(z)has distinct eigenvalues). AsG⁽^L⁾(z)satisfies the matrix-quadratic equation of (295),ν_isatisfies

deth

J₁⁽^H⁾+ (E0+zEL+J₀⁽^H⁾)ν_i+EHν_i² i

=0, (316)

and the associated right eigenvectoru_iis the solution of h

J₁⁽^H⁾+ (E₀+zE_L+J₀⁽^H⁾)ν_i+E_Hν_i²

i·u_i =0. (317)

Note that bothν_iand vectorsu_iare functions ofz.

By substituting (284) into (316) and by some basic manipulations we get det

ν_iD₀+ν_izD_L+ν_i²D_H

| {z }

⊕ s_Hσ_H+ν_iS_H

| {z }

M₂

⊗I

=0, (318)

from which it follows thatM₁⊕M₂has a zero eigenvalue. Letδ_jandσ_kdenote the eigenval-ues ofM₁andM₂, respectively. Since the eigenvalues ofM₁⊕M₂areδ_j+σ_k, to have a zero eigenvalue there must existj⁰ andk⁰such thatδ_j⁰ =−σ_k⁰. The eigenvector ofM₁belonging toδ_j⁰ is denoted byθ_j⁰, the one of M₂ belonging toσ_k⁰ is denoted by ψ_k⁰. Let us introduce φ_i =θ_j⁰⊗ψ_k⁰⊗₁.

Now we show thatφ_iis an eigenvector ofG⁽^L⁾(z)associated withν_i, thus it satisfies (317):

J₁⁽^H⁾+ (E₀+zEL+J₀⁽^H⁾)ν_i+EHν_i² i·φ_i

= [I⊗M₂⊗I+M₁⊗I⊗I]·(θ_j0⊗ψ_k0⊗₁)

=θ_j0⊗(σ_k⁰ψ_k⁰)⊗₁+ (δ_j⁰θ_j0)⊗ψ_k0⊗₁=_0.

(319)

7.4 departure process analysis of priority qeues 121 since the service process of the low priority class is stopped during the busy period of the high priority class.

The proof will be similar to the one of Lemma7, using the same techniques as in [47] again.

Letλ_k andv_k be the eigenvalue and the corresponding right eigenvector ofG_H₀ (for sim-plicity we assume thatG_H₀ has distinct eigenvalues). AsG_H₀ satisfies the matrix equation of (300),λ_ksatisfies

and the associated right eigenvectorv_kis the solution of

" By substituting (284) into (321) and by some basic manipulation we get

det

from which it follows thatN₁⊕N₂has a zero eigenvalue. Letα_jandβ_hdenote the eigenval-ues ofN₁andN₂, respectively. Since the eigenvalues ofN₁⊕N₂areα_j+β_h, to have a zero Next, we show thatµ_kis an eigenvector ofZ. Observe thatµ_kis an eigenvector ofZif and only if it is an eigenvector ofG⁽^L⁾(z)|_z₌_λ_k since

122 analysis of the mmap[k]/ph[k]/1 priority qeue

As Lemma7states that the eigenvectors of matrixG⁽^L⁾(z)and matrixE₀+zE_L+E_HG⁽^L⁾(z) are the same, it is enough to prove thatµ_kis an eigenvector of(E₀+zEL+EHG⁽^L⁾(z))|_z₌_λ_k:

(E₀+λ_kEL+EHG⁽^L⁾(λ_k))·µ_k

(D₀⊗I⊗I+λ_kDL⊗I⊗I+

∑

∞ i=0

λⁱ_k(DH⊗I⊗I)(G^˜⁽_i^L⁾⊗I)

·(ζ_j⁰⊗ξ_h⁰)

= [N₁/λ_k⊗I]·(ζ_j0⊗ξ_h0) =α_j⁰/λ_k(ζ_j⁰⊗ξ_h0).

(325)

AsG_H₀ andZhave the same eigenvectors, the same matrix diagonalizes them, consequently they commute.

Note that the theorem can be generalized to the case when the eigenvalues are not distinct, but it requires the detailed discussion of the combination of the multiplicities of the eigenval-ues ofM1andM2(N1andN2) that we neglect here.

Based on this commutative property, the next theorem enables the efficient computation ofG_H₀ andZ.

Theorem 34. MatricesG_H₀ andZsatisfy the following coupled matrix-quadratic equations:

0= J₁⁽^H⁾+ (E₀+J₀⁽^H⁾)Z+E_LG_H₀Z+E_HZ², 0= J₁⁽^L⁾+ (E₀+J₀⁽^L⁾)GH₀ +EHZGH₀+ELGH₀2.

(326)

Proof. To obtain equations forZ, we multiply (290) byG_H₀ⁱfrom the right, sum it fromi=1 to∞, and add (289) to it. By using (315) we get

0= J₁⁽^H⁾+ELZGH0+ (E0+J₀⁽^H⁾)Z+EH

∑

∞ i=0

∑

i k=0

G_k⁽^L⁾G_i⁽₋^L⁾_kGH0

i. (327)

The last term becomesE_HZ²by swapping the sums and exploiting thatG_H₀andZcommute, providing the first matrix quadratic equation. The second matrix quadratic equation can be obtained from (300), when the definition ofZis applied.

Interestingly, the two matrix equations show perfect symmetry. While the solution of cou-pled Sylvester equations has an extensive literature (e.g. [27, 60,89] are recent methods), there are no methods available for coupled matrix quadratic equations (according to our best knowledge). A very simple method is given by Algorithm6. (We are sure that more efficient solution methods can be developed as well, but even this simple method performs very well in our numerical examples.)

Note that Algorithm6needs only successive solution of matrix quadratic equations, that is much more efficient than the direct application of the results of Section7.4.2both in time and in space requirement, since the infinite series ofG_i⁽^L⁾matrices and the solution of (300) are not needed.

Although the results above have been derived in an algebraic way, matricesG_H₀ and Z have important probabilistic interpretations. Let us denote by(nH,nL)the set of states in which there aren_H high andn_L low priority customers in the queue. Then, entry (a,b) of matrixG_H₀ is the conditional probability that starting from state ain(0, 1)the first visit to (0, 0)occurs in stateb. Similarly, the entry (a,b) of matrixZis the conditional probability that

7.4 departure process analysis of priority qeues 123

Algorithm 6Solving the coupled matrix quadratic equations of (326) procedureGH₀,Z= SolveCoupled(E₀,EH,EL,J₀^(H),J₁^(H),J₀^(L),J₁^(L))

k←0 GH₀(0)←I repeat

Solve0=J₁^(H)+ (E0+J₀^(H)+ELGH₀(k))Z^(k+1)+EHZ^(k+1)²forZ^(k+1) Solve0=J₁^(L)+ (E0+J₀^(L)+EHZ^(k+1))GH₀(k+1)+ELGH₀(k+1)2

forGH₀(k+1)

k←k+1

untilkGH₀(k)−GH₀(k−1)k<eandkZ^(k)−Z^(k−1)k<e returnGH₀(k), Z^(k)

end procedure

starting from stateain(1, 0)the first visit to(0, 0)occurs in stateb. Using these probabilistic interpretations it is easy to see thatZandG_H₀ commute: let us investigate the busy period generated by a high and a low priority customer, thus the system is in(1, 1)initially. Since the probability that the first passage to(0, 0)occurs in statebis not affected by the order of service, we immediately have thatG_H₀Z=ZG_H₀.

7.4.4.2 The boundary distribution

MatricesT_i,i ≥0, that are necessary to compute boundary probabilitiesy_0,i, are defined by infinite summations (see (302)). To obtain matricesT_iefficiently we introduce closely related matricesAndefined by

A_n=

∑

∞ i=n

G_i⁽^L⁾G_H₀ⁱ⁻ⁿ, n≥0. (328)

Theorem 35. MatricesA_ncan be obtained recursively as the solutions of discrete Sylvester-type matrix equations

0= (−E0−J₀⁽^H⁾−EHG₀⁽^L⁾)⁻¹ ELAn−1+EH n−1 k

∑

G_k⁽^L⁾A_n−k

−An

+ (−E₀−J₀⁽^H⁾−E_HG₀⁽^L⁾)⁻¹E_HA_nZ

(329)

forn>0, and forn=0we have thatA₀= Z.

Proof. The case ofn=0is trivial (see (315) for the definition ofZ).

To derive equations forn > 0, let us multiply equations (290) byG_H₀ⁱ⁻ⁿfrom the right and sum them up fromi=nto∞. We get

0=EL

∑

∞ i=n

G_i⁽₋^L⁾₁GH₀i−n

| {z }

An−1

+(E₀+J₀⁽^H⁾)

∑

∞ i=n

G_i⁽^L⁾GH₀i−n

| {z }

+EH

∑

∞ i=n

∑

i k=0

G_k⁽^L⁾G_i⁽₋^L⁾_kGH₀i−n.

124 analysis of the mmap[k]/ph[k]/1 priority qeue

The last term can be transformed further by swapping the order of the two summations, lead-ing to

where the commuting property of matricesZ andG_H₀ was utilized. Collecting all parts of the equation gives Sylvester equations, establishing the theorem.

Observe that the Sylvester equation for Andepends only on matrices Ai,i < n, that are already available in stepn.

A consequence of Theorem35is that matricesT_ican be expressed explicitly.

Corollary 18. MatricesT_i,i≥0,can be obtained as

Finally, it is possible to solve the sum of matricesAnexplicitly as well, that will be useful later.

Proof. Summing equations (329) from n = 1 to ∞, swapping the sums, and applying the definition ofAˆ provides the corollary after some simple algebraic manipulations.

Corollary 20. The sum of matricesTi, denoted byTˆ =_∑^∞_i₌₀Ti, is given by

Tˆ =E0+J₁⁽^H⁾+ELGH0+EL+EHA.ˆ (334) Proof. The corollary follows from Corollaries18and19.

7.4 departure process analysis of priority qeues 125

7.4.4.3 Calculating the joint moments of the departure process

The steps to obtain the quantities necessary for the departure process analysis are as follows.

1. First solve matricesG_H₀ andZwith Algorithm6.

2. Solve the matrix quadratic equations to obtain matricesGˆ⁽^L⁾andRˆ⁽^L⁾.

3. Calculate matricesGifori =0, 1, 2based on (289) and (290), matricesRi fori= 0, 1 based on (291) and (292), matrices A_i fori = 1, 2based on (329), and matricesT_i for i=1, 2based on Corollary18.

4. Compute matrixAˆ by solving (333) and matrixTˆ from (334).

After this preparation the required stationary probabilities are calculated.

• Vectorsy_0,0,y_0,1andy_0,2 are computed directly from (304) and (301) using theTi ma-trices computed above.

• To derive vectory⁽₀^H⁾=_∑^∞_i₌₀y_0,i, equations (301) are summed up fromi=1to∞, and (304) is added leading to

∑

∞ i=1

∑

i k=0

y_0,kT_i−k

| {z }

from (301)

+y_0,0T₀−y_0,0J₀⁽^L⁾

| {z }

from (304)

∑

∞ k=0

y_0,k

| {z }

y⁽₀^H⁾

∑

∞ i=k

T_i−k

| {z }

Tˆ

−y_0,0J₀⁽^L⁾,

(335)

implyingy₀⁽^H⁾ =y_0,0J₀⁽^L⁾Tˆ⁻¹.

Using the quantities calculated so far, it is possible to obtain all the six components of vector v, namelyx_0,0,x_1,0,x_1,1₊,x₂₊_,0₊,x_0,1andx_0,2₊. The joint moments of the departure process are given by (314).

Part III

Q U E U E I N G N E T W O R K S

8

Q U E U E I N G N E T W O R K A N A LY S I S B A S E D O N T H E J O I N T M O M E N T S

8.1 traffic based decomposition for the analysis of qeueing networks

In document E F F I C I E N T M AT R I X-A N A LY T I C S O L U T I O N O F M U L T I -T Y P E Q U E U E I N G S Y S T E M S W I T H C O R R E L AT E D T R A F F I C gábor horváth Dissertation submitted to the Hungarian Academy of Sciences for the degree of Doctor of (Pldal 124-0)