210
European Journal of Operational Research 51 (1991) 210222 NorthHolland
Theory and Methodology
Analytical perturbations in Markov chains H . B a r u h a n d T. A l t i o k Rutgers University, New Brunswick, N J, USA
Abstract: The problem of recomputing the steadystate probabilities in a Markov chain is considered, after a small deviation is introduced to the original infinitesimal generator. Approximate expressions are developed to calculate the corresponding perturbed vector of probabilities. Computational stability and accuracy issues associated with the approximation are discussed. As an application, the perturbation approach is incorporated into an optimization scheme to identify workloads that maximize the output rate in a production line, where qualitative and quantitative measures are introduced to control the accuracy of the perturbation. Keywords: Work load allocation, perturbation theory, bowl phenomenon
1. Introduction An important issue in the analysis of Markov chains with large state space is the computational difficulty associated with solving the set of balance equations. For many applications, such as telecommunications, computers, and production systems, the order of the probability vector is expected to be very large. The high dimensionality of the system places restrictions on the methods of solution for the steadystate probabilities. Existing methods of solution can be divided broadly into two categories. In the first category, the solution is based on Gaussian elimination or LU decomposition (e.g., Funderlic and Mankin, 1981; Harrod and Plemmons, 1984; Kaufman, 1983). These methods take into consideration the nature of the transition rate matrix. It is, in general, very sparse, and for certain applications in blockdiagonal form, which can be exploited to facilitate the solution for the transition probabilities (e.g., Altiok and Ranjan, 1989). An indepth treatment of sparse matrices can be found in Schendel (1989). The second category consists of iterative methods to obtain the steadystate probabilities, such as GaussSeidel iteration, or the power method (e.g., Altiok and Stidham, 1983; Jennings and Stewart, 1975; Stewart, 1978). The motivation behind such methods is the sparsity of the transition rate matrix. The rate of convergence is dependent on the ratio of the moduli of the two highest eigenvalues of the rate matrix. Under certain circumstances, this ratio is almost one, such that convergence takes too long, or, in some cases, is not achieved at all. It should also be noted that the rate of conx;ergence is inversely proportional to the order of the rate matrix. Improvement of the convergence rate by aggregation and disaggregation is considered by Haviv (1987). In both categories mentioned above the problem addressed is of computing the steadystate probabilities of a given Markov chain. In several circumstances, one needs to change the governing parameters slightly and to recompute the associated steadystate probabilities. Consider, for example, the problem of Received September 1988; revised July 1989 03772217/91/$03.50 © 1991  Elsevier Science Publishers B.V. (NorthHolland)
H. Baruh, T Altiok / Analytical perturbations in Markov chains
211
finding the optimal values of the parameters of a Markov chain, such that an objective function of interest is optimized. A typical application of this is the optimal allocation of resources in production lines, which is the application considered in this paper. At each step of the optimization, values of a subset of the parameters differ slightly from those in the previous step. This implies that during each iteration the steadystate probabilities need to be computed at least once, especially in projected gradient methods. If the computational difficulties that already exist in finding the steadystate probabilities in largescale Markov chains are considered, the complexity of the problem at hand becomes obvious. The question then arises as to whether it is possible to take advantage of the fact that only small changes have been made in the governing parameters when recalculating the steadystate probabilities. In other words, given a certain set of parameters and associated probabilities, we wish to compute the steadystate probabilities associated with a set of slightly perturbed parameters using a minimal amount of computational effort. We will make use of matrix perturbation theory to achieve this objective. The issue of perturbing the governing parameters of a stochastic process and computing the associated Probabilities based on their previous values has received some attention in the literature (e.g., Ho and Cao (1985) using analytic expressions, and Ho, Eyler and Chien (1979) using simulation). Schweitzer (1968) derives exact expressions for the steadystate probabilities after the system parameters are varied. Meyer (1980) uses this result to calculate bounds on the deviation from the original solution, and Haviv and van der Heyden (1984) and Golub and Meyer (1986) refine these bounds. The difficulty associated with calculating bounds on the probabilities is that it requires substantial computational effort. Note that a perturbation analysis for determining the steadystate probabilities can also be used for sensitivity analyses for cases when the governing parameters are not known accurately. Group inverses and eigenvector derivatives are proposed to analyze the deviations in Markov chains in Meyer and Stewart (1988). A perturbation approach has also been used for a special case of the rate matrix in Hunter (1986). In this paper, we use first and secondorder perturbations of Markov chains to compute probabilities associated with varying sets of parameters. Computational issues associated with the perturbation approach are discussed. As an application, we use the perturbation approach within an optimization procedure to allocate workloads in a production line, such that the output rate is maximized.
2. Matrix perturbation Consider an ergodic Markov chain with the following steadystate balance equations
O'rp = O,
pTb = 1,
(1)
where Q is the transition rate matrix of order n and rank n  1, p is the steadystate probability vector in the form P = [ P l P2"" "Pn] T, in which Pi ( i = 1 , 2 . . . . . n) are steadystate probabilities, and b = [1 1 • • • 1] 7. The transition rate matrix Q is generally sparse and diagonally dominant, and the number of nonzero entries in it is, in general, of order O(n). In certain cases, it can be expressed in blockdiagonal form, which helps computationally in obtaining a solution to equation (1), as will be demonstrated later. When the system order is in the hundreds of thousands, the computational problems associated with solving equation (1) are enormous. Not only one encounters numerical instabilities, because the sum of a very large number of nonnegative terms is equal to one, but computer memory and storage problems also exist. To solve equation (1), there basically are two approaches: The first approach treats the system as a set of linear equations and uses an LU decomposition (Stewart, 1973) or Gaussian elimination. To this end, one approach is to convert (1) into a nonsingular system by replacing one of the redundant rows of QT by the vector [1 1 • • • 1] x, and the corresponding element in the vector 0 by 1. Denoting the resulting augmented rate matrix by R, we can write
Rp = a,
(2)
where a = [0 0 • • • 0 1 0 • • • 0] x. The computational effort here is of order O(n3). Some approaches have
11. Baruh, 72 Altiok / Analytical perturbations in Markov chains
212
been proposed to reduce this effort further by taking advantage of the blockdiagonal form of the rate matrix when it exists (e.g., Neuts, 1981; Altiok and Ranjan, 1989; Yeralan and Muth, ~1987; Jafari and Shantikumar, 1986). The second category consists of iterative techniques that take advantage of the sparsity of Q, such as the power method or GaussSeidel. The convergence rate mostly depends on the ratio of the two highest eigenvalues of Q (Stewart, 1973). Because Q is singular and its diagonal entries are negative numbers (equal to the additive inverse of all the offdiagonal elements), all of its eigenvalues are expected to have negative real parts, except the zero eigenvalue. The procedure in the power method is to multiply QT by a factor a, such that the moduli of the eigenvalues of a Q r are spaced between  1 (approximately) and 0. A common way to select a is by a < 1 / f l , where fl = max I Q~I (i = 1, 2 . . . . . n). Then, one adds p to both sides, which shifts the moduli of the eigenvalues to the range of 0 (approximately) to 1. The resulting relation is (aQ'r + I ) p = p ,
(3)
which can now be solved iteratively for p. As is well known, the iteration converges to the eigenvector corresponding to the eigenvalue of a Q r + I with the highest modulus, which in our case is 1, hence converging to the steadystate probability vector. Assuming that the moduli of the eigenvalues of Q are evenly spaced, one can estimate the ratio ~0 of the moduli of the first two eigenvalues as ~0 = 1 / ( 1

1/n) = 1 + 1/n.
(4)
We will quantitatively substantiate this assumption later. To compute the first eigenvector of a Q T + 1 accurately, ~0k must be large, where k is the number of iterations. Note that as n ~ o¢, lim(1 + 1 / n ) " = e. The number of iterations to achieve convergence to the exact probabilities increases with n. Sometimes, if the ratio of the two highest eigenvalues is very close to one, convergence cannot be achieved and roundoff errors dominate. The computational effort is of the order O ( n * ) computations per iteration, where n* is such that the number of nonzero entries of Q is of order O(n * ). Let us now consider the case where the governing parameters are slightly modified, and the steadystate probabilities associated with the new system are sought. Denoting the new rate matrix by Q T + AT and probability vector by p + Ap we can write (QT + A Q T ) ( p + Ap) = O.
(5)
We wish to explore the possibilities of solving equation (5) by making use of the solution of equation (1), and by using a minimum number of computations. To this end, we assume that AQ v is small and can be expressed as AQ T = eQ1, where the elements of Q] are of the same order as Q, and consider a perturbation expansion of Ap as A p = epl + e:P2 + e3p3 + "'" •
(6)
Introducing equation (6) into (5), and separating into like powers of e, we obtain
e0:
QTp = O,
el : QTpl :   Q l P ' e2 : QTp2 =  Q]P]'
(7) g : QTPi =  Ql Pi l,
As can be seen, evaluation of equation (7) to determine p does not provide additional savings, because for each perturbation order a set of n equations has to be solved. However, one can exploit the form of the
H. Baruh, T. Altiok /Analyticalperturbations in Markov chains
213
matrices Q and Q1 to make use of the perturbation assumption. In many cases the number of governing parameters is relatively low, and the transition rate matrix QV and changes in it can be expressed as Q T = cIEI + C2E2 + . . . + cmEm,
(8a,b)
Q 1 = elEI + e2E2 + . . . + e m E m,
where m denotes the number of governing parameters and c r and er (r = 1, 2 . . . . . m) are their values for QT and Q1, respectively. E~ (r = 1, 2 . . . . . m) are matrices whose forms are known in advance and whose entries generally are either zeros, ones, or simple fractions. One can use this property to lower the number of computations to approximate Ap. Substitution of equation (8b) into equations (7) gives for the first and secondorder perturbations m
QTPl =  ~ eiEiP,
QXP2 = 
i=1
(9)
E ejEjPl. j1
Expanding Pl and P2 as m
Pl
=
P2 = E
erhr, r=l
r=l
m
(10)
E eresgr~, s=l
where h r and grs have the properties hV~b = O, g ~ b = 0 (r, s = 1, 2 . . . . . m), we can obtain the values for h~ and gr, by solving QTh r =  E r p ,
QTg~,=E~h,,
(11)
r = 1, 2 . . . . . m,
r, s = l , 2
. . . . . m.
Equations (11) have a distinct advantage that once h~ and grs are calculated, given any change over the original system parameters, the secondorder perturbation of the system probabilities can be computed using equations (10). It is easy to see that equations (10) require computations of order O(n), which is very suitable for our original objectives of varying the governing parameter several times and computing the associated probabilities. It remains to compare the computational effort associated with the perturbation solution with the effort in obtaining the exact solution, and to assess the feasibility of the perturbation solution. To this end, we can write the relations that need to be solved in order to obtain a perturbation solution as OTp=0, QT[gll
O'r[h I h 2 . . . g12
"'"
glm g 2 1
h,,,l=[E, "'"
p E2p
gram] =  [ f i b ,
... Elh2
EmP], "'"
Elnm
E2nl
"'"
Emnm].
(12)
If a solution such as an LU decomposition is sought, three decompositions are needed. If Eqs. (12) are written in terms of the augmented matrix R, and it is possible to store the inverse of R, then only one LU decomposition is required, together with some multiplications of matrices by vectors (order O(n2)). However, it may not be possible to provide the storage capacities for this, because while R is a sparse matrix, its inverse is, in general, fully populated. In any case, if an LU decomposition is the approach chosen to calculate the probabilities, it is more advantageous to use the perturbation equations. For systems in which the rate matrix can be expressed in blockdiagonal form, such as the production line described later on in this paper, this effort can further be reduced. When a solution to equation (12) is sought using the power method, computation of the perturbation quantities require rn 2 + m + 1 solutions. This implies that if the probabilities are to be calculated more than m 2 + m + 1 times the perturbation approach is useful. However, as stated before, the convergence and stability of the power method and GaussSeidel changes from system to system, reducing the feasibility of such methods. Ultimately, selection of the method of solution for the exact probabilities that one should use, given a certain Markov chain, depends on the nature of the problem. We also need to determine the extent to which the perturbation approach is going to yield accurate results. This, of course, is also problem dependent. To analyze the accuracy of the approximation,
H. Baruh, T. Altiok / Analyticalperturbationsin Markov chains
214
qualitative and quantitative measures need to be developed. In the next section we analyze qualitative bounds. In Section 4 we consider an application, where we incorporate the perturbation approach into an optimization procedure to find the workload allocation for maximum output rate in production lines. Also in that section, quantitative measures are introduced to control the accuracy of the perturbation approach.
3. Perturbation bounds In a perturbation analysis it is desirable to keep track of the magnitude of deviation from the exact solution as the system parameters are varied, in our case of Markov chains. We would like to have bounds on the stationary distribution of the Markov chain due to slight changes in the transition rate matrix. These bounds can in turn be used to indicate the extent of deviation from the original probabilities and to evaluate the accuracy of the perturbation solution. We will make use of the results in Meyer (1980) to analyze the relationship between the norm of Ap and changes in the Markov chain parameters. To this end, we consider equation (2) and denote the change in the augmented rate matrix by AR. We can find Ap using the relation (R + A R ) ( p + Ap) = a, which yields
Ap= (I+
R  1 A R )  ~ R 1 ARR,
(13)
which is very similar to the expression for Ap derived in terms of the rate matrix QT (Meyer, 1980):
Ap=[I+(QT)
# AQT] 1(QT) # AQTp,
(14)
where (QT)# is the group inverse of QT. Note that n  1 rows of AR are the same as the corresponding n  1 rows of AQ T. The rate matrix can be expressed as (Meyer, 1980) QT = p c , p  l , where
and P and C are nonsingular matrices of orders n and n  1, respectively. It follows that (QT)n = pDp1, where D=
oT
.
(16)
It is necessary to determine (QT)# as shown above because Q is singular. In essence, the expression for (QT)# is similar to a singular value decomposition (Stewart, 1973). Given (16), it is shown in Meyer (1980) that II a p II/11P tl < ~¢/(1  ~¢),
(17)
where x = II AQ T II II(QT) # II is the condition number of Q (Stewart, 1973). When equation (13) is used or approximated to compute h p , equation (17) still gives perturbation bounds, except that the condition number now becomes ~ = II AR II II R1 II. Meyer (1980) indicates that using equation (14) to find the norm of Ap leads to more reliable results because, in general, the norm l] R1 II is much larger than [I(QT) # IIIn the perturbation approach outlined in the previous section, as the transition matrix is changed, equation (17) can be used to calculate an upperbound for II Ap II and this value can be compared with the computed value of Ap, approximated as d = e p l I 82p2. If II d II/11P II > to/(1  i¢) then it can be concluded that the perturbation solution is no longer correct. However, if II d II/11P II < x/(1  x), this does not guarantee the accuracy of the perturbation solution. It should also be noted that if an iterative approach is used to calculate the steadystate probabilities, the perturbation bounds cannot be calculated using equation (17). More elaborate perturbation bound based on (17) can be found in Haviv and van der
14. Baruh, T Altiok / Analytical perturbations in Markov chains
215
Heyden (1984), where use is made of Y~Pi = 1 in calculating the bounds. In some cases, bounds for each element of the probability vector can be determined. Our primary interest with perturbation bounds is to investigate the accuracy of the perturbation solution using a minimum number of computations. Evaluations of some of the bounds described above may require as much computation as the solution of (1). We will therefore opt to develop additional qualitative and quantitative bounds and accuracy checks, which are more practical and require less computational effort. However, generation of such bounds is dependent on the problem at hand, so we describe them in the next section within the context of an application.
4. An application As an application of the perturbation solution for the steadystate probabilities we consider the problem of optimally allocating workloads to work stations in a production line to maximize the output rate. The issue of workload allocation is a significant one since it has a direct impact on performance measures in production lines. For instance, bottlenecks are direct consequences of a careless workload allocation. The maximum output rate can be achieved by a perfect balance; however, for practical reasons, perfect balance is very difficult to achieve in real life. Consider a production line consisting of m stations is series with finite capacity buffer storages in between. Items follow the same sequence of stations: that is, once the operation of an item is finished at the kth station the item goes directly to the (k + 1)st station, where it enters into service immediately if the server is free; otherwise, it joins the queue if there is a place in the buffer. If there is no place in the buffer of the (k + 1)st station, the item has to wait at the kth station, which therefore becomes blocked. Blocking continues until a space opens up in the buffer of the (k + 1)st station. We assume that there is always a supply of raw materials available to the first station of the line. Because the last station releases the finished products into a storage of infinite capacity, there is no blocking at the last station. Let N~ be the capacity of the jth station. Processing times are assumed to be exponentially distributed with a mean of sj at station j ( j = 1, 2 . . . . . m). It is also assumed that the expected total work content to be completed is fixed, which can be expressed as the sum of the mean processing times. Further details can be found in Altiok and Stidham (1983). The problem is to compute the optimal values of mean processing times to maximize the output rate, which can be formalized as follows: maximize
J = [1pm(O)]/Sm
subject to
sl+s2+
"'" +sin=K,
(18)
S 1 , S 2 , . . . , S m > O,
where pro(O) is the steadystate probability of having the mth station idle, and K is the total work content allocable to the work stations. If is also possible to have a different objective function, such as profit maximization. The quantity p,,(0) can be found through the analysis of the underlying Markov chain. In the above production line, in order to have a Markovian description the state vector has to keep track of the number of jobs in each station and of their status to determine whether the station is blocked or not. Hence, the state vector can be defined as
~'~= [kl; k2, 12; k3, 13; " ' ' ;
km
1, /ml;
Ira] T,
(19)
where k, = 1 if station i is blocked, k i = 0 if station i is not blocked and l, is the number of jobs currently residing at station i. The corresponding steadystate probability will be p ( k l ; k2, /2; ..; lm)" Considering the above, one can compute pro(0) as
pro(O) = Y'~Y', "'" • p ( k , ;
k 2, /2; . . . . l,,_,; 0).
(20)
H. Baruh, T. A ltiok / Analytical perturbations in Markov chains
216
The transition rate matrix Q can be generated considering all the possible values of I2 and all the possible transitions between the states. The way Q is generated, however, directly affects the solution of the balance equations. In the case when the states in the state vector I2 are ordered lexicographically, the resultant Q matrix becomes a block tridiagonal matrix in the form A~
A
A2 B A
C B
C
Q=
(21) A
N,,,  2 N.,1
N,.
B
C
A
B
Ca
C3
C2
where Nm represents the buffer content of the last work station. The above form of Q also leads to an efficient recursive solution of the balance equations. Let the probability vector associated with row i be P(i), such that the steadystate probability vector p can be partitioned as p = [pT(l ) pT(2 ) •  P T (N,,)] V .
(22)
It is then possible to solve equation (1) by recursively solving for P(i) (i = 1, 2..... Nm) corresponding to each row of QT, from which the individual probabilities can be found. The recursive procedure is as follows: P(1)=(AT)'(AT)P(0)=Z(1)P(0),
Z ( 1 ) = ( A T) I ( _ A T ) ,
e(2) = (A~)'[(BT)Z(1)
A~] e(0) = Z(2)P(0),
e(j)
 1)  c ~] e ( o ) = z ( / ) e ( o ) ,
= (#)
'[(B')z(j
Z(j)=(AT)I[(BT)Z(j1)cT],
Z(2) = ( A T )  ' [ (   B T ) Z ( 1 ) AT],
j = 3 , 4 . . . . . Nm  l ,
P(Nm) = (CT)1 [ (  c T z ( N m  1)] P(0)  Z(Nm) P(0),
Z(Nm) = ( c T )  ' [ (   c T z ( N , ,  1)].
(23)
Finally, P(0) can be found from  B T p ( Nm  1 )  c T p ( Nm  2)  c T p ( Nm) = 0
(24)
 B ~ Z ( N , .  l ) e ( O )  C Z ( N , .  2 ) e ( O )  C?Z(Nm)e(O) = O.
(25)
or
A normalization process should be incorporated into this last equation for the probabilities once "all the vectors are evaluated. Note that A and C2 need to be nonsingular for the procedure to be valid, which is the case for our production line problem. Additional computational savings can be realized when calculating the perturbation parameters by Eqs. (12), because the matrices (AT) a, Z(1), Z(2) . . . . . Z(Nm), B T Z ( N m  1), c T Z ( N m  2), and C3xZ(Nm) need to be computed only once and then used for each order of the perturbation solution. In order to solve the optimization problem we employ a gradient method for nonlinear programming problems with equality constraints (see, for example, Phillips et al., 1976). The procedure can be summarized as follows:
H. Baruh, T. A ltiok / Analytical perturbations in Markov chains
217
(1) Begin with an admissible set of work loads, compute the associated probabilities and corresponding value of the objective function. (2) Compute next partial derivatives of the objective function w.r.t, the work loads, and use these in conjunction with the constraints to find the optimal directions, denoted by vi (i = 1, 2 . . . . . m). Because the partial derivatives cannot be computed in closed form, they need to be approximated. Denoting the objective function by J = (s 1, s 2 . . . . . Sin), we can write OJ/Os,= [J(s,, s2
. . . . .
Si qa . . . . . S i n )  J ( s ) ] / 6 ,
i = 1 , 2 . . . . . m,
(26)
where s = [s 1 s 2 ... sm] T and ~ is the step size, determined in general by trial and error. We observe from the above that the objective function needs to be evaluated m + 1 times in each step of the optimization, which requires the steadystate probabilities to be computed m + 1 times as well. This makes it infeasible to use a LU decomposition or the solution outlined in equations (23)(27) to calculate the probabilities throughout the optimization and necessitates use of the perturbation approach. (3) After the optimal direction is determined, the new workloads s ' are computed using the relation s ' = s + v6, where v = [vl u2 Um]T. Then, the new steadystate probabilities and the corresponding value of objective function are computed. Sometimes J ( s ' ) is less than J ( s ) , in which case the interval 6 is made smaller, and the output rate is recomputed. It follows that at this stage the balance equations need to be solved at least once more, bringing the total to at least m + 2 solutions for the steadystate probabilities for each step of the optimization. We now incorporate the perturbation solution to the optimization.. To accomplish this, one has two options. The first choice is to calculate the exact probabilities at each step of the optimization using a global method such as LU or equations (23)(27), or an iterative approach like the power method, and to then use the perturbation equations to compute the partial derivatives ~ J / ~ s , (i = 1, 2 . . . . . m ) . If m is greater than two, computational savings are realized for a global approach. Of course, one can opt to use the power method to compute the partial derivatives of J, since a good initial estimate is available in the form of the original probabilities. However, for largeorder systems it is very likely that iterative methods to solve the balance equations lead to convergence problems even when the initial estimates are close to the actual solution. In the second choice, which is the approach proposed here, after the initial exact probabilities are computed by a global method, the perturbation relations are used to compute both the partial derivatives and the exact probabilities throughout the optimization. This approach eliminates the need to calculate the exact probabilities at each iteration, thus reducing the computational effort even more. However, if the governing parameters sl, s 2. . . . . s,, become very different from their original values, the approximate expressions for the probabilities and the objective function lose their accuracy. This may cause convergence to a nonoptimal set of governing parameters. Under such circumstances one must recompute the exact probabilities using the current set of s 1, s 2. . . . . s,, and begin the perturbation analysis all over. It follows that the total computational effort for this approach depends on how m a n y times the optimization problem needs to be reset by recomputation of the exact probabilities by a global method. The decision to recompute the exact probabilities can be made only if the accuracy of the perturbation is monitored. In the previous section, a qualitative way of analyzing this accuracy was described, where the norm of the probability vector was considered. For the problem of maximizing the output rate in production lines, we develop quantitative guidelines to control the accuracy of the perturbation. We build a number of checks into the optimization and observe them at every iteration. If a flag is triggered, either the step size is reduced or the steadystate probabilities are recomputed using an exact method. Following is the list of checks: (a) Because the steadystate throughput of each station will be the same, the objective function J can be calculated using the throughput of any of the m stations. Equation (18) is one way, on which the optimization is based. However, because of the approximate nature of the solution, the calculated values for J are not the same. Denoting these quantities by/1 (k), J2(k) . . . . . J,,(k ), where k denotes the iteration " ' "
H. Baruh, T. Altiok /Analyticalperturbations in Markoo chains
218 number and
J l ( k ) = [1  p a ( 0 , k )  p l ( b ,
k)]/Sl(k),
J 2 ( k ) = [a  p 2 ( 0 , k )  p z ( b , k ) ] / s 2 ( k ) ,
Jm(k) = [1 pm(O,
(27)
k)]/sm(k),
where pi(0, k) (i  1, 2 . . . . . m) is the probability of having the ith station empty in the kth iteration, and pi(b, k) (i = 1, 2 . . . . . m) is the probability of the ith station being blocked. Note that Jm(k) is the objective function defined by equations (18). A flag is triggered when (i) J ~ ( k + l ) < J , ( k ) , i = 1 , 2 . . . . . m  l , or (ii) J~(k )/Jj( k ), i, j = 1, 2 . . . . . m, is greater than a user defined threshold. (b) The sum of the probabilities may not be equal to 1, due to the inaccuracy of the perturbation solution. (c) The optimization method generates an optimal direction based on the partial derivatives of J. This direction can be considered as a derivative of a kind, so that it should vanish as the optimal solution is approached. The procedure is stopped if this value falls below a threshold value. Similarly, it may be necessary to alter the value of 3 in (26), as explained in step (3) of the optimization. If 3 becomes lower than a threshold, the iteration is stopped. (d) The total number of iterations can be used to stop the procedure. (e) A range for the value of pm(O) can be estimated, and if the perturbation solution exceeds this range a flag is triggered. (f) To preserve the accuracy of the perturbation solution, the workloads should deviate very much from their initial values. The iterations may be stopped when sj(k)/sj(O) ( j = 1, 2 . . . . . m) is not within an allowable range. (g) The optimal allocation of the workloads has been shown to have a special form, known as the 'bowl phenomenon' (Hillier and Boring, 1966), which allocates less work to the stations in the middle. This property can be used to monitor the values of the Workloads. F o r example, for a threestation production line, one should always have s I > s2, and s 3 > s 2, and define thresholds for s2/s 1 and s2/s 3. 5. Results As an illustrative example, we consider a production line with three work stations, each having exponential processing times with means Sl, s2, and s 3, respectively. The two buffers have capacities of N 1 and N 2. The state vector for this system is of order n = ( N 1 + 1 ) * ( N 2 + 2) + N 2 + 1 (Altiok and Stidham, 1983). We select the parameters as N~ = 10, N 2 = 5, such that n = 83 and Esj = 25. First, we consider the problem of computation Of the steadystate probabilities, in an effort to determine whether to use the power method, or an L U decomposition or the approach outlined in equations (23)(27). To investigate convergence of the power method, we calculate the eigenvalues of aQV + 1 (equation (3)). Table 1 compares the modulus of the second eigenvalue ~k2 (the first eigenvalue is X1 = 1) versus values of a, where the initial distribution of the workloads is taken as even. Using the formula a = 1/fl, where Table 1 Modulus of the second highest eigenvalue for varying values of a a 1.50 2.00 2.50 2.78 2.80
~.2 0.99013 0.98685 0.98356 0.98173 0.98158
H. Baruh, T. Altiok / Analytical perturbations in Markov chains
219
Table 2 Secondorder perturbation, expected total work load = 25. The initial expected work loads are 8.3333, 8.3333, 8.3333 si
s2
s3
J
J1
J2
J3
9.2836 9.1000
0.10118 0.10160
0.10204 0.10206
0.10184 0.10197
0.10181 0.10193
8.7461 8.8658
0.10176 0.10180
0.10172 0.10179
0.10175 0.10179
0.10175 0.10179
Ca&u~te exactprobabihties 7.9352 8.0981
7.7812 7.8019
Reca&u~te exactprobabilities 8.2832 8.2044
7.9707 7.9298
Conuergence
/3 = max IQii I yields the value of a = 2.78. The corresponding value of the second eigenvalue is X 2 = 0.98173, which is very close to the value of 1 / w estimated by (4) as ~0 = 1 + 1/n, so that 1/o~ = 1  1/n = 1  1 / 8 3 = 0.9880. Note that using a value of a larger than 2.80 yielded systems whose largest eigenvalue had a modulus higher than 1. Also, it was observed for all values of a that the moduli of the eigenvalues were not evenly distributed between 0 and 1, slowing convergence even more. We thus conclude that, for the particular case considered, an iterative approach is not satisfactory, and consequently compute all initial probabilities and perturbation parameters h r and grs using an LU decomposition or the approach outlined in equations (23)(27). It should also be noted that the eigenvalue analysis described above was conducted for systems of varying orders and for other buffer capabilities. In all cases, the modulus of the second eigenvalue was very close to (sometimes slightly smaller and sometimes larger) 1  1/n. In addition, several very closelyspaced eigenvalues were observed, which slows convergence and creates roundoff errors.
Table 3 Secondorder perturbation, expected total work load = 25. The initial expected work loads are 8.6667, 7.6667, 8.6667 Sl
s2
s3
J
Jl
J2
.]a
8.9205 8.8432 8.8476
0.10174 0.10180 0.10179
0.10175 0.10177 0.10178
0.10173 0.10174 0.10175
0.10175 0.10178 0.10178
8.8684
0.10180
0.10179
0.10179
0.10179
Ca&u~te exactprobabi~ties 8.3251 8.1964 8.2378
7.7544 7.9603 7.9146
Reca~u~te exactprobabi#ties 8.1884
7.9432
Convergence
Table 4 Secondorder perturbation, expected total work load = 25. The initial expected work loads are 8.0000, 9.0000, 8.0000 Sl
S2
S3
J
J1
J2
J3
8.6580 9.3991 9.6051
0.10149 0.10081 0.09994
0.10176 0.10311 0.10330
0.10174 0.10287 0.10295
0.10168 0.10242 0.10247
9.0379 8.8558
0.10168 0.10180
0.10169 0.10180
0.10178 0.10203
0.10178 0.10201
Calcu~te exactprobabirties 7.9631 7.8684 7.7239
8.3789 7.7324 7.6710
Recalcu~te exactprobabi~ties 8.0527 8.1974
Convergence
7.9094 7.9468
1t. Baruh, T. Altiok / Analytical perturbations in Markov chains
220
Table 5 Secondorder perturbation, expected total work load = 15. The initial expected work loads are 5.0000, 5.0000, 5.0000
SI
$2
$3
J
J1
J2
J3
5.3419 5.4122
0.16962 0.16950
0.16997 0.17008
0.16989 0.16996
0.16988 0.16993
5.2809 5.3180
0.16965 0.16966
0.16963 0.16965
0.16964 0.16966
0.16964 0.16966
Ca&u~te exactprobabi#ties 4.8567 4.8422
4.8013 4.7456
Reca&u~te exactprobabi#ties 4.9298 4.9147
4.7893 4.7673
Convergence
Table 6 Secondorder perturbation, expected total work load = 15. The initial expected work loads are 5.3333, 4.3333, 5.3333
sl
s2
s3
J
J1
J2
J3
5.3842 5.2021 5.2159
0.16927 0.16949 0.16952
0.16925 0.16977 0.16979
0.16923 0.16971 0.16972
0.16926 0.16982 0.16985
5.3681 5.2999 5.3181
0.16961 0.16966 0.16966
0.16964 0.16966 0.16967
0.16963 0.16966 0.16966
0.16962 0.16966 0.16966
Ca~u~te exactprobabi~ties 5.0890 4.9513 4.9342
4.5268 4.8466 4.8499
Reca&u~te exactprobabi~ties 4.8764 4.9260 4.9153
4.7554 4.7742 4.7667
Convergence
We next consider optimization of the workloads using the secondorder perturbation. The value of & used in obtaining approximate partial derivatives of the objective function, is chosen as 6 = 0.004. Tables 2 and 3 show the results of the optimization procedure for different sets of initial values for s 1, s 2, and s 3 (even distribution, bowl phenomenon). The column labeled J contains the actual value of the objective function. As can be observed, the results are similar. Table 4 gives the results when the initial distribution
Table 7 Secondorder perturbation, expected total work load = 15. The initial work loads are 4.6667, 5.6667, 4.6667
sl
s2
s~
J
J1
J2
J3
4.8564 5.0578 5.2766 5.5223 5.8067
0.16539 0.16782 0.16918 0.16882 0.16600
0.16543 0.16821 0.17065 0.17274 0.17446
0.16543 0.16822 0.17065 0.17267 0.17416
0.16542 0.16811 0.17027 0.17171 0.17229
5.6073 5.3816 5.2502
0.16833 0.16958 0.16958
0.16833 0.16966 0.16968
0.16836 0.16996 0.17039
0.16836 0.16994 0.17031
5.3183 5.3146
0.16966 0.16966
0.16966 0.16966
0.16966 0.16966
0.16966 0.16966
Ca&u~te exactprobabi#ties 4.6871 4.7017 4.7051 4.6865 4.6229
5.4566 5.2405 5.0182 4.7912 4.5704
Reca~u~te exactprobabi~ties 4.7361 4.8716 4.9927
4.6566 4.7468 4.7570
Reca~u~te exactprobabilities 4.9350 4.9215
Convergence
4.7467 4.7640
221
H. Baruh, T. Altiok / Analytical perturbations in Markov chains
Table 8 Firstorder perturbation, expected total work load =15. The initial expected work loads are 5.0000, 5.0000, 5.0000 sl
s2
s3
J
J1
J2
J3
5.3423 5.6235
0.16962 0.16820
0.17044 0.17156
0.17127 0.17418
0.17116 0.17357
5.3998 5.1755
0.16955 0.16938
0.17029 0.17237
0.17033 0.17256
0.17036 0.17269
5.4858
0.16919
0.17028
0.17068
0.17085
5.2370
0.16958
0,17045
0.17056
0.17065
Calculate exact probabilities
4.8567 4,7697
4.8010 4.6068
Recalculate exact probabilities
4.8911 5.0136
4.7091 4.8108
Recalculate exact probabilities
4.8128
4.7014
Recalculate exact probabilities
4.9615
4.8015
No convergence after 4 tries
Table 9 Firstorder perturbation, expected total work load = 15. The initial expected work loads are 5.3333, 4.3333, 5.3333 Sl
S2
S3
J
JI
.]2
J~
Calculate exact probabilities 5.0886 4.5274 4.8599 4.7282
5.3841 5.4118
0.16927 0.16951
0.16999 0.17220
0.16972 0.17111
0.16963 0.17078
Recah'ulate exact probabilities 5.0186 4.8341
5.1472
0.16927
0.17018
0.17037
0.17047
Recalculate exact probabilities 4.8247 4.7087
5.4666
0.16929
0.17037
0.17085
0.17098
Recalculate exact probabilities 4.9762 4.8100
5.2138
0.16952
0.17041
0.17053
0.17063
No convergence after 4 tries
is taken as an inverted bowl, which represents a poor choice. The results are still satisfactory, b u t convergence to the optimal workloads is slightly slower. We then repeat the o p t i m i z a t i o n process for a total work load of K = 15. Tables 5  7 show the results for different initial values of s 1, s 2, a n d s 3. Similar to the previous case, more realistic initial values give better results. However, the effects of selecting a worse starting p o i n t are more p r o n o u n c e d . This is because the ratios of the o p t i m a l workloads are larger when the expected total work c o n t e n t is lower, i.e., s l / s 2 is greater in the case of K = 15 than the case K = 25. F o r the same reason, we also observe that for all cases of initial admissible workloads, convergence to the o p t i m a l workloads is a little slower in Tables 5  7 t h a n in Tables 2  4 . Finally, we consider o p t i m i z a t i o n of the workloads using a firstorder p e r t u r b a t i o n . The same parameters as the previous case are selected, i n c l u d i n g 8 = 0.004. The results are listed in Tables 8 a n d 9. As can be seen, the results are n o t as accurate as the secondorder case, a n d convergence to the optimal values is slower. These results can be improved a little if a smaller initial value is chosen for 8.
6. Conclusions A p e r t u r b a t i o n technique is developed for c o m p u t i n g the steadystate probabilities in a M a r k o v chain after a small deviation is i n t r o d u c e d to the original infinitesimal generator. Expressions for the first a n d
222
1t. Baruh, T. Altiok / Analytical perturbations in Markov chains
secondorder perturbations are developed and the computational savings resulting from using the perturbation solution instead of the exact solution are analyzed. The perturbation approach is incorporated in an optimization scheme to identify work loads that maximize the output rate in production lines, where qualitative and quantitative measures are developed to control the accuracy of the perturbation.
Acknowledgements This paper has been supported in part by the Center for Computer Aids in Industrial Productivity, New Jersey Department of Higher Education/Rutgers University, by the National Science Foundation grants NSFECS8808772 and NSFMEA8501877, and also by the North Atlantic Treaty Organization grant CRG. 900580.
References Altiok, T., and Ranjan, R., "Analysis of production lines with general service times and finite buffers: A two node decomposition approach", Engineering Costs and Production Economics 17 (1989) 155165. Altiok, T., and Stidham, S., "The allocation of interstage buffer capacities in production lines," l i E Transactions 15 (1983) 292299. Funderlic, R.E., and Mankin, J.B., "Solution of homogeneous systems of linear equations arising from compartment models," S l A M Journal of Scientific and Statistical Computing 2 (1981) 375383. Golub, G., and Meyer, C.D., "Simultaneous computation of stationary probabilities with estimates of their sensitivity", Dept. of Mathematics, North Carolina State University, 1986. Harrod, W.J., and Plemmons, R.J., "Comparison of some direct methods for computing stationary distributions of Markov chains", S l A M Journal of Scientific and Statistical Computing 5 (1984) 453469. Haviv, M., "Aggregation/disaggregation methods for computing the stationary distribution of a Markov chain", S l A M Journal of Numerical Analysis 24 (1987) 952966. Haviv, M., and Heyden, L., v.d., "Perturbation bounds for the ~tationary probabilities of a Markov chain", Advances in Applied Probability 16 (1984) 804818. Hillier, F., and Boling, R., "The effect of some design factors on the efficiency of production lines with variable operation times", J. Industrial Engineering 7 (1966) 651658. Hunter, J.J., "Stationary distributions in perturbed Markov chains", Linear Algebra and Its Applications 82 (1986) 201214. Ho, Y.C., and Cao, X.R., "Performance sensitivity to routing changes in queueing networks and flexible manufacturing systems using perturbation analysis", 1EEE Journal on Robotics and Automation RA1 (1985) 165172. Ho, Y.C., Eyler, A., and Chien, T.T., "A gradient technique for general buffer storage design in a production line", International Journal of Production Research 17 (6) (1979) 537580. Jafari, M., and Shanthikumar, G., "Finite state spatially nonhomogeneous quasibirth and death processes", Department of Industrial Engineering Report, Rutgers University, 1986. Jennings, A., and Stewart, W.J., "Simultaneous iteration for partial eigensolution of real matrices", Journal of Inst. Math. Appl 15 (1975) 351368. Kaufman, L., "Matrix methods for queuing problems", S l A M Journal of Scientific Statistical Computing 4 (1983) 523552. Meyer, C.D., "The conditions of a finite Markov chain and perturbation bounds for the limiting probabilities", S l A M Journal of Algebra and Discrete Methods 1 (1980) 273283. Meyer, C.D., and Stewart, G.W., "Derivatives and perturbations of eigenvectors", S l A M Journal of Numerical Analysis 25 (1988) 679691. Neuts, M.F., MatrixGeometric Solutions in Stochastic Models  An Algorithmic Approach, The Johns Hopkins University Press, Baltimore, MD, 1981. Phillips, D.T., Ravindran, A., and Solberg, J.J., Operations Research, Wiley, New York, 1976. Schendel, U., Sparse Matrices: Numerical Aspects with Applications for Scientists and Engineers, Ellis Horwood, Chichester, 1989. Schweitzer, P., "Perturbation theory and finite Markov chains", Journal of Applied Probability 5 (1968) 401413. Stewart, G.W., Introduction to Matrix Computations, Academic Press, Orlando, FL, 1973. Stewart, W.J., "A comparison of numerical techniques in Markov modelling", Communications of the A C M 21 (1978). Yeralan, S., and Muth, E.J., "A general model for a production line with intermediate buffer and station breakdown", l i e Transactions 19 (1987) 130139.