are nonsingular, which implies that inverses of those matrices can be determined.
The computation of the matrix R is by means of the iterative procedure [5].
The sequence {Rk}k is entry-wise nondecreasing and converges monotonically to a nonnegative matrix R. This follows the fact that B−1 is a nonnegative matrix.
The number of iterations needed for convergence increases as the spectral radius of R increases. We terminate the iteration and return with the solution of R when
kRk+1− Rkk∞6 ε, where ε is a given small constants.
3.2 An algorithm for matrix decomposition
Ramaswami’s formula [7]
Consider computing π such that πQ = 0. That is,
h
‧
It also gives
Then we have
h
which is equivalent to
h
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Then, we have B∗1 = B1 · I.
Next, to determine H and V0, we solve the following equations:
V0− V1H = B,
V1 = C, and
−V0H = A.
From the first two equations, it yields h
π∗ πn+1 i
B0 B1 B−1 V0
= 0.
Then, by solving the following equations
π∗(B0− B1V−10 B−1) = 0 (3.9) and
π0· 1 + π1· 1 + · · · + πn(I − R)−1· 1 = 1, we get π∗ = [π0, π1, ..., πn].
LU factorization
Considering the matrix Q1. Here, we assume that π∗Q1 = 0. The equations are of the homogeneous system. We use LU factorization to obtain π∗ in the following steps.
Step 1: Let the first column of Q1 be replaced by the column vector (1, · · · , 1, (I − R)−1· 1)T.
‧
Then, the modified Q1 is rewritten as a new matrix Q3, and we have π∗Q3 =h
Then, we have
Step 3: Applying Gaussian elimination, we transform Ω and Ω into a zero matrix.
‧
Then it gives
Zn =
Theorem 1. Zn is a nonsingular matrix.
Proof : By Step 1, we know
The solution of π∗ is unique by Matrix Geometric Solution, and Q3 is a non-sigular matrix. We transpose the matrix Q3 to QT3. By Step 3, we determine Zn. Q3 is a nonsigular matrix, so is Zn.
Theorem 2. (Roger and Charles [9]) Let Z ∈ Mm×m, a set of m × m matrices.
There exists permutation matrices D, E ∈ Mm×m, a lower triangular matrix L ∈ Mm×m, and an upper triangular matrix U ∈ Mm×m such that
Z = DLUE.
If Z is nonsingular, one may take E = I and Z may be written as Z = DLU.
Proof :
‧
If rank Z=k, Z has a k-by-k nonsingular submatrix, which may, by permutation of rows and columns, be permuted into the upper left corner. Now apply Theorem D in Appendix B to the upper left corner and apply Theorem LU in Appendix A to achieve a factorization. If Z is nonsingular, Theorem D in Appendix B indicates that permutation on the right is unnecessary in order to apply Theorem D in Appendix B, which verifies the second factorization and completes the proof.
Step 4:
and according to above Theorems 1 and 2, we can infer that as follows Remark 1 and Remark 2. Let Zn({i}), i = 1, · · · , 4(n + 1) be formed with the first i rows
‧
It gives the LU factorization of Zn as follows:
The following algorithm is given for Li and Ui:
Algorithm 1: LU factorization Input U0 = BT00∗∗
for i = 1 : n
do Li = CU−1i−1 do Ui = BTii∗− LiFi end
After completing the LU factorization, the vector π can be obtained via block forward and backward substitution:
Algorithm 2: Forward and backward substitution Input y0 = [1, 0, 0, 0]T
for i = 1 : n
do yi = −Liyi−1 end
do πn = U−1n yn
‧
According the above algorithm, we obtain the stationary probability π∗ = [π0, π1, ..., πn].
Remark 2. If Zn is a 4(n + 1) × 4(n + 1) matrix and singular with some 1 ≤ j ≤ 4(n + 1) such that
det(Zn)({j}) = 0,
then by Theorem 2 there exists a permutation matrix D ∈ M4(n+1)×4(n+1) matrix such that
After completing the LU factorization, the vector π can be obtained via block forward and backward substitution:
Algorithm 1: Forward and backward substitution Input y0 = [DT({4})][1, 0, 0, 0]T
for i = 1 : n
do yi = −Liyi−1
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
end
do πn = U−1n yn for i = n − 1 : −1 : 0
do πi = U−1i (yi− Fi+1πi+1) end
In the above Algorithm, DT({4}) is the first four rows and columns composing a 4 × 4 matrix.
According to the above algorithm, we obtain the stationary probability π∗ = [π0, π1, ..., πn].
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Chapter 4
Inter-Departure times
4.1 Departure process
Characterizing the departure process involves developing an infinitesimal generator for the inter-departure times. The elements needed in this development are the departure-point stationary probabilities d = (d0, d1, d2, · · · ). They are related to the continuous time stationary probabilities π = (π0, π1, π2, · · · ) by the following relationships. Here, we denote the total arrival average rate by λ due to superposi-tion of the two arrival streams, and
λ = p[( 1
λ1p + ( 1 λ1 + 1
λ2)(1 − p)]−1+ p[(1
γ1p + ( 1 γ1 + 1
γ2)(1 − p)]−1.
d0 = π1A10/λ, (4.1)
d1 = π2A21/λ, (4.2)
d2 = π3A32/λ, (4.3)
di−1 = πiAi(i−1)/λ, f or i = 3, 4, · · · , n. (4.4)
The departure process has three different partitions:
First, when there is no customer in the system, a departure has to wait until at least one customer has occurred followed by a service.
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Second, when there is at least 1 and at most (n − 1) customers in the system, the inter-departure time can be a function of single processing customer’s remaining service or it could evolve from a customer arrived and its completion of service.
Third, when a departing customer leaves the system with at least n remaining customers, the minimum of the remaining service time of one of these customers or a complete service time of another customer becomes the inter-departure time.
We find that when there remain at least n customers in the system, the inter-departure time characteristics are the same for all these cases. The infinitesimal generator matrix for the departure process Gn1 is given by
Gn1 =
0 1 2 · · · n − 1 n+
0 B00 C 0 · · · 0 0
1 0 B11 C · · · 0 0
2 0 0 B22 · · · 0 0
... ... ... ... . .. ... ... n − 1 0 0 0 . .. B(n−1)(n−1) Cb
n+ 0 0 0 . .. 0 SSin
,
where only two of these sub-matrices are given by SSin = S1⊕ · · · ⊕ S1,
C = C · 1.b
The probabilities of the departure-process system starting in the various states d = (d0, d1, · · · , dn+) are made up of the departure point probabilities, with
dn+ = (
∞
X
a=n+1
πaA/λ)1.
This series can be written in closed form as dn+ = (
∞
X
a=n+1
πnRa−nA/λ)1
= πn[(I − R)−1− I]A/λ)1.
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
4.2 Moments of inter-departure times
The stationary inter-departure time is of the phase type distribution characterized by [d, Gn1]. Thus, from Neuts [8], the moments of inter-departure times random variable X are given below
E[Xk] = k!(−1)kd(Gn1)−k1, k = 1, 2, . . . , V ar[X] = E[X2] − (E[X])2.
4.3 Lag
kcorrelations between successive depar-tures
The stationary probabilities of the states of the arrival process, θ, are obtained from its Markovian arrival process representation. We define
Abn(n−1) = θ ⊗ (S1o⊕ · · · ⊕ S1o).
To get the lagk correlations of the output inter-departure times, it is necessary to develop the generator matrix cGn segmented into two matrices: the internal tran-sitions(without departures) Gn1 and the matrix containing the departure transition Gn2 such that cGn = Gn1+ Gn2.
The matrix Gn2 is characterized as
Gn2 =
0 1 2 · · · n − 1 n+
0 0 0 0 · · · 0 0
1 A10 0 0 · · · 0 0
2 0 A21 0 · · · 0 0
... ... ... ... . .. ... ...
n − 1 0 0 0 . .. 0 0
n+ 0 0 0 . .. (1 − t)Abn(n−1) tSSout ,
‧
國立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
where t is the probability that departure leaves the system with at least n customers remaining. Define t and SSout as
t =
P∞
a=n+1πa1 πn1 +P∞
a=n+1πa1 = dn+ dn−11 + dn+, SSout = S1o⊕ · · · ⊕ S1o.
Given the Gn1 and Gn2 matrices, from Bodrog, Horvath, Telek [2], and Telek, Horvath [11], the lagk correlation is computed by
lagk= (λ)2d(−Gn1)−1((−Gn1)−1Gn2)k(−Gn1)−11 − 1 2(λ)2d(−Gn1)−21 − 1) .
‧
國立 政 治 大 學