國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
following. The parameters of arrival processes are given by (λ1, p, λ2) = (10, p, 10), (γ1, p, γ2) = (20, p, 5), and consider ρ = λ/(nµ) from 0.04 to 0.88 with respect to different p’s.
Figure 5.8: Condition number determined by P1
Queue empty rate
In this section, we consider the queue empty rate by comparing MA, RA and LU methods wit forty, fifty, sixty, seventy and eighty servers.
The parameters of arrival processes are given by (λ1, p, λ2) = (10, p, 10), (γ1, p, γ2) = (20, p, 5), and consider ρ = λ/(nµ) from 0.04 to 0.88 with respect to different p’s.
Figure 5.9: The queue empty rate computed by LU method
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Then we use the RA method to compute the queue empty rate. The parameters of arrival processes are given by (λ1, p, λ2) = (10, p, 10), (γ1, p, γ2) = (20, p, 5), and from 0.04 to 0.88 with respect to different p’s.
Figure 5.10: The queue empty rate computed by RA method
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Chapter 6
Conclusion
In this thesis, we present a new computing scheme for computing the stationary probabilities of a phase-type queueing model with multiple servers. The matrix geometric solution procedure has been compared by using Ramaswami’s formula and blocks LU factorization. With LU factorization, an efficient algorithm for solving stationary probabilities is provided to deal with the complex computation of large matrices due to a large number of system states. Through a number of smaller sub-matrices, the state balance equations of a phase-type M AP/M/n queue are solved. Numerical examples are given to demonstrate the proposed matrix geometric solution procedure. Performance measures of these models are also illustrated with a number of approximation and simulation results. As the traffic is light, we find that the stationary probabilities obtained from our approaches and simulations are almost the same. At last, we use two different methods (Matlab and Promodel) to compute the moments of inter-departure times and the variance. We can find the values of two different methods are almost the same.
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Bibliography
[1] Bitran, G.R., Dasu, S., Analysis of the P P h/P h/1 queue. Operations Re-search, Vol. 42, No. 1, pp.158–174, 1994.
[2] Bodrog, L., Horvath, A., Telek, M., Moment characterization of matrix ex-ponential and Markovian arrival processes. Annals of operations Reseach, to appear, 2008.
[3] Chuan, Y.W., Luh, H., Solving a two-node closed queueing network by a new approach, International Journal of Information and Management Sciences, Vol.
16, No. 4, pp. 49–62, 2004.
[4] Curry, G.L., Gautam, N., Characterizing the departure process from a two server Markovian queue: A non-renewal approach, Proceedings of the 2008 Winter Simulation Conference, pp. 2075–2082, 2008.
[5] El-Rayes, A., Kwiatkowska, M., Norman, G., Solving infinite stochastic process algebra model through martix-geometric methods, Proceedings of 7th Process Algebras and Performance Modelling Workshop (PAPM99), J. Hillston and M.
Silva (Eds.), pp. 41–62, University of Zaragoza, 1999.
[6] Gene H. Golub, Charles F. Van Loan, Matrix Computations, 3rd Edition, The Johns Hopkins University Press, 1996.
[7] Latouche, G., Ramaswami, V., Introduction to Matrix Analytic Methods in Stochastic Modeling, ASA-SIAM Series on Statistics and Applied Probability (SIAM), Society for Industrial Mathematics, Philadelphia, PA, 2000.
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
[8] Neuts, M.F., Matrix-Geometric Solutions in Stochastic Models, The John Hop-kins University Press, 1981.
[9] Roger, A.H., Charles, R.J., Matrix analysis, 4th Edition,The Press Syndicate of the University of Cambrige, 1990.
[10] Sikdar, K., Gupta, U.C., The queue length distributions in the finite buffer bulk-service M AP/G/1 queue with multple vacations, Sociedad de Estadistica e Investigacion Operativa, Vol. 13, No.1, pp. 75–103, 2005.
[11] Telek, M., Horvath, G., A minimal representation of Markov arrival processes and a moments matching method. Performance Evaluation, Vol. 64, pp. 1153–
1168, 2007.
[12] Whitt, W., The queueing network analyzer, The Bell system Technical Journal, Vol. 62, No. 9, pp. 2779–2814, 1983.
[13] The MathWorks Company, MATLAB The Language of Technical Computing:
Using MALTAB, Version 6, 2002.
[14] Promodel Corp., Promodel User Guide, Promodel Corp., 2001.
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Appendix A
Theorem LU.(Thm 3.5.2 in [9] )
Suppose that A ∈ Mn×n and that rank A = k. If detA({1, · · · , j}) 6= 0, j = 1, · · · , k, then A may be factored as
A = LU
with L ∈ Mn×n lower triangular and U ∈ Mn×n upper triangular. Furthermore, the factorization may be chosen so that either L or U is nonsingular; both L and U may be chosen nonsingular if and only if k = n, that is, if and only if A is nonsingular.
Proof :
We first show that, under the assumption on leading minors, A({1, · · · , k}) may be factored as L({1, · · · , k}) U({1, · · · , k}), with both nonsingular. It is possible to solve for the relevant entries of L and U, one by one. Let L = [lij] and U = [uij].
Set u11 = 1, and let li1= ai1, i = 1, · · · , k.
Solve for
u1j = a1j
l11, j = 2, · · · , k.
Continue. Set u22 = 1 and let li2= ai2− li1u12, i = 2, · · · , k. Solve for u2j = a2j − l21u1j
l22
, j = 3, · · · , k.
Continue, letting successive diagonal entries of U be 1 and then solving for the next column of L({1, · · · , k}) and then the next row of U({1, · · · , k}).
Each time there is one equation in one unknown to be solved. This equation will be solvable since each liiis nonzero (because det L({1, · · · , i}) × det U({1, · · · , i}) = det A({1, · · · , i})). this completes the factorization of A({1, · · · , k}).
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Partition A. Since rank A = k = rank A11, we see that the rows of [A21 A22] are unique linear combinations of rows of [A11 A12], that is
A21= BA11 and A22 = BA12,
for some uniquely determined B ∈ Mn−k,k. Now partition the desired L and U, noting that nonsingular L11 and U11 have been
U12= L−111A12 and L21= A21U−111.
Then
A22 = L21U12+L22U22= A21U−111L−111A12+L22U22 = BA11A−111A12+L22U22= A22+ L22U22
To complete the factorization, it is necessary and sufficient that L22U22= 0 We may, for example, choose L22(respectively U22) to be any nonsingular lower (respectively upper) triangular matrix in Mn−k we like and choose U22(respectively L22 ) to be 0. Since L11 and U11 are nonsingular, either L or U may be chosen to be nonsingular. If k = n, L = L11 and U = U11 will be nonsingular; if k < n, not both L and U can be nonsingular because A is singular. This completes the proof.
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Appendix B
Theorem D.(Thm 3.5.6 in [9] )
Let Zk ∈ Mk×k be nonsingular. Then there is a permutation matrix D ∈ Mk×k such that
det(DTZk)({1, · · · , j}) 6= 0, j = 1, · · · , k
Note that DTZk is just a recordering of the rows of Zk. Proof :
The demonstrations is by induction on k, If k =1 or 2, the result is clear by inspection; suppose that it is valid up to and including k −1. Consider a nonsingular Zk ∈ Mk×k matrix and delete its last column. The remaining k − 1 cloumns are linearly independent and hence contain k − 1 linearly independent rows. Permute these rows to the first k − 1 positions and apply the induction hypothesis to the nonsingular upper (k − 1)-by-(k − 1) submatrix. This determines a desired overall permutation. Since DTZk is nonsingular, the proof is complete.
‧
Following is the code of the program.
function Rmatrixn(lam1,lam2,gam1,gam2,mu,p,n)
%%%%%%% input data lam1=10;
lam2=10;
gam1=5;
gam2=20;
mu=input(’\n Please input the service rate at all phase, \mu=’);
p=input(’\n Please input the probability that the customer leaves the system after the first phase of the service time, \p=’);
n=input(’\n Please input the number of server, \n=’);
%%%%%%%%%%%%%%%%%% The range of p%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% Define the basic matrices %%%%%%%%%%%%%%
T1=[-lam1 (1-p)*lam1;lam2 -lam2];
T10=[p*lam1;0];
T10p=kron(T10,[1,0]);
T2=[-gam1 (1-p)*gam1;gam2 -gam2];
T20=[p*gam1;0];
T20p=kron(T20,[1,0]);
Tp=kron(T10p,eye(2))+kron(eye(2),T20p);% go into server transposed matrix S1=[-mu];
S10=[mu];%departure server transposed matrix An=kron(n*mu,eye(4));
Bn=[-lam1-gam1-n*mu,(1-p)*gam1,(1-p)*lam1,0;gam2,-lam1-gam2-n*mu,0,(1-p)*lam1;
lam2,0,-lam2-gam1-n*mu,(1-p)*gam1;0,lam2,gam2,-lam2-gam2-n*mu];
C=Tp;
F=[kron(T1,eye(2))+kron(eye(2),T2)+C,ones(4,1)];
%%%%%%%%%%%%%% mean arrival rate%%%%%%%%%%%%%%%%%%%%%%
lam=(inv((1/lam1)*p + ((1/lam1)+(1/lam2))*(1-p))+inv((1/gam1)*p + ((1/gam1)+(1/gam2))*(1-p)))*p %total arrival rate
%%%%%%%%%%%%%% Compute R matrix %%%%%%%%%%%%%%%%%%%%%%
R2=zeros(4,4);
R1=-C*inv(Bn);
i=0;
delta=10;
while delta > (10^-8)
R2=-C*inv(Bn)-R1^2*An*inv(Bn);
%%%%%%% Start Matrix geometric method %%%%%%%%%
t1 = cputime;
Qn=[zeros(4*(n+1),4*(n+1)),ones(4*(n+1),1)];% steady state
Qn(4*(n+1)-3:4*(n+1),4*(n+1)-3:4*(n+1)+1)=[Bn+R2*An,inv(eye(4)-R2)*ones(4,1)];
for k=1:n
Qn(4*(k+1)-3:4*(k+1),4*k-3:4*k)= kron(k*mu,eye(4)) ;
Qn(4*k-3:4*k,4*k-3:4*k)=[-lam1-gam1-(k-1)*mu,(1-p)*gam1,(1-p)*lam1,0;
gam2,-lam1-gam2-(k-1)*mu,0,(1-p)*lam1;
lam2,0,-lam2-gam1-(k-1)*mu,(1-p)*gam1;
0,lam2,gam2,-lam2-gam2-(k-1)*mu] ; Qn(4*k-3:4*k,4*(k+1)-3:4*(k+1))=C ;
‧
sol = [zeros(1,(n+1)*4),1] / Qn % Use the matrix geometric procedure to solve pi_0 pi_1 pi_2... pi_n sum(sol) %The queue empty rate by the matrix geometric
t2 = cputime;
t2-t1% Cpu time of the matrix geometric method
%%%%%%%%%%%%% Start Ramaswami’s formula%%%%%%%%%
t3 = cputime;
A01=[An] ;
B_0(4*(k+1)-3:4*(k+1),4*k-3:4*k)= kron(k*mu,eye(4)) ;
B_0(4*k-3:4*k,4*k-3:4*k)=[-lam1-gam1-(k-1)*mu,(1-p)*gam1,(1-p)*lam1,0;
while delta > (10^-8) R3=-An*inv(Bn)-R1^2*C*inv(Bn);
sol = [zeros(1,4*(n+1)),1] / [P,[ones(4*(n-1),1);inv(eye(4)-R2)*ones(4,1);zeros(4,1)]]% Use Ramaswami’s formula to solve pi_0,...pi_n sum(sol) %The queue empty rate of Ramaswami
t4 = cputime;
t4-t3 % CPU time by Ramaswami
%%%%%%%%%%%%%%%% Starting of verifying Matrix with a replaced column %%%%%%%%%%%
Qn=[zeros(4*(n+1),4*(n+1))];% steady state
Qn(4*(n+1)-3:4*(n+1),4*(n+1)-3:4*(n+1))=[-lam1-gam1-(n)*mu,(1-p)*gam1,(1-p)*lam1,0;gam2,-lam1-gam2-(n)*mu,0,(1-p)*lam1;
lam2,0,-lam2-gam1-(n)*mu,(1-p)*gam1;0,lam2,gam2,-lam2-gam2-(n)*mu]+R2*kron(n*mu,eye(4));
for k=1:n
Qn(4*(k+1)-3:4*(k+1),4*k-3:4*k)= kron(k*mu,eye(4)) ;
Qn(4*k-3:4*k,4*k-3:4*k)=[-lam1-gam1-(k-1)*mu,(1-p)*gam1,(1-p)*lam1,0;
Z=[J1, Qn( 1:4*(n+1),2:4*(n+1) )];
sol = [1,zeros(1,4*(n+1)-1)] / [Z];%
sum(sol) %The queue empty rate
%%%%%%%%%%%%%%% Start LU factorization%%%%%%%%%%%%%%%
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
t5 = cputime;
J=(J1)’;% Transpose
K=[J1, Qn( 1:4*(n+1),2:4*(n+1) )]’;%Use Gaussian elimination temp=[J];
for i=4*(n+1):-1:9 for j=i-1:-1:1
temp(j)=temp(j)-(K(i-4,j)/K(i-4,i))*temp(i);
end temp(i)=0;
end
temp;%Use Gaussian elimination to make first row to be zero
%%%%%%%%%%%% Compute Z_n %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=[temp;K(2:4*(n+1),1:4*(n+1))];%Use Gaussian elimination
%%%%%%%%%%%%% Algorithm 1: LU factorization %%%%%%%%%%
U(1:4,1:4)=S(1:4,1:4);
for i = 2:(n+1)
L((i-1)*4+1:(i-1)*4+4,(i-2)*4+1:(i-2)*4+4) = S((i-1)*4+1:(i-1)*4+4,(i-2)*4+1:(i-2)*4+4)*
/(U((i-2)*4+1:(i-2)*4+4,(i-2)*4+1:(i-2)*4+4));
U((i-1)*4+1:(i-1)*4+4,(i-1)*4+1:(i-1)*4+4) = S((i-1)*4+1:(i-1)*4+4,(i-1)*4+1:(i-1)*4+4)-L((i-1)*
4+1:(i-1)*4+4,(i-2)*4+1:(i-2)*4+4)*S((i-2)*4+1:(i-2)*4+4,(i-1)*4+1:(i-1)*4+4);
end
L;%get L1,L2,...,Ln U;%get U1,U2,...,Un
%%%%%%%%%%%% Algorithm 2: Forward and backward substitution%%%%%%
y=zeros(4*(n+1),1);
for i=2:(n+1) y(1:4,1)=[1;0;0;0];
y(4*i-3:4*i,1)=[zeros(4,1)-L((i-1)*4+1:(i-1)*4+4,(i-2)*4+1:(i-2)*4+4)*y(4*(i-1)-3:4*(i-1),1)];
end
y;%solve y1 y2...yn x=zeros(4*(n+1),1);
x(4*(n+1)-3:4*(n+1),1)=1/(U(4*(n+1)-3:4*(n+1),4*(n+1)-3:4*(n+1)))*y(4*(n+1)-3:4*(n+1),1);
for i=(n+1):-1:2
x(4*(i-1)-3:4*(i-1),1)=1/(U(4*(i-1)-3:4*(i-1),4*(i-1)-3:4*(i-1))) * [y(4*(i-1)-3:4*(i-1),1)-S(4*(i-1)-3:4*(i-1),4*i-3:4*i)*x(4*i-3:4*i,1)];
end
x % Use LU factorization to get pi_0,pi_1,...,pi_n sum(x) %The queue empty rate by the matrix geometric t6 = cputime;
t6-t5 %Cpu time of the matrix geometric method
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Appendix D
List of frequantly used symbols in the thesis symbol its mean
n the number of servers.
S1 an arrival stays in the server.
S1o an arrival finishes the service.
µ the service rate.
λ1 arrival rate in the first phase of the first stream.
λ2 arrival rate in the second phase of the first stream.
γ1 arrival rate in the first phase of the second stream.
γ2 arrival rate in the second phase of the second stream.
p the probability of an arrival to the queue.
B the internal phase changes for the composite arrival process.
C an arrival goes into the system.
A an arrival finishes the service, and departures the system.
θ stationary probability.
λ total arrival rate.
Tm correspond to phase transitions.
Tmo correspond to the rate as arrivals enter the systems phase transitions.
πk the vector of probabilities of k customers in the system.
Q the generator matrix of a continuous time Markovian process.
R means of the iterative procedure.
π∗ h
π0 π1 · · · πn−1 πn i V an upper triangular matrix.
W a lower triangular matrix.
H a submatrix of W.
V0 a submatrix of V.
P B0 − B1V−10 B−1.
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
symbol its mean
Ω a 4 × 4 matrix which the first row is 1, else 0.
Zn a block tridiagonal matrix.
D permutation matrices.
E permutation matrices.
L a lower triangular matrix.
U an upper triangular matrix.
F a submatrix of U.
d0 the departure-point stationary probabilities.
Gn1 the generator matrix without departures for the departure process.
Gn2 the generator matrix with departures for the departure process.
Gcn Gn1+Gn2
t the probability of departure leaves the system with at least n customers remaining.
lagk the lagk correlation.