A DECOMPOSITION METHOD FOR MULTISTATION PRODUCTION SYSTEMS
BY
BARIS TAN
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1994
To
Ebru
and
my parents
ACKNOWLEDGMENTS
I would like to take this opportunity to express my
sincere appreciation to Dr. Sencer Yeralan for his guidance
and his advice in this work.
Special thanks are extended to Dr. Eginhard J. Muth,
Dr. Sherman X. Bai, Dr. Barney Capehart, and Dr. Keith Doty.
iii
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS ................... .................... iii
ABSTRACT ..................................... .. vi
CHAPTER
1 INTRODUCTION............................. ... 1
1.1 Overview, Motivation and Contribution 1
1.2 Background of Production Systems ..... 3
1.3 Organization of this Dissertation.... 5
2 PAST WORK .................................... 8
3 MATRIX POLYNOMIAL FORMULATION OF PRODUCTION LINES
WITH LIMITED BUFFER CAPACITY............... 17
3.1 Introduction .......................... 18
3.2 The Matrix Polynomial Solution Procedure 20
3.3 Properties of the Solution Procedure. 25
3.4 Observations .... .................... 35
3.5 Comparison to Past Work .............. 36
3.6 An Example .............................. 43
4 THE ONE STATIONONE BUFFERINSPECTION STATION
MODEL....................................... 52
4.1 Introduction ......................... 52
4.2 Model Description and Assumptions .... 53
4.3 Model I Reliable Workstation ....... 56
4.4 Model II Unreliable Workstation.... 57
4.5 ClosedForm Expressions for the .....
Performance Measures of the
Workstation ........................ 59
4.6 The Relationship Between the Matrix
Polynomial Solution Procedure and
the Balance Equations .............. 63
4.7 Equivalence of an Unreliable Station
with a Reliable Station ............ 68
4.8 Modeling the One StationOne Buffer .
System as an M/G/1/N Queue ......... 73
4.9 The Relationship Between the Matrix .
Polynomial Solution and the M/G/1/N
Solution ............................ 76
CHAPTER
5 AN APPROXIMATION METHOD FOR MULTISTATION
PRODUCTION LINES WITH LIMITED BUFFER
CAPACITY ................................... 78
5.1 Introduction ........................... 78
5.2 Feedback Control Revisited .......... 82
5.3 Mathematical Programming Formulation. 86
5.4 Properties of the Subproblem ......... 87
5.5 The Convergence Properties of the
Solution Procedure ................. 100
5.6 Numerical Work and Observations ...... 104
5.7 Improvements to the Decomposition
Approximation ....... ................ 107
5.7.1 Dynamically changing the ..
scalar K ................ 114
5.7.2 An approximate formula for
the function K(ai) ....... 115
5.8 Alternate Forms for the Decomposition
Algorithm ......................... 124
5.9 Summary ............................. 126
5.10 Extension I: A Workcell with Rework.. 127
5.11 Numerical Results for the Workcell
with Rework ................ ........ 129
5.12 Extension II: A Production System ...
with Parallel Stations ............. 133
5.13 Numerical Results for the Production
System with Parallel Stations ...... 136
6 SUMMARY AND FUTURE RESEARCH ................. 138
APPENDIX
1 NUMERICAL PROCEDURE TO FIND THE ROOTS OF THE
MATRIX POLYNOMIAL.......................... 141
2 NUMERICAL PROCEDURE TO FIND THE STEADYSTATE
PROBABILITY VECTOR......................... 143
3 EXTENSION TO PHASETYPE ARRIVALS............. 145
4 NUMERICAL RESULTS FOR APPROXIMATE SOLUTIONS OF
PRODUCTION LINES ...... ..................... 154
BIBLIOGRAPHY........................................... 166
BIOGRAPHICAL SKETCH.................................... 173
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
A DECOMPOSITION METHOD FOR MULTISTATION PRODUCTION SYSTEMS
By
BARIS TAN
April, 1994
Chairman: Sencer Yeralan
Major Department: Industrial and Systems Engineering
Markovian models for unreliable multistation production
systems with finite interstation buffers are considered. A
matrix polynomial formulation is developed to find the
steadystate probabilities of a singlebuffer production
line model. Perhaps the most significant result is that the
computational effort in finding the steadystate
probabilities does not depend on the buffer size. The
computational effort only depends on the size of the
submatrices of the transition rate matrix, or more
specifically, on the number of eigenvalues of the associated
matrix polynomial. The solution procedure is thus
convenient for numerical work, especially where a single
buffer model is to be solved repeatedly. Such is the case
in decomposition methods for queueing network models of
manufacturing systems. The value of the methodology is not
in solving a single class of problems but in its flexibility
that applies to a wide range of models. This flexibility
allows various system assumptions to be adopted. For
example, server variability, station breakdown, and
inspection with possible scrapping, may all be implemented.
The properties of the matrix polynomial formulation are
exploited. It is shown that the results obtained in the
past for specific models are the special cases of the
general results presented in this study.
A decomposition method which uses the matrix polynomial
solution procedure is developed to find the performance
measures of heterogeneous multistation production systems.
The method decomposes the system into one stationone buffer
subsystems. The output of each subsystem is fed into the
next one according to some control algorithms. Convergence
properties of the decomposition method are investigated.
The uniqueness and the convergence of the solution procedure
are proven for multistation production lines. The
convergence properties are then used to improve the solution
procedure. Results obtained for the serial arrangement of
workstations are extended to more general arrangements.
vii
CHAPTER 1
INTRODUCTION
1.1 Overview, Motivation and Contribution
Production lines have been the subject of numerous
studies in industrial engineering. A production line is a
serial arrangement of workstations. Items receive an
operation at each successive station. Although topologies of
workstations more general than series arrangements are
increasingly being used, the understanding of production line
behavior is fundamental to the design and operation of such
complex arrangements. Moreover, complex arrangements of
workstations may be decomposed into simple subsystems of
production lines with one or two workstations. These kind of
decomposition approaches seem to yield good approximations.
This dissertation is motivated to provide methodologies
which can eventually be used in practical solutions. Exact
solutions to multistation production systems may forever be
out of the reach of analytical methods. Most analytical
work is limited to two or three station systems. Although
decomposition seems to be a promising approach to evaluating
multistation systems, rigorous studies of its properties are
limited. For example, does a given decomposition algorithm
1
converge to a unique solution or only to a local solution?
How does the convergence rate depend on the system
parameters? Should the decomposition algorithm be modified
for highly unbalanced systems when there are persistent
bottlenecks? Moreover, most decomposition approaches
iteratively solve a specific subproblem consisting of an
equivalent single workstation. It is difficult to model
different operating disciplines such as online inspection
and scrapping or as breakdown and repair mechanisms.
This dissertation, with the aim of ultimately being
able to provide practical solutions to realistic problems,
takes a twopronged approach:
1. Develop an expedient and flexible methodology to
solve the subproblems consisting of a singlebuffer
equivalent station.
2. Rigorously study the properties of the decomposition
approach in order to further streamline the global solution
methodology.
Task 1 leads to a general model and methodology whose
computational effort is independent of the buffer capacity.
Analytical properties and insights are exploited to assist
the solution procedure. The model is based on a matrix
polynomial formulation of the steadystate equations. The
flexibility of this model is essential in studying
production systems with heterogeneous stations.
Task 2 gives rise to a new decomposition approach which
is applicable to a wide range of production systems. The
convergence properties of the model are studied. Again, the
analytical properties and insights obtained are used to
improve the solution procedure.
These two tasks are the major contributions of this
dissertation. Namely, an efficient new method to evaluate
singlebuffer stations, and a new decomposition approach
with predictable convergence behavior are presented.
Although this dissertation takes a rigorous analytical
approach to proposing and evaluating the new methodologies,
its inherent value is in facilitating the development of
useful tools, mostly in terms of software, which are
accurate, efficient, fast, reliable, and robust.
1.2 Background of Production Systems
The two most important performance measures of a
production system are the production rate and the amount of
workinprocess inventory. The production rate is defined
as the average number of items produced per unit time in the
long run. The common measure used to describe workin
process is the expected number of items in each of the
buffers. For operating and design purposes, it is important
to have mathematical models to express the relationship
between the system parameters and the performance measures.
Analytical models prove to be rather complicated due to
station interference among workstations, namely the
phenomena of starvation and blocking. A nonoperating
station with operational neighbors reduces the level of the
downstream buffer and increases the level of the upstream
buffer. For example, a failed station drives the downstream
buffer empty, while driving the upstream buffer full. If
no other breakdowns or repairs take place, a failed station
will either block the upstream station or starve the
downstream station. That is, it forces one or more of its
neighbors to become idle. The adverse effects of blocking
and starving are more pronounced when the buffers have lower
capacities. Blocking and starving are the primary
mechanisms responsible for the reduction of a production
line's efficiency or the production rate.
In general, production lines can be viewed as a series
of queueing systems with finite or infinite waiting rooms
between successive servers. The production rate is a
function of the number of stations, the capacities of
interstation buffers, the distribution of service times,
times to breakdown, and times to repair. A general
production system is viewed as a network of workstations and
buffers.
1.3 Organization of this Dissertation
In Chapter 3, we seek to exploit the special structure
of QuasiBirthDeath Processes (QBDP). More specifically,
we develop a matrix polynomial formulation of the state
homogeneous finitestate QBDPs. The properties of QBDPs are
then used to assist in the solution procedure for computing
the steadystate probabilities. In contrast to past work, we
seek a solution by exploiting the spectral characteristics
of the submatrices of the associated transition matrix. The
spectral properties of the associated matrix polynomial give
insight into the model. Comparisons to past work are
discussed.
To illustrate the flexibility of our approach and
formulation, in Chapter 4 we present a "one stationone
bufferinspection station" model. The model developed in
this chapter will serve two purposes: to illustrate the
general results and properties uncovered in Chapter 3 on a
specific example and to establish a flexible subsystem which
will later be used as a building block in the decomposition
technique to study multistation production systems. This
model is solved by using the matrix polynomial procedure.
Closedform expressions for the performance measures are
derived from the spectral characteristics of the matrix
polynomial. The approximation method for multistation
production system decomposes the system into one stationone
buffer subsystems. Therefore it is important to investigate
this model to get insight and to improve the decomposition
method. This model is also used as a vehicle to illustrate
the relationship between the matrix polynomial solution
procedure and the balance equations. An equivalent reliable
system of an unreliable system is then presented. The
properties of the balanced lines are investigated by using
the equivalent system. It is shown that the one stationone
buffer system can also be solved by using M/G/1/N queues.
The relationship between the matrix polynomial solution
procedure and the solution procedure to solve M/G/1/N queue
is discussed.
In Chapter 5 we present an approximation procedure to
find the performance measures of multistation production
lines and networks. The approximation procedure decomposes
the line into one stationone buffer subsystems which are
discussed in detail in Chapter 4. The input to each
subsystem is adjusted to mimic the output of the upstream
subsystem. The approximation method uses control algorithms
such as proportional control to adjust the input rates into
the subsystems so that the output of each subsystem
converges to the production rate. The results of Chapter 3
and Chapter 4 are used to calculate the performance measures
of individual stations in the network or in the production
line. The most significant theoretical contribution of this
chapter is a formal study of the convergence properties of
the decomposition procedure. Convergence and uniqueness of
the decomposition solution procedure are proven. The
solution procedure is improved by the results obtained from
the convergence properties. It is shown that the
decomposition method is also very effective to analyze
highly unbalanced production systems with persistent
bottlenecks. Results obtained for the serial arrangement of
workstations are extended to more general arrangements. A
model of a workcell with rework and a production system with
routing are given at the end of this chapter.
In Appendix I we present a general numerical procedure
to find the eigenvalues of a matrix polynomial. In Appendix
II a numerical procedure to find the steadystate
probability vector from the matrix polynomial is presented.
The model presented in Chapter 4 is extended to phasetype
arrivals in Appendix III. The numerical results are given
in Appendix IV.
CHAPTER 2
PAST WORK
Production lines play a major role in mass production
systems. Since the pioneering work on production lines by
Koenigsberg (1959) there has been increasing interest in
modeling production lines. Of past work, those by Muth and
by Buzacott deserve further credits. The goal of most of the
studies is the evaluation of the principal performance
measures, namely the production rate and the expected buffer
levels. Since most models focus on steadystate solutions,
computing the steadystate probabilities suffices. Most
approaches lead to Markovian models. Because of
computational complexities, the majority of production line
models are limited to two and threestage models. Although
production lines can be viewed as queues in series with
finite or infinite waiting rooms in front of each server, the
known results in queueing theory are not readily applicable
to production lines.
Most of the past research on production lines assumes
discrete materials flow which leads to Markov chain models.
Most of the studies deal with a particular model that is
specified by a set of assumptions. Almost all models cite
the following assumptions.
1 Station service times, station breakdown times, and
station repair times are statistically independent and all
process parameters are time invariant.
2 Station service times, station breakdown times, and
station repair times are geometrically or exponentially
distributed or have phasetype distributions.
3 Service on an item starts as soon as both the item
and the station are available. Similarly, repair starts
immediately upon breakdown. There are no deliberate pauses
or delays.
In addition, many models further assume that
4 No items are lost.
5 Stations cannot break down when idle or blocked.
Exceptions to these assumptions will be mentioned for
each individual model.
The literature on production lines may be classified
according to whether or not station breakdown is included.
Individual studies in each category can be further classified
according to the number of stations, the presence of
interstation buffers, the model used, the method of analysis,
and the solution procedure. Some of the studies deal with
production line models that have exact analytic solutions.
These studies are limited to two or threestation production
lines. For longer lines, approximate solution procedures are
developed.
For an extensive review of the pertinent literature, the
reader is referred to Dallery and Gershwin (1991). They give
an extensive review of both discrete and continuous lines and
analytical results for general manufacturing systems.
Hunt (1956), Hatcher (1969), AviItzak and Yadin (1965),
and Rao (1975) presented models of two and threestation
production lines with no station breakdown. Hillier and
Boling (1967) examined the effect of line unbalance on the
production rate of production lines with exponential or
Erlang service times and identical buffers. For three
station production lines, they showed that placing a slower
station in the middle improves the production rate. This
ordering is referred to a "bowl phenomenon."
Muth (1979a) introduced a solution procedure to analyze
production lines with general service times and no
interstation buffers. His approach is based on the analysis
of the holding time distribution at the first station. The
holding time is the total time spent by a part on the
station. It is the sum of the processing time and some
blocking time. He used a passage time model in which the
distribution of the interdemand time is obtained as the
solution of an integral equation. He gave specific examples
for two and threestation production lines with uniform and
Erlang distributed service times.
Some earlier work sought empirical solutions to the
performance measures of production lines. Anderson and
Moodie (1969) and Knott (1970) derived empirical formulae for
the utilization of production lines with identical servers.
Lately Blumenfeld (1990) developed an analytical formula for
throughput of a serial production line with variable
processing times and limited buffer capacity. His formula is
based on combination of theoretical results given by Hunt
(1956), Alkaff and Muth (1987) and simulation studies and
regression analyses given by Anderson and Moodie (1969) and
Knott (1970).
Vladzievskii's paper (1952) is perhaps the earliest work
on production lines that are subject to breakdown. After his
work, Sevast'Yanov (1962), Koenigsbergs (1959), Buzacott
(1967), Sheskin (1974), Okamura and Yamashina (1977), Schick
and Gershwin (1978), Gershwin and Berman (1981), Berman
(1982), and Sastry (1985) gave analytical solutions to 2
station production line subject to failure. Sastry and Awate
(1988) presented a model for a production line that includes
rework and inspection. Freeman (1964) used a regression
model based on a set of simulations of a threestation
production line with deterministic service times and
exponentially distributed times to repair and times to
breakdown.
Many production lines can be modeled such that the
associated onestep Markovian transition rate matrix is in
block tridiagonal form.
BO C
A B C
A B C
R = (2.1)
A B C
A BM
Models by Yeralan and Muth (1987), Shantikumar and Tien
(1983), and Hajek (1982) all give rise to transition
probability matrices as in (2.1). A system whose transition
probability matrix is in the block tridiagonal form as in
equation (2.1) is called a QuasiBirthDeath Process (QBDP).
A detailed discussion of QBDPs is given by Neuts (1981).
Neuts studied Markov chains with matrix geometric invariant
vectors. For a class of Markov chains, Neuts showed that the
subvectors of the steadystate probability vector are related
to one another in a matrix geometric fashion. Neuts also
gave a probabilistic interpretation to the rate matrix
associated with the matrix geometric invariant vectors.
Yeralan and Muth (1987) gave eight models which differ
from each other in model assumptions concerning breakdown and
transfer mechanism. They presented three solution procedures
for computing the steadystate probabilities of their models.
These are the recursive solution procedure, the matrix
geometric procedure, and the scalar geometric procedure. A
general state nonhomogeneous discrete flow model is given by
Yeralan and Muth. They present a very general view of the
twostation production line. They described production rate
and average buffer levels as functions of the matrices in
equation (2.1) by establishing a relationship with Neuts'
(1981) matrix geometric systems.
Muth (1979a) showed that the phenomenon of station
breakdown and server variability can be combined.
Muth (1979b) studied general properties of unreliable
product lines without making specific assumptions on service
time distributions. He analyzed several models which differ
in their failure mechanisms.
Muth (1979c) presented the reversibility property. He
stated that the production rate of a line remains unchanged
if the station ordering is reversed.
Ho, Eyler, and Chien (1979) presented a simulation based
study that examines the effect of having one additional
buffer space at a certain location along the production line.
Ammar and Gershwin (1980) presented some fundamental
properties for the conditions of production line (production
network) equivalence.
There are a few models of production lines with
continuous materials flow. These models are sometimes used
to approximate the discrete flow model, but they also have
important applications in the chemical industry.
An early study by Buzacott (1971) presented a mostly
qualitative discussion of a hydraulic model.
Sevast'Yanov (1962), Fox and Zerbe (1974), Murphy
(1974), Wijngaard (1979), Gershwin and Schick (1980), (1983),
De Koster (1987) analyze two and threestation continuous
material production lines with exponentially distributed
failure and repair times.
Yeralan, Franck and Quasem (1986) presented a model of
production lines with continuous materials flow similar to
Wijngaard's model (1979). They provided a general model in
which no restrictive assumptions are made about the
distributions of the station breakdown and repair times when
stations are blocked or starved. Their model is decoupled
into three, one describing the blocked states, one describing
the starved states, and one describing the remaining
(interior) states. They give the production rate and the
expected level of the buffer in closed form.
Exact analytical methods to investigate the performance
of production lines with more than three stations and with
finite capacity interstation buffers are unavailable. There
has recently been a flurry of activity in developing
decomposition algorithms to study multistation production
lines. See Perros (1989) and Dallery and Gershwin (1991) for
a complete list of references. Sevast'Yanov (1962), Gershwin
(1984), Dallery, David, and Xie (1988), Dallery, David, and
Xie (1989) present decomposition methods that approximate the
production line as a series of twostation lines separated by
a buffer. Recently, Dallery and Frein (1993) presented a
unified view of the decomposition methods based on
approximating the production line as a series of twostation
lines. They proved uniqueness and convergence of the
decomposition method for a production line with reliable
exponential servers. Shantikumar and Tien (1983) and Altiok
and Stidham (1983) present recursive solution procedures
based on the balance equations.
Buzacott and Kostelski (1987) used both the matrix
geometric procedure and the recursive algorithm procedure to
compute the steadystate probabilities. They concluded that
the recursive algorithm procedure is better than the matrix
geometric procedure. Liu and Buzacott (1989) developed a
decomposition method based on the concept of zerobuffer
equivalence.
Although the decomposition methods often yield results
comparable to those obtained by simulation, their
applicability is limited to production lines whose stations
are similar to those described by the approximating two
station model. Small variations in the assumptions governing
the operating policies of the stations and of the buffer may
have significant impacts on the performance of a twostation
production line with a finite buffer (Yeralan and Muth,
1987). Thus, in order to construct a decomposition method
for production lines with heterogeneous stations, a simple
but efficient method capable of evaluating the performance
and operating characteristics of a wide range of twostation
models is required. The method developed in Chapter 3 is
proposed as a general approach to evaluate a wide range of
singlebuffer production line models.
CHAPTER 3
MATRIX POLYNOMIAL FORMULATION OF PRODUCTION LINES WITH
LIMITED BUFFER CAPACITY
In this chapter, we consider a generic subsystem of a
production line with a finite interstate buffer. A Markovian
state space model is developed. The steadystate
probabilities of the model are sought. It is known (Neuts
1981) that the transition rate matrix has a nested block
tridiagonal structure, as given in equation (2.1).
Following the spectral characteristics of the associated
transition rate matrix the steadystate equations are
formulated as a matrix polynomial equation. The spectral
properties of the associated matrix polynomial give insights
into the model and suggest expedient solution procedures for
the steadystate probabilities. The model describes an
entire production system when there is only one buffer. It
is shown that for singlebuffer production lines, the
computational effort of the solution procedure is
independent of the buffer size. The formulation and
solution procedure developed is shown to generalize various
past approaches, including the sumofproducts formulation
(Schick and Gershwin 1978, Gershwin and Berman 1981, Sastry
and Awate 1988), balance equations formulation (Okamura and
Yamashina 1977), and M/G/l/N queueing models. More
specifically, it is shown that the critical system
parameters, or latent values, correspond to the eigenvalues
of the matrix polynomial.
3.1 Introduction
Many production lines can be modeled such that the
associated onestep Markovian rate matrix is in block
tridiagonal form, as given in equation (2.1).
For example, models by Yeralan and Muth (1987),
Shantikumar and Tien (1983), and Hajek (1982) all give rise
to transition probability matrices as in equation (2.1). A
system whose transition probability matrix is in the block
tridiagonal form as in equation (2.1) is called a Quasi
BirthDeath Process (QBDP). A detailed discussion of QBDPs
is given by Neuts (1981).
In the development below we seek to exploit the special
structure of QBDPs. More specifically, we develop a matrix
polynomial formulation of the statehomogeneous finitestate
QBDPs. The properties of QBDPs are then used to enhance a
solution procedure for computing the steadystate
probabilities. In contrast to past work, we seek a solution
by exploiting the spectral characteristics of the
submatrices of the associated transition matrix. The
spectral properties of the associated matrix polynomial give
insight into the model. The practical advantage to our
approach is in the reduction of the computational effort in
solving the steadystate probabilities. Most notably, once
the eigenvalues and eigenvectors of the associated matrix
polynomial are obtained, the steadystate probabilities for
any buffer capacity are computed by solving a fixed number
of linear equations, i.e., by solving for the weights.
Decomposition methods perform well, often yielding
results comparable to those obtained by simulation. Past
work commonly starts with a given subsystem with fixed
assumptions. The system is then described as a series of
such subsystems. The rigid assumptions on the subsystem
limit the applicability of decomposition techniques. For
example, if the subsystem conserves parts, so must the
entire system. Modeling a production line with possible
scrapping then requires a new subsystem. It is known that
small variations in the assumptions governing the operating
policies of the stations and of the buffer may have
significant impacts on the performance of a twostation
production line with a finite buffer (Yeralan and Muth
1987). In order to construct a decomposition method for
production lines with heterogeneous stations, a simple but
efficient method capable of evaluating the performance and
operating characteristics of a wide range of singlebuffer
subsystem is required. The method developed in this chapter
is proposed as a general approach to evaluate a wide range
of singlebuffer subsystems or models.
3.2 The Matrix Polynomial Solution Procedure
We consider a production line model such that the
transition rate matrix is in block tridiagonal form as given
by equation (2.1). It is known that the steadystate
probability vector, n, satisfies the following equations.
SR = (3.1)
x u= 1 (3.2)
where R is the rate matrix, u = (1,1,1,...,1)T, and 0 =
(0,0,0,...0)T. The steadystate probability vector is unique
if R is irreducible. Any vector satisfying (3.1) but not
(3.2) is called a nonnormalized steadystate probability
vector. Let g = (g01,gl,g2,..,gM) be such a vector
partitioned according to the partitioning of R. Then
gR = 0 (3.3)
Equation (3.3) is now written in partition form.
g0B0 + glA = 0 (3.4)
giiC + giB +gi+1A = 0 i = 1,2,...,M1 (3.5)
gM1C + gMBM = 0 (3.6)
Equations (3.4) and (3.6) are called the boundary
equations and equation (3.5) is called the interior
equation. In computing steadystate probabilities, our
approach exploits the fact that the interior equations are
repeated Ml times.
A matrix difference equation in the form (3.5) is known
to have a solution in the following form (Lancaster,
Gohlberg and Rodman, 1982):
gi = Xie for i = 0,1,...,M (3.7)
where X is a scalar and e is a vector of the same dimension
as gi. Substituting equation (3.7) into equation (3.5)
yields,
kileC + kieB + Xi+leA = 0 (3.8)
Provided that X # 0 we have
%2eA + XeB + eC = 0 (3.9)
Let us now define the matrix polynomial L(X).
L(X) = X2A + XB + C (3.10)
L(k) is a matrix polynomial whose eigenvalueeigenvector
pairs solve equation (3.5). From equations (3.9) and (3.10)
we get
eL(X) = 0 (3.11)
Any eigenvalueeigenvector pair (Xj,ej) of L(k)
provides a general solution to equation (3.5). The
eigenvalues X are the roots of the characteristic polynomial
det(L(1)) = 0. It is known that det(L(X))=0 has 2n
solutions where n is the size of the submatrices A, B, and
C. The eigenvalues, 11 to k2n are ordered as follows.
l1 ........ ,N finite and nonzero
XN+1 **... N+Z = 0
kN+Z+l,.,k2n = o or +x
where N is the number of finite and nonzero eigenvalues and
Z is the number of zero eigenvalues. Note that Z is the
dimension of the null space of C and the number of
infinitely large (o) eigenvalues 2nNZ, is equal to the
dimension of the null space of A.
A method to find the coefficients of the characteristic
polynomial is given in the Appendix I. This method is
especially convenient for numerical work. The roots of the
characteristic polynomial are obtained by Muller's method,
since all roots, real or complex roots are of interest. The
possibility of having complex eigenvalues is discussed in
Section 3.4.
Any linear combination of the general solutions also
solves the matrix difference equation (3.5). All possible
solutions for equation (3.5) can thus be expressed as
follows:
N
gi = Wjjej i=0,l,...,M (3.12)
3=1
where wj is the weight associated with the jth eigenvalue
eigenvector pair.
The case X=0 only affects go. In addition to the
components in equation (3.12), any eigenvector of C
corresponding to the zero eigenvalue may be included. That
is
N N+Z
g0 Owjj + wjj (3.13)
>1 JN+1
where ejC = 0 for j = N+1,...,N+Z. Note that equation
(3.12) includes equation (3.13), provided that the sum is
computed from j=l to j=N+Z. Similarly, for X=oo, i.e.,
1'=0 for the matrix polynomial L'(k') = XdL(Xl), where d is
the degree of L(1), additional components may have to be
included in gM. Parallel to equation (3.13)
N 2n
gM= Zw j+ Z wjej (3.14)
]=1 j>N,21
where ejA = 0 for j=N+Z+l,...,2n.
More formally, we define a finite Jordan pair (XF,JF)
of L(X) as
XF = [X(kl), X(X2), X(13),...X(Xp)]
JF = diag[J(X1), J(X2), J(X3),..J(Xp)]
where (X(,j),J()j)) is a Jordan pair for every finite
eigenvalue of Xj of L(1) and p is the number of different
eigenvalues of L(X). For a thorough review of matrix
polynomials, the reader is referred to Lancaster, Gohlberg,
and Rodman (1982). In order to solve equation (3.6) we have
to consider an additional Jordan pair (Xo,Jo) of L(X) for
II=oo. (Xoo,J) is obtained from the matrix polynomial L(1)
corresponding to II=o We call (Xo,Jx) an infinite Jordan
pair of L(X). As mentioned (XO,JO) is also a Jordan pair
of the matrix polynomial L'(k') = XdL(k1) corresponding to
X'=0. Equation (3.12) can be expressed in matrix form as
gi = WJFiXF i=0,1,2,...,M1 (3.15)
gM = WJFMXF + WoXJo (3.16)
In the above equation W is a row vector of size N+Z. The
elements of W are the weights corresponding to each
eigenvalueeigenvector pair. Similarly, Wo is a row vector
of size 2nZN. The elements of W. are the weights
corresponding to eigenvectors with eigenvalues IX=oo.
Equation (3.15) with a unspecified vector W is referred
to as the general solution to the matrix difference equation
(3.5). It is known that by choosing various vectors W,
W(JF)iXF spans the entire solution space of the matrix
difference equation (3.5). That is to say, any sequence of
vectors gi, i=0,1,...,M which satisfies (3.5) can be written
in the form of equation (3.15) by selecting the appropriate
W. In order to specify the particular solution, that is a
solution which also satisfy equations (3.4) and (3.6), the
proper weight vectors, W and W. must be calculated by
using the boundary equations (3.4) and (3.6).
3.3 Properties of the Solution Procedure
We now present some properties of the solution
procedure. We make two easily justifiable assumptions.
Assumption 1: The process is ergodic, that is, R is
irreducible.
Assumption 2: The matrix (A + B + C) is an irreducible
transition rate matrix.
If Assumptions 1 and 2 did not hold, the system could
be decoupled into two or more subprocesses. Each such
process would then be examined separately.
Note that Assumption 1 does not imply Assumption 2.
The system, described by R, may be ergodic although A+B+C is
reducible. Consider for example
[a 0 (a + c) 0 c 0
A = B = C =
S[(c + bo) bO (a + bM) bM
bo (c + bo) bM (c + bM)J
where A+B+C=O and thus reducible but R is ergodic.
It should be mentioned that the quantities such as eAu,
eBu, eCu, or some other derivatives which appear frequently
in the following observations, are related to the rate flow
also known as the balance equations. For example, niAu is
the effective rate at which the buffer level changes from i
to il. It will be shown that there exists relationships
between the spectral characteristics of the matrix
polynomial and the balance equations.
Corollary 1:
Au # 0 and Cu # 0 u = [1,1,...,1]T.
Corollary 1 holds since otherwise the process would not
be ergodic and R would be reducible.
Property 1:
Au + Bu+ Cu = 0 where 0 = (0,0,0,...,0)T and is the
same size as the left hand side.
Since R is a transition rate matrix, the elements of
each row of R sum to 0, and thus, the rows of A, B, and C
must also sum to 0.
Property 2:
One is always an eigenvalue of L(k).
Property 2 holds because the rows of L(1) = A + B + C
sum to zero as given in property 2.
Property 3:
gkCu = gk+1Au for k=0,1,2,...,M1 (3.17)
This property is commonly known as the balance of flow
equations. More precisely, equation (3.17) states that the
effective rate at which the buffer level changes from k to
k+1 is equal to the rate at which the buffer level changes
from k+l to k. The proof is by mathematical induction. From
equation (3.4) we have g0BOu = g1Au and since from equation
(2.1), (BO + C)u = 0, we obtain
g0C g = glAu (3.18)
Thus the property holds for k=0. Combining equation
(3.5) for i=l and property 1 yields
g0Cu + glBu + g2Au = 0.
Substituting glAu for g0Cu from equation (3.18) gives
gI(A+B)u + g2Au = 0.
From Property 1, we have (A+B)u = Cu, and thus
g1Cu = g2Au (3.19)
Thus the property holds for k=l. Now assume that the
property holds for kl, i.e.,
gk1Cu = gkAu (3.20)
Again, using property 1 and equation (3.5) for i=k, we
obtain
gkiC1 gk(Au + Cu) + gk+iAu = 0.
Now using equation (3.20) gives
gkCg = gk+lAu (3.21)
Thus the property holds for k. This completes the proof by
induction.
Property 4:
XjejAu ejCu = 0 provided that Xj#l (3.22)
We refer to equation (3.22) as the balance equations in
component form (BECF). Property 4 states that the BECF hold
for all eigenvalueeigenvector pairs with eigenvalue
different from one. Notice that if the eigenvalue has value
one, the corresponding BECF may or may not hold. A more
thorough study for X=l is presented later.
Suppose (Xj,ej), j = 1,2,...,2n are the eigenvalue
eigenvector pairs. Then from equation (3.9) and property 1,
we obtain
kj2ejAu .jej(Au + Cu) + ejCu = 0, or,
Xj(Xj l)ejAu (Xj l)ejCu = 0
(Xj 1) (jejAu ejCu) = 0 (3.23)
Thus provided that hj # 1, XjejAu ejCu = 0, which
completes the proof.
Property 5:
All eigenvalues Xj=l with multiplicity m form a single
Jordan block. That is, there is one genuine eigenvector and
m1 generalized eigenvectors corresponding to eigenvalue
one.
If there exist more than one Jordan block corresponding
to eigenvalues Xj = 1, then A+B+C would be reducible. This
contradicts with the Assumption 2.
Theorem 3.1
The eigenvalue, X1 = 1 has multiplicity m 2 2 if and
only if elAu = elCu where el is the genuine eigenvector
corresponding to eigenvalue one.
Proof:
In order to prove Theorem 3.1, we first enumerate the
eigenvalues. Let l1=1 have multiplicity m 2 2, i.e.,
1 = X2 = ...= Xm = 1 with the corresponding eigenvectors
el, e2, ...., em. The vector el is the genuine eigenvector,
and the vectors e2,...,em are the generalized eigenvectors.
The defining equation for e2 is (see Lancaster, Gohlberg,
and Berman, 1982)
lim[endL(X) + e2L()] =0 (3.24)
1L dX
or,
e1(2A + B) + e2(A + B + C) = 0. (3.25)
Postmultiplying equation (3.25) by u and using property
1 gives elBu = 2elAu. Using property 1 again yields,
elAu = e1Cu (3.26)
Thus if X=1 has multiplicity then elAu = elCu.
We prove that X=1 has multiplicity of two or more if
elAu = elCu by contradiction, more precisely by showing that
the following three statements contradict.
(i) X, = 1 with eI(A+B+C) = 0
(ii) Xj r 1 for j # 1 (multiplicity is 1) or
equivalently e2(A+B+C) + el(2A + B) # 0 Ve2 e Rn
(iii) elAu = elCu
Let Q = A+B+C and g = el(CA). Note that since el is
uniquely determined by Q, g is a known vector and
el(2A+B) = q.
Statement (ii) implies that
v = e2Q # el(CA) = Ve2 e Rn (3.27)
The vector v is any vector in the vector space, SQ,
defined by the rows of Q. Since rank(A+B+C) = n1, SQ is an
n1 dimensional space. Furthermore Q u = Q, which implies
v u = e2Qu = 0 or v I u (v is orthogonal to u), therefore
u I SQ.
From (iii) el(CA)u = g u = 0, that is g I u. Since
SQ I u and g I u, g is in the (n1) dimensional space SQ.
This contradicts with equation (3.27), since v = q for some
e2#O. This completes the proof.
Theorem 3.2:
If there is only one eigenvalue equal to one, then the
corresponding weight w1 is equal to zero.
Proof:
In order to prove Theorem 3.2, property 3 is used in
conjunction with equation (3.12), yielding
N+Z N.Z
lwj(kj)ke = wj(j)+ ejAu
Let li=l and kj1l for j=2,3,....,2n
N*Z N+Z
wleu + 2wj(kj)keju = wleAu + wj(kj)+lejAu
j=2 =2
(3.28)
Property 4 assures that the sums on the left and
righthand sides cancel, yielding,
wl(elAu elCu) = 0 (3.29)
Thus if elAu = elCu, then wl need not be zero. But
Theorem 3.1 assures that if elAu = elCu there are two or
more eigenvalues with value one. Therefore if there is only
one eigenvalue with value one then elAu # elCu, thus w1=0.
Theorem 3.3:
The eigenvalue, Xk = 1 has multiplicity m 2 3 if and
only if elAu = elCu and elAu = e2(CuAu) where el is the
genuine eigenvector and e2 is the generalized eigenvector.
Proof:
In order to prove Theorem 3.3, we first enumerate the
eigenvalues. Let X1 = 1 have multiplicity m 2 3 with the
corresponding eigenvectors el, e2, ...., em. The defining
equation for e3 is
lim el d2 + e2dL() + e3L() = 0 (3.30)
112 d2 d)
or,
elA + e2(2A + B) + e3(A + B + C) = 0. (3.31)
Postmultiplying equation (3.31) by u and using property
1 gives elAu + e2(2A + B)u. Using property 1 again yields,
elAu = e2(CA)u (3.32)
The proof elAu = elCu and elAu = e2(CA)u => =l has
multiplicity of three or more is obtained in a manner
similar to the one used in the proof of Theorem 1.
Theorem 3.4:
If there are two eigenvalues equal to one, then the
weight w2 corresponding to the generalized eigenvector is
equal to zero.
Proof:
Property 5 states that all the eigenvalues equal to one
form a single Jordan block. Let J1 be the Jordan block
corresponding to two eigenvalues equal to one, i.e.,
J1 = [0
Note that since we are dealing with left eigenvectors,
the off diagonal elements of J1 are below the diagonal.
Let l1=X2=1 and Xj#l for j=3,4,.... ,N+Z. Equation
(3.13) can be written explicitly as
1 0 0 0 0 el
1 1 0 0 0 e2
0 0 h3 0 O 3
gi = [wl, 2, 3, .. wN ] 0 0 0
0 0 0 0 4 e4
L 0 0 0 .N+ZJ LeN Z
N+Z
gi = (W1 + iw2)el + w22 + Z wjXjej (3.33)
3=3
In order to prove Theorem 3.4, property 3 is used in
conjunction with equation (3.33), yielding
(1 + iw2)el + w22 + ~wjjj Cu =
J=3
(w + (i + 1)w2)el +2 w22 ej Au (3.34)
3=3
Theorem 3.1 and property 4 are used to simplify
equation (3.34) yielding
w2elAu = w2e2(CA)u (3.35)
Thus if elAu = e2(CuAu) then w2 need not be zero. But
Theorem 3.3 assures that if elAu = elCu and elAu = e2(CuAu)
there are three or more eigenvalues with value one.
Therefore if there are only two eigenvalues with value one
then elAu # e2(CuAu), thus w2=0.
3.4 Observations
It is observed that, the numerical effort excluding
concerns of accuracy and precision in computing the steady
state probabilities is independent of the buffer size for
singlebuffer models. In practical terms, however, in order
to maintain numerical precision on a computer, the
computational effort would increase with buffer capacity, as
more digits of precision or variables with more range are
needed. Multiplication and division of doubleprecision
numbers take longer than operations on single precision
numbers. This was verified in all numerical work.
It follows from Theorem 3.1 that if only a single
eigenvalue with value one exists, it will receive a zero
weight. In the case of multiplicity, one of the eigenvalues
equal to one receives zero weight. The fact that the
balance equations hold in component form offers a way to
verify the numerical work in computing the eigenvectors.
It is clear that eigenvalueeigenvector pairs with
eigenvalues 1 l>1 have the most effect on the subvectors gi
for i+M and with eigenvalues 1I1<1 have the most effect on
the subvectors gi for i+0 Typically, eigenvalues with
large absolute values will only affect the solution at one
of the extremes of the solution vector, either corresponding
to emptybuffer states, or to fullbuffer states. This
property may be exploited in facilitating numerical
computations.
For general A, B, C submatrices with the property
(A+B+C)u=O, A20, C20 and B={bij} where bij= ,=
it can be shown that the eigenvalues of L(k) = X2A+XB+C need
not be real. However, in all the numerical experiments with
about 15 different production line models with various
assumptions on blocking, items transfer, and breakdown
mechanisms, similar to those given in Yeralan and Muth
(1987), it was observed that the eigenvalues were real. Not
having complex eigenvalues seems intuitive, since complex
eigenvalues are generally associated with oscillatory
behavior.
3.5 Comparison to Past Work
Gershwin and Berman (1981) presented a Markov process
model of a production line consisting of two unreliable
stations with random processing times and finite storage
buffers. Service, failure, and repair times for station i
are assumed to be exponential random variables with rates
Ii, Pi, ri; i=1,2, respectively.
In this section we present the matrix polynomial
solution for the same problem and we discuss the results
obtained by Gershwin and Berman (1981) and their
relationship to matrix polynomial solution procedure.
Gershwin and Berman first derived the difference
equations that describe the behavior of the internal states.
Then they propose a solution to these equations in the sum
ofproducts form as given in equation (3.36).
p(n,al,a2) = cXn1Y2 (3.36)
where p(n,al,a2) is the probability of a state, n is the
number of pieces in the buffer, and ai (i=1,2) is a binary
variable indicating whether station i is operational (ai=l)
or under repair (ai=0). Then they derived a nonlinear set
of equations which is obtained by inserting the proposed
solution into the internal equations. This set of
nonlinear equations are reduced to a fourth degree
polynomial in Y1. Each root of this fourth degree
polynomial determines a triplet, (X,Y1,Y2). Since the
internal equations are linear, any linear combination also
solves the internal equations. Therefore they conclude
that, for internal states
4
p(n,al, a2) = YcjX YlY22 (3.37)
where cj, j=1,2,3,4, are parameters determined by the
boundary equations.
In our matrix polynomial formulation, the model has
4M+8 states: (WOI), (DOI), (WiW), (WiD), (DiW), (DiD),
(BMW), (BMD) for i=0,1,2,...,M where M is the buffer
capacity. Station states W, D, I, and B refer to working,
brokendown, idle, and blocked stations respectively. Then
the steadystate probability vector, xi, is
7 i = ["WiW, WiD, 7DiW, lDiD o
The relationship between the steadystate probabilities
defined in our formulation and Gershwin and Berman's
formulation is written as
pWiW p(i,1,1)
"WiD p(il,0)
1 DDiW p(i,0,l)
"DiD p(i,0,0)
We first transformed the transition rate matrix for
this model into the standard form given in equation (2.1) by
elementary row operations. These operations constitute
expressing the steadystate probabilities of states (WOI),
(DOI), (BMW), and (BMD), KWOI, EDOI, 7BMW, and 7BMD as
P2r 1
KWOI = (WOW + DOW) "na = ( Pl IK + 42 "DkJW )
C1t rl
01 1
XBMW= ( WMW + xWMD) 7BMD= (P2 7BMW+ CLl7WMD)
JL2 r2
The corresponding submatrices A, B,C, BO, and BM are
four by four matrices as given below.
P2 0 0 0
0 0 0 0
0000
A=
0 0 P2 0
0000]
"u, 0 0 0
0 0 0 0
0 0 0 0
(l + i2 + P + P2) P2
B r2 (P, + p, + r
B =
0 r1
(l + P + PP2) P2
r (P + + r2)
r, +p2 0
0 r,
(2 + P + P2) P2
r2 + p, (pi + .LI + q)
BM =
0 00
0 r*1
PI
2) 0
(l2 + r' + r.
r2
P1
0
(.2 + r, + r2)
r1
Pi
0
(IL2 + r+ r2)
TZ
0
P1
) P2
(r, +r2)
01
P1
P2
(r1 + r2)
0
P1
P2
(r + r2)
At this point, before any calculations, it is known
that L(k) will have 2n = 8 eigenvalues. In addition, n
rank(C) = 2 of the eigenvalues will be 0 and nrank(A) = 2
of the eigenvalues will be infinitely large. The remaining
2n22 = 4 eigenvalues will be finite and nonzero.
Moreover, one of these eigenvalues will be one, leaving
three eigenvalueeigenvector pairs in the solution for the
interior states. In other words the determinant of the
associated matrix polynomial
L(1)= 12A + XB + C
is a sixth degree polynomial given by
IL(1)I = X2 (1)(%3 + a212 + al, + a0).
where the coefficients a0, al, and a2 can easily be computed
from the submatrices.
Comparing equation (3.36) to equation (3.7) shows that
the vector e can be written as
e = [Y1Y2, Y, Y2, 1] (3.38)
where Y1 and Y2 are the parameters defined by Gershwin and
Berman as given in equation (3.36) and X in this equation
corresponds to X, i.e. an eigenvalue of the matrix
polynomial L(X).
In Gershwin and Berman's paper, if the nonlinear set of
equations for internal states is reduced to a fourth degree
polynomial in X (instead of in Y1), the resulting fourth
degree polynomial is equivalent to IL(X) I/2. We know from
Theorem 3.1 that the eigenvalue with value one receives zero
weight for all production line models. Again Gershwin and
Berman show that this is true for their system too.
The eigenvector corresponding to the eigenvalue with
value one is easily computed as
el r1r2 rl r2 1 (3.39)
PJP2 P1 P2
Note that it is in the form given in equation (3.38).
Theorem 3.1 states that if elAu = elCu then eigenvalue one
has multiplicity. In this model this condition is
equivalent to
ri r2
1 2 2 (3.40)
Pi + r P2 + r2
In other words, as Gershwin and Berman noted, when the
standalone efficiencies of the stations are equal, one is a
double root of the matrix polynomial. In Gershwin and
Berman's solution procedure, this case must be treated
separately, whereas in Matrix Polynomial Solution Procedure,
it is in the framework defined by Jordan pair of the matrix
polynomial.
Now consider Theorem 3.3, Theorem 3.3 assures that the
eigenvalue X=1 has multiplicity m 2 3 if and only if
elAu = elCu and elAu = e2(CA)u. For this system, these
conditions are equivalent to
rl 2 3(Pl + P2 + rl + r2) + [LI + 52
r1 r 2 __
SPi + ri 2 P2 + r2 (2 ~)( + ) + 4 + 2 Llrl 6 P
rl r2 p 1r t 2r2 ri
(3.41)
We cannot offer any physical interpretation to the
latter part of equation (3.41). Gershwin and Berman do not
mention the possibility of eigenvalue X=l having
multiplicity m 3.
The matrix polynomial solution procedure is more
convenient for computerized solutions since its general
framework automatically accommodates multiplicity of any
eigenvalue.
In Gershwin and Berman's solution procedure, this
condition which let eigenvalue X=l have multiplicity m 2 3
must have been considered separately.
As it can be understood from the comparison of two
methodologies, the matrix polynomial solution procedure is a
general procedure of solving a system of difference
equations that describes the behavior of the internal states
of a Markov process model of a production line.
The versatility of the solution procedure is not only
in solving a specific system but solving a wide range of
onebuffer systems. The only effort required to solve a
production line model is in determining the submatrices A,
B, C, B0, and BM, given in the standard form of equation
(2.1).
3.6 An Example
We illustrate the solution procedure with a simple
numerical example. Consider a twostage production line
model with station breakdown and server variability as
sources of randomness. Service time, time to breakdown, and
time to repair are assumed to be exponentially distributed.
The system consists of a buffer with capacity M and two
workstations. Station 1 may breakdown with rate a and is
repaired with rate P. Station 2 never breaks down. The
service rate for station 1 and 2 are a and b respectively.
Station 1 Buffer Station 2
Twostage production line
Figure 3 
We define the system state by the triplet (i,j,k)
where i denotes the state of station 1, i=W( working), D
(down), or B (blocked); j is the number of items in the
buffer, j=0,1,...,M; and k denotes the state of station 2,
k=W (working) or S (starved). The states are ordered as
((WOS),(DOS),(WOW),(DOW), ... ,(WMW),(DMW),(BMW)}. The
transition rate matrix is given below.
(a + a)
b
0
a a
P 0
0 (a + a + b)
b P
. (a + a + b)
b
a
(3 + b)
0
Elementary row operations are performed to transform
the lowerright submatrix of P into the standard form given
in equation (2.1). The operation constitutes expressing the
a
steadystate probability of state (BMW) as times the
b
steadystate probability of the state (WMW). We obtain,
[b 0
0 b
B =(a + a +
0
Bo '(C + a)
B (a + b)
BM =
[a 0]
C = ]
0 0
b) a
(3 + b)]
a1
(3 + b)]
Part 1
Let, for example, a=4, 0=5, a=4, b=l and M=2. Then
[1 0 9 4
A = B =
0 1 5 61
B = 5 BM = 6
5 5 5 6
4 0
0 0
At this point, before any calculation we know that L(X) has
2n=4 eigenvalues. Property 2 states that another eigenvalue
is one. Moreover, nrank(C)=l of the remaining three
eigenvalues is zero. Since A is full rank, there does not
exist an eigenvalue that is infinitely large.
Note that the standalone production rate of
workstation 1 is ap/(a+3) = 20/9 and of workstation 2 is
b = 1. Therefore the production rate of the entire line is
bounded by min(20/9, 1) = 1.
Using the equation L(k) = 0, where L(k) is computed
from equation (3.10), we obtain the eigenvalues Xl=l,
%2=12, %3=2, %4=0. We substitute the eigenvalues into
equation (3.11) to obtain the eigenvectors el=(5,4),
e2=(3,2), e3=(l,1), e4=(0,1).
Thus
5 4 1 0 0 0
3 2 0 12 0 0
XF = JF =
XF 1 1 JF 0 0 2 0
0 1 0 0 0 0
The particular solution is sought by using the boundary
conditions. Substituting equation (3.13) into equation
(3.4), and equations (3.13) and (3.14) into equation (3.6),
we obtain
15 4 15 4
2 2 2 124 2 124
W = (0 o) and W = (0 0).
1 1 16 16
5 5 0 0
where the w is the weight vector to be calculated, i.e., w =
(wl, w2, w3, w4). We need to solve a simultaneous set of
four linear equations, i.e.,
5 4 15 4
2 2 2 12 2 12
W = (0 0 0 0),
1 1 16 16
5 5 0 0
yielding wl=0, w2=0.0000643, w3=l, w4 = 0.2, and hence,
using equations (3.13) and (3.14), go = (2595 3108),
gl = (5220 5160), g2 = (10800 10080), g3 = (25920 17280).
The nonnormalized probability for state (BMW) is obtained as
(4/1)25920=103680. We can verify the solution by checking
that gP = 0. The normalized steadystate probabilities are
calculated by normalizing the values obtained. The result
is
t=(0.014,0.017,0.028,0.028,0.059,0.055,0.141,0.094,0.564)
The corresponding production rate is 0.969 and the expected
buffer level is 1.7115.
Suppose now that the buffer capacity is increased to
three. Since the eigenvalues and eigenvectors do not
change, we proceed to compute the new weights. We have,
from equations (3.7) and (3.15),
15 4
2 125 2 *125
32 32
0 0) thus yielding,
5 4 15 4
2 2 2 125 2 125
w = (0 0 0 0).
1 1 32 32
5 5 0 0
With the new weights, subvectors gi are readily
available from equations (3.14) and (3.15). The steady
state probability vector is n then found by normalizing g.
x=(0.00696,0.0083,0.0139,0.0139,0.0280,0.0277,0.0580,
0.0541,0.139,0.093,0.556)
The production rate corresponding to increased buffer
size from is easily computed from this new steadystate
probability vector. The new productionrate is 0.9847 and
the expected buffer level is 2.6468.
The effect of increasing buffer size is shown in Table
3.1. The buffer size is denoted by M and the production
rate is denoted by r. Figure 3.2 shows the effect of
increasing buffer size. As mentioned, the production rate
is bounded by one.
49
Table 3.1
Effect of increasing buffer size
Buffer size
1
2
3
4
5
6
7
8
9
10
11
12
13
14
production rate
0.8651
0.9363
0.9690
0.9847
0.9924
0.9962
0.9981
0.9991
0.9995
0.9998
0.9999
0.9999
1.0000
1.0000
Production Rate versus Buffer Size
0.98
0.96
0.94
0.92
0.9
Pro 0.88
/ +t~
 I
0 2 4 6 8 10 12 14
Buffer size
Figure 3.2 Effect of increasing buffer size on the
production rate. (c=4, 0=5, a=4, b=l)
Part 2
Now, we will consider the balanced case for the same
production line. Theorem 3.1 assures that if elAu = elCu
then the line will be balanced and the eigenvalue with value
one will have multiplicity of two or more. For this model,
this condition is equivalent to
a = b
In other words if the standalone production rates of the
workstation 1 and 2 are equal then the eigenvalue with value
one will have multiplicity of two or more. In this example,
we will give a numeric example to illustrate the properties
given in the preceding sections.
Let a=2, 0=8, a=5, and b=4. Note that 5.(8/10) = 4,
that is, the line is balanced. Calculating the determinant
of the associated matrix polynomial L(1) yields four
eigenvalues: 1I = 0, X2 = 3.75, X3 = X4 = 1. Thus, as
Theorem 3.1 assures, the eigenvalue with value one has
multiplicity of two. Following the same steps as in Part 1
gives
0 1 0 0 0 0
0.8321 0.5547 0 3.75 0 0
Xp = Jp =
0.9701 0.2425 0 0 1 0
1 2.9701 0 0 1 1
Again the weight vector w is calculated by solving a
simultaneous set of four linear equations, yielding
wl=0.1198, w2=0.022, w3=0.9928, and w4=0. Note that the
weight associated with the generalized eigenvector, w4 is 0
as Theorem 3.4 assures. The normalized steadystate
probabilities are given by
n = (0.1531, 0.0570, 0.1539, 0.0375, 0.1569, 0.0355,
0.1681, 0.0280, 0.2102)
The production rate of this system is 3.1596 and the
expected buffer level is 1.0047. Note that since the line
is balanced the expected buffer level is very close to the
half of the buffer size.
CHAPTER 4
THE ONE STATIONONE BUFFERINSPECTION STATION MODEL
4.1 Introduction
The model developed in this chapter will serve two
purposes: to illustrate the general results and properties
uncovered in Chapter 3 on a specific example, and to
establish a flexible subsystem which will later be used as a
building block in decomposition techniques to study
multistation production systems.
Most of the recent decomposition methods are based on
decomposing the production line into twostation subsystems
(Svast'yonov, 1962, Gershwin, 1984, Choong and Gershwin,
1987, Dallery, David, and Xie, 1988, Gun and Makowski,
1989). In our approximation method, the production system
is decomposed into one stationone buffer subsystems. In
the literature, the decomposition methods presented by
Hillier and Boling, 1967, Takahashi, Miyahara and Hasegawa,
1980, Altiok, 1982, 1989, Perros and Altiok, 1986 are also
based on decomposing the production line into one station
one buffer subsystems.
First, the matrix polynomial formulation of the general
one stationone bufferinspection station model is given.
The closedform expressions of the performance measures are
derived in terms of the Jordan pair of the associated matrix
polynomial. Next, its relation to the sumofproducts
solution is discussed by examining the balance equations.
Finally, a passage time model in which the service
completion time distribution is obtained is used to model
the one bufferone stationinspection station system as a
M/G/1/N queueing system.
4.2 Model Description and Assumptions
Consider a production line with finite interstation
buffers fed by a stream of discrete parts. The input
process is homogenous Poisson process. When the buffer is
full, the input process is switched off, and when there is
space available in the buffer, the input process is started
again. For exponential interarrival times, because of the
memoryless property of the exponential distribution, this is
same as the lostarrivals model where the arrival process is
never stopped, but parts arriving while the buffer is full
are lost. Let the service time of each station be
exponentially distributed. The output of the each station
is inspected. Those parts that fail inspection are routed
to be reworked. In subsystem i, let Bi be the buffer
capacity, pi be the service rate of the station, and p be
the probability that a part passes inspection. Let us
approximate the input and output of workstation i as a
Poisson stream. Let ai be the rate at which parts are
brought into the buffer. Let 6i be the rate which parts
passing inspection are shipped out from the subsystem, and
let yi be the rate at which parts failing inspection are
sent to the rework station. In addition, let us approximate
station interference (blocking) by two parameters. Let bi
be the probability that the output of subsystem i into the
next (downstream) buffer is blocked. We approximate the time
subsystem i remains blocked by an exponentially distributed
random variable. Let Pi be the rate at which the blocking of
subsystem i is removed. Let us define a stochastic process
Y(t) to describe the blocking mechanism of the downstream
stations.
S 1 if the downstream buffer is full
Y( 0 if the downstream buffer is not full
A sample realization of the process Y(t) is given in
Figure 4.1.
Y(t)
1 
0 > t
TB
Time to blocking removal
Figure 4.1 A sample realization of the process Y(t)
The blocking probability b is the probability that
downstream buffer is full, i.e., Y(t) = 1.
b = P [Y(t) = 1]
Let TB be the continuous random variable which is the
length of the period which the output of the station is
blocked. Similarly, the blocking removal rate is defined as
P = 1 / E[TB].
For simplicity, the distribution of the blocking
removal time TB is approximated by an exponential
distribution with mean 1/0. This is a first order
approximation. The blocking removal time distribution may
be described more accurately by a phasetype distribution.
In this study, only the firstorder approximation is given,
but the extensions are possible.
Inspection is performed while the part is held at the
workstation. When the downstream buffer is full and a part
fails inspection, the workstation holding the defective part
is blocked. Figure 4.2 below depicts such a system.
to rework
from b p  
..... .... .... ... ........ .....
upstream
to downstrea
Workstation and work buffer
Figure 4.2 The one stationone bufferinspection
station model and its associated
parameters
Both unreliable and reliable workstations for this
production system are considered.
4.3 Model I Reliable Workstation
The state of the subsystem is expressed by a pair (i,j)
where i is the number of parts in the buffer and j is the
state of the station. There are 2M+3 possible states:
(0,I), and (i,W) and (i,B) for i=0,1,2,...,M where M is the
buffer capacity. Station states I, W, and B refer to idle,
working, and blocked stations, respectively. The submatrices
of each subsystem, provided that the states are ordered as
above, are given below. The subscripts of the rates are
omitted. A row operation is performed to bring the
submatrices to conform to those given in equation (2.1).
The row operation constitutes writing the steadystate
probability of state (0,I) as a function of the steadystate
probabilities of states (0,B) and (0,W).
A pp + pb 0p a 0
A = C =
P 0 0 a
B [(it + a) p[pb 1
0 (a + )J
B a Ipb pipb ppb
B0 = [ P a BM = [
S(c + 0) 0 P
(4.1)
The bar refers to the complement of the parameter,
e.g., p = 1 p. Since rank(A)=1, the model has three
finite eigenvalues. One of the eigenvalues is 1. Since
rank(C)=2 there is not an eigenvalue that is infinitely
large.
4.4 Model II Unreliable Workstation
The versatility of the matrix polynomial solution
technique is demonstrated by altering the model and the
architecture of the production system. In addition to the
model assumptions given in model I, now let the stations be
subject to breakdown. Assume that time to breakdown and time
to repair are exponentially distributed. Further, assume
that the station may breakdown only when working on a part.
Upon breakdown, the repair operation immediately begins.
Upon repair, the station continues with the servicing of the
part. Let pi and ni respectively be the rates at which the
station fails and is repaired. The model of each subsystem
now has 3M+4 states: (0,I), and (i,W), (i,B), (i,D) for
i=0,1,2,...,M where M is the buffer capacity. The
additional station state D refers to a brokendown station.
The corresponding submatrices, suppressing the indices of
the rates, are given below. The model with station
breakdown has four finite eigenvalues, none of which is
zero.
ppb + p 0 0 a 0
A= 0 0 0, C= 0 a 0 .
0 0 0 0 0 a
a PIpb .L pb (
BO = P a P 0
= o ,
a p P L ptpb (P j (pb
B = 0 B= 0 0 ,
0q o r e 11 0 T
(4.2)
Once the Jordan pair of the matrix polynomial is
obtained, i.e., the eigenvalues and the eigenvectors are
determined, performance measures such as the production
rate, the expected buffer level, the blocking probability
and the blocking removal rate can be calculated from the
steadystate probabilities as functions of the Jordan pair.
In the next section, the closedform expressions for the
production rate, and the expected buffer level are derived.
4.5 ClosedForm Expressions for the Performance Measures of
the Workstation
Let the nonnormalized steadystate probability vector
for the workstation be
=i = [gwi, gBit gDi]
where gji is the probability that the process in state
(i,j), i is the number of parts in the buffer, and j is W
(working), B (blocked), or D (down).
First the rate matrix is transformed to the standard
block tridiagonal form given in equation (2.1) by elementary
row operations. These operations consist of writing the
probability that the workstation is idle as
Pidle = _0.D (4.3)
wheir(p + pb) P '
where D = p + 0 .
a I
Then the sum of nonnormalized steadystate
probabilities, Sb is
M
Sb = gi g + Pidle (4.4)
1=0
where u = [1,1,1]T.
Assuming that there is only one eigenvalue equal to
one, let JF be the matrix obtained from JF by eliminating
the Jordan block corresponding to the eigenvalue equal to
one, i.e.,
JF = diag(J(X2), J(K3),*, J(Xp))
where Xj@l for j>l, and p is the number of different
eigenvalues of L(1). Note that (IJF)1 always exists.
Similarly, let XF be the matrix obtained from XF by
excluding the eigenvectors associated with the eigenvalue
one, and let W be the weight vector obtained from W by
excluding the weights associated with the eigenvalue equal
to 1.
Since the weight associated with the eigenvalue equals
to one is zero, the sum of the nonnormalized steadystate
probabilities for the case the eigenvalue one does not have
multiplicity, can be written in closedform by substituting
equations (3.14), (3.15) and (4.3) in equation (4.4) that
yields
M
Sb Z= gi + Pidle= (W JF)1 JF_ l)XF + WKXJ,)u + WXpD.
1=0
(4.5)
If the eigenvalue one has multiplicity, equation (4.5)
must include the effects of nonzero weights associated with
the eigenvalues equal to one.
Let the normalized steadystate probability vector be
7i = [(Wi'r Bi' "Di] = Sb1 9i
Let the normalized steadystate probability vector for
the steadystate probabilities that the station is working,
down and blocked be s = [sw, sB, sD].
M M
s = Sb Z 9i
1=0 1=0
W(I JF)I JF)XF + W,J.X,
s =    (4 .6)
( W(I F)I F FM1)4XF+ WJ.X,Xu + WXFD
Now, the production rate y is calculated as
y=sG where G=[lpb, 0, ]T (4.7)
Similarly the blocking probability b' is written as
b'= 7(M)u = Sb1g u = Sb[WJFMXF + WJX]u (4.8)
where Sb is given in (47). The blocking removal rate of the
subsystem P' is
P'= b'lrMH = (Sb')~WJFMXF+ W,,JX~)H (4.9)
where H=[p(p + pb), ]T.
Similar steps yield a closedform expression for the
expected buffer level E.
M M
E = 2 i7ju = S1b ig u
1=0 1=0
WJFI JF)2(I Y IFM((M + 1)I MJF))XF+ M WKJX,
(W(I JF)(I JF41)XF + WJXj~u + WXFD
(4.10)
Note that, if the Jordan pair of the matrix polynomial
L(1) can be written in closed form, i.e., the roots of
det(L(X)) can be expressed in closed form, then these
performance measures can also be expressed in closed form.
For this model, det(L(X)) is a fourth degree polynomial
and one is a root of this polynomial. Therefore it can be
factored as (Xl)Q(X) where Q(X) is a third degree
polynomial. Since the roots of a third degree polynomial
can be expressed in closed form, in principal, the
production rate of the one bufferone station model can be
written in closed form. An experiment with the symbolic
mathematics software Maple V R.2 gave the production rate in
closed form, but the expression was too lengthy to be of
practical use or to derive insights from.
The closedform expressions for the performance
measures of the rework station can easily be obtained from
the equations above by changing the vectors D, G, and H and
by using the Jordan pair of the matrix polynomial associated
with the rework station. Note that the rework station is
reliable and there is no inspection, i.e., p=0,. Therefore
D  (4.11)
[ a_T
and G = H = aD.
4.6 The Relationship Between the Matrix Polynomial Solution
Procedure and the Balance Equations
In this section, the balance equations for the one
stationone buffer with inspection model are derived. The
relationship between the sumofproducts solution and the
matrix polynomial solution is discussed. The state
transition diagram for the one bufferone stationinspection
station model is given in Figure 4.3.
The balance equations for the states at which the
buffer level is zero are
(xri(0) = P7B(O) + (pbp + Pp)7w(O) (4.12)
(a + P)rB(O) = (4pb)7w(0) (4.13)
(a + r)x7D(O) = (pxW(O) (4.14)
(a + + (p)xW(O) = TiD(O) + 17rB(l) + (l pb + Pp)w(1) + 7)i(0)
(4.15)
where 7j (i) is the steadystate probability that the process
is at the state (i,j) where i is the number of parts in the
buffer and j is W (working), B (Blocked), or D (down).
Figure 4.3 State transition diagram for the one station
one bufferinspection station model
The balance equations for the interior states, i.e.,
the states at which the buffer level is less than M and
greater than zero are
(a + I.. + (p)tW(n) =
7tnD(n) + anw(n 1) + pxB(n + 1) + (ptpb + Ptp)~(n + 1)
(4.16)
(a + P)B(n) = a7B(n 1) + ((.pb)xw(n) (4.17)
(a + rl)7D(n) = (pw(n) + acD(n 1). (4.18)
Similarly the balance equations for the states at which
the buffer is full are
0B(M) = a7tB(M 1) + (tpb)7t(M) (4.19)
r17D(M) = 9pW(M) + an (M 1) (4.20)
(a + P1 + p)xc(M) = ED(M) + caxM(M 1) (4.21)
The solution of the difference equations (4.12) to
(4.21) gives the steadystate probabilities. Gershwin and
Berman (1981) proposed a sumofproducts solution for a
similar set of difference equations for a twostation
unreliable production line with random processing times. In
the spirit of their approach, for our model, the following
solution for the interior steadystate probabilities is
proposed;
tw(n) = wkXW (4.22)
EB(n) = wkXYB (4.23)
7D(n) = wkXYD (4.24)
where w, X, YW, YB, YD are the scalars to be determined.
The scalars X, YW, YB, YD can be calculated by substituting
the equations (4.2224) into the interior equations.
Consequently the scalar weights, w, are calculated from the
boundary equations. If there exist more than one solution,
a linear combination of the possible solutions gives the solution.
Substituting the equations (4.224.24) into the equation
(4.17) yields,
_t pbhYw
YB = (4.25)
X(a + p) a
Furthermore, substituting the equations (4.224.24) into the
equation (4.18) yields,
YD = (4.26)
?(cL + Tr) a
Now, substituting the equation (4.24), (4.25), and (4.26)
into the equation (4.16) gives the following fourth degree
polynomial in X,
P(X) = A4X4 + A3k3 + A2 2 + Alj + AO
where
A4 = (a + )( + P)(pb + ip) + iLpbp]
A3 = aa + b4jh (a + h)(m + a + j)] bnpb + mpgc2a + h + bfa ampbb
A2 = aa + b(a + h)a jha + a2hnFp + mpg+ a2a + h + bfa(m + a + j)
A1 = a 20c + + + + (P
A0 = a3. (4.27)
It can be shown that P(X) is the determinant of the
matrix polynomial L(X). Therefore the roots of P(X) are the
eigenvalues of the matrix polynomial L(k). Since one is
always an eigenvalue of L(X) as Property 2 in Chapter 3
assures, one of the roots of P(X) is also one. The
equations 4.25 and 4.26 give a triplet {Yw, YB, YD) for each
root of P(%). It can be shown that the vector [Yw, YB, YD]
is the eigenvector of the matrix polynomial L(1) associated
with the eigenvalue X. The weights corresponding to
different values of X are then calculated by using the
boundary equations (4.12)(4.15) and (4.19)(4.21).
The matrix polynomial solution procedure is a general
procedure of solving a system of difference equations that
describe the internal states of a production system.
4.7 Equivalence of an Unreliable Station with a Reliable
Station
In this section, we investigate the relationship
between the M/G/1/N solution and the matrix polynomial
solution of the one bufferone stationinspection station
model. The general idea of treating unreliable stations by
substituting equivalent reliable stations has received some
attention (Gaver 1962, Muth 1979a, Alkaff 1985). We direct
our attention to this approach. As a result we obtain an
alternative method to calculate the performance measures of
one stationone buffer system. The relationship between the
M/G/1/N solution and the matrix polynomial formulation is
discussed.
The unreliable station with blocking can be replaced by
a reliable station with no blocking where service completion
time includes the effects of the service time, blocking time
and repair time. The reliable workstation with no blocking
is referred to as the equivalent workstation. The state
transition diagram of the service completion process is
depicted in Figure 4.4. The service completion time
distribution of the equivalent workstation is the same as
the first passage time distribution from the state at which
part arrives to the station, to the state at which the
service is completed.
part
arrives
B f
service completion
Figure 4.4 State transition diagram of the service
completion process
The rate matrix for the process depicted in Figure 4.4
 (P Pb p (p p b p + p1
0 P 0 P
71 0 71 0
0Q0 0 0
The states are ordered as W (Working), B (Blocked), D
(Down), and S (Service completion). Let matrix A be the
matrix obtained from Q by eliminating the row and column
corresponding to state S. The Laplace transform of the
cumulative density function of the first passage time to the
state S, F(s) is
F(s) = [1 0 0](sI A)1 1
1
s + ++ (p
F(s) [1 0 0] 0
s
r
 jab p
pbp
s+
0
9 1
0 1
s + i1
This equation yields
F(s) L(s + rj)(s(l b(l p)) + P)
F(s) 2
s(s + P)(S2 + s(7j + (l + p) + ir)
The Laplace transform of the probability density
function of the service completion time density, f(s), is
f(s) sF(s) = (s + )(s(l b(l p)) + ) (4.28)
(s + P)(s2 + s(r + (L + (p) + [an)
Therefore the service completion time distribution is a
phasetype distribution.
The expected service completion time, E[T], is derived
from equation (4.28) as
E[T] = lim df(s) oi + tb(l p)'r + P (4.29)
E[T] = lim (4.29)
s0 ds rtp
We define the service completion rate of the equivalent
station, is, as
1 1
E[T] b(l p) 1 + (4.30)
Note that if the blocking probability b is zero or the
blocking removal rate 0 goes to infinity, then the service
completion rate, ps is equal to the standalone production
rate of the station, i.e.,
S = 1 l if b = 0 or p >o.
(P +r
The following properties of the service completion
rate, ps, are derived from equation (4.30).
Properties of the Service Completion Rate:
1 ps is a concave increasing function of the
processing rate, p.
Hts =
23s
as
Prq(y(+ (p)
(PtI + pb(l p)rl + (PP)
2b(l p)202( + (p)
(Pri + pb(l p)rn + (pp)
1im A s = 0 IlimtS
mnk m
> 0
< 0
b(l p
b(l p)
2 p is a convex decreasing function of the blocking
probability, b.
01 (1 p)
(POi + pb(l p)rn + (P))
2 pt3 3(1 p)2
< 0
(On + pb(l p)rj + (po)
lim PSs =
b0 7r + (9
lim [S =
b 1l
n + + L (i p)
PrI + (PP + LTI(l P)
3 4s is a concave increasing function of the
blocking removal rate, P.
12r12b(1 p)
(OBn + Iib(l p)'i + (pP)
2t2 22(1 + (p)b(l p)
> 0
< 0
]ii s = 0o lbr IIs = 1
11 +
iS =
8b
ajLS
5J
s Ls
ap
8LS 
Fi 
~W
(13i + pb( P)' + (PP)b
> 0
Balanced Workstation. One stationone buffer
inspection system is said to be balanced, when the input
rate to the station a, is equal to the service completion
rate of the equivalent workstation, its, i.e.,
a = =S = (4.31)
P + tib(l p)r + (pp
Result. The matrix polynomial L(k) has two eigenvalues
equal to one, i.e., X=l has multiplicity of two, if and only
if the station is balanced.
Theorem 3.1 states that k=l has multiplicity of two, if
and only elAu = elCu. For this system, this is equivalent
to equation (4.31). That is the matrix polynomial L(k) has
two eigenvalues equal to one, i.e., %=l has multiplicity of
two, if and only if the station is balanced.
4.8 Modeling the One StationOne Buffer System as an M/G/1/N
Queue
When the unreliable workstation is replaced with a
reliable station, the service completion time of the
equivalent station is a phasetype distribution whose
Laplace transform is given in equation (4.28). Therefore
one stationone buffer model is equivalent to M/G/1/N queue
with switchoff arrivals where the service time distribution
is phase type, interarrival time distribution is exponential
and queue length is limited to N where N=M+1, and M is the
buffer size of the station. Note that because of the
memoryless property of the exponential distribution, for the
one bufferone station model, the switchoff arrival process
is equivalent to the lost arrival process.
All the performance measures can be derived from the
probability mass function of the number of customers in the
system. Let p'(n) be the probability that there are n
customers in the M/G/1/N queue. Similarly, let p(n) be the
probability that there are n customers in the M/G/1 queue.
The solution to the M/G/1/N queue can be derived from the
solution of M/G/1 queue by normalization (see Kleinrock,
1975), i.e.,
p(n)
p'(n) (n) n = ,...,N
1 pP(N)
(1 p)P(N)
p' (N)N (4.32)
1pP(N)
where P(N)= p(n) and p= is the traffic intensity and
is is the service completion rate given in equation (4.31).
The solution of the M/G/1 queue can be obtained by using
moment generating functions. The moment generating function
of p(n), GN(Z), is
G ( 1 = p)f (c(1 z))(1 z) (4.33)
G,(z)= (4.33)
f(a(l z)) z
where f(.) is the Laplace transform of the probability
density function of the service completion time as given in
equation (4.28).
Substituting (4.28) into (4.33) yield
GN(Z) = R(z)/P(z) (4.34)
where R(z) is a third degree and P(z) is a fourth degree
polynomial as given below.
R(z) = c RI(z) R2(z) R3(z)
where
c = [bp(l + a(l p)) + ap(n + p)]
R1(z) = z 1
R2(z) = az a r
R3(z) = za(l b(l p)) + ab(l p) a P
and
P(z) = p4 4 + p3z3 + p2z2 + plz + P0
where
p4 = a3
P3 = a2(3a + ri + p + (t + (p)
P2 = (a + P)(a + r)a (pIa + a 2(pb + ip) + (2a + r1 + P)a(p + a + (p)
pi = (a + P)[(pn (a + r)([ + a + )] (pb + ip)(2a + 1 + P)a auipbp
PO = (a + r)[(a + )(ppb + ip) + Lipbp] (435)
(4.35)
The probability function p(n) is obtained from (4.34)
by inverse transformation.
4.9 The Relationship Between the Matrix Polynomial Solution
and the M/G/1/N Solution
Comparing P(z) given in equation (4.35) to the
determinant of the matrix polynomial L(k) given in equation
(4.27) shows that the denominator of the moment generating
function GN(Z), P(z), is also the determinant of the matrix
polynomial L(X1). To understand this relationship, note
that P(z) is a fourth degree polynomial and R(z) is a third
degree polynomial. Therefore equation (4.35) can be
converted into partial fractions which has four terms in the
form of cj/(zzj), where cj and zj are scalars defined by
the model parameters, a, t, p, 7i, p, b, and P. The
corresponding term in the probability mass function is
Cj.(l/zj)n. The central role played by the latent values zj
is the same as those by the eigenvalues Xj of the matrix
polynomial. Due to the definition of the moment generating
function, it turns out that Xj = 1/zj. Similarly P(z) is
the determinant of the matrix polynomial L(1l).
Property 2 in Chapter 3 assures that X=1 is an
eigenvalue of L(X). To show that z=l is a root of P(z),
i.e., Z1 = 1, consider
lim P(z) = f(0) 1
z*l
Since f(s) is the Laplace transform of the proper density
function of the service time distribution, f(0)=1.
Therefore lim P(z) = 0.
l 1
Theorem 3.2 assures that the eigenvalue with value one
receives zero weight. In the M/G/1/N solution, since z=l is
also a root of the numerator R(z), they cancel each other.
In other words, it also receives zero weight.
In Appendix 4 an extension to the one bufferone
stationinspection station model is given. In this chapter,
it is assumed that the interarrival times of the parts
arriving to the station are exponentially distributed. In
Appendix 4, the one bufferone stationinspection station
model is modified to accommodate phasetype interarrival
time distribution. The matrix polynomial formulation of the
new model is given. An alternative solution procedure based
on the M/G/1/N solution is also presented.
CHAPTER 5
AN APPROXIMATION METHOD FOR MULTISTATION PRODUCTION LINES
WITH LIMITED BUFFER CAPACITY
5.1 Introduction
A prominent application of the methodology developed in
Chapters 3 and 4 is in constructing decomposition algorithms
for multistation production lines or systems. The latter,
i.e., a production system, refers to a collection of
stations interconnected as a general network.
We first discuss the properties of the decomposition
approximation for multistation production lines. The
uniqueness and convergence properties of the decomposition
approximation are studied. The insights gained in studying
the uniqueness and convergence properties are exploited to
modify the search algorithm to attain the highest possible
rate of convergence. Extensions to general production
systems are illustrated by a simple workcell with rework and
a production system with parallel stations.
First consider a production line partitioned into
segments containing a single buffer followed by a single
workstation as shown in Figure 5.1. Extensions to general
networks of workstations are given at the end of this
chapter.
Bi 1 bi
il Pi2
b .. 2
Figure 5.1 Decomposition of production lines
Each segment is evaluated by solving a subproblem.
Strictly speaking, the one bufferone stationinspection
station segments are not evaluated independently of the
system. The input process from the upstream station, and
the blocking process of the downstream station need to be
specified in the subproblem. These processes are unknown at
the beginning of the solution procedure.
In most of the production line models, it is assumed
that the first station is never starved. Assuming
exponentially distributed service times, time to breakdown
and time to repair, the production rate of the first station
can be given in closed form as a function of the station
parameters, the blocking probability b2, and the blocking
removal rate 02, of the second station. The state
transition diagram for the first station is depicted below.
The first station has only three states: W (Working), B
(Blocked), and D (down).
W B
D
Figure 5.2 Statetransition diagram of the first
station
The steadystate probabilities for the states that the
station is working, blocked, or down, KW, xB, or KD are
calculated from the statetransition diagram as
[ 1 1 p1102 P1
[EW XB t"D b= 1 2 1
1+p + _l 2 0 2 711
P2 ll
Then the production rate of the first station yl, is
calculated as
Y1 = tlTW b (5.1)
P1l 2 + + 1
P2 ll
where I, (p, q1 are the service rate, failure rate, and the
repair rate of the first station respectively. This
approach may be extended to cases where the distributions
associated with the first station are other than the
exponential, or where the items transfer mechanisms, station
breakdown and repair mechanisms are different.
Figure 5.3 depicts the decomposition of the first and
the second station of the production line.
l 1l Y1 a2! 2 72
S( 1f B22 2
b 21 b3
P21 3
Figure 5.3 Decomposition of the first station
Each subproblem gives the output rate yi, the
probability that it blocks the upstream subsystem bi, and
the rate at which its blocking effect is removed Pi. The
solution of one subsystem is used in the calculation of the
other subsystems, i.e., the blocking probability and the
blocking removal rate of a station are used in the
calculation of the upstream station. The process of
successive approximations is repeated until the solutions
converge. In this study, it is attempted to capture only
the firstorder effects of each subsystem. These effects
are then related to the upstream and downstream subsystems.
A wealth of extensions are possible but are not given in
this study.
Although the underlying ideas in this approximation
technique are simple, numerical experience proves it to be
very powerful. Furthermore, there are many conventions to
select the unknown parameters and use the result of one
subproblem to solve the next, which makes the approach a
general flexible framework rather than a rigid recipe.
5.2 Feedback Control Revisited
The underlying mechanism used in the decomposition
approximation technique is feedback control. In general
terms, the deviation of the system from its desired
performance is used as a restoring force to guide the system
back to its desired performance as shown in Figure 5.4.
reference
G x (t)
Figure 5.4 Feedback control
In our case, the desired performance corresponds to
satisfying the conservation of materials throughout the
production system. Again extensions are possible. When
scrapping is allowed, the output of a station plus the items
scrapped must be equal to the number of items the station
receives. That is the materials are not conserved within
the production system but are conserved in the larger system
consisting of the subsystem and the scrapping.
From an algorithmic model point of view, each subsystem
is conceptualized as a "black box" with multipleinputs and
multipleoutputs. Given the inputs, the outputs are
computed by the matrix polynomial solution procedure
developed in Chapters 3 and 4. The inputs to this system
are the rate at which parts are brought into the buffer, the
probability that the output of subsystem is blocked, and the
rate at which blocking of the output of the subsystem is
removed. The outputs are the rate at which parts passing
inspection are shipped out from the subsystem, the rate at
which parts failing inspection are routed to the rework
station, the probability that the subsystem blocks the input
stream, and the rate at which blocking of the subsystem is
removed.
Two mechanisms are used to interlink the subsystems.
One of them is to feedback some of the output variables of
the subsystem i to the subsystem i1. This ensures the
global information exchange between the subsystems. For
this purpose, the probability that output of the subsystem i
blocks the input stream bi, and the rate at which blocking
is removed Pi are used in the calculation of the subsystem
i1. This mechanism is analogous to statefeedback
mechanism in controlsystems theory. The second mechanism is
implemented by changing the input parameter according to the
deviation between the production rate of the station and a
reference. The only input variable that can be controlled
for this purpose is the rate at which the parts are brought
to the subsystem, ai. Note that the other inputs of the
subsystem, bi+1 and Pi+l, are the outputs of the downstream
subsystem.
b Subsystem > Yi
i+l Pi
Figure 5.5 "Black box" representation of each subsystem
If items are conserved, the output of each subsystem
must be the same, i.e.,
Y = Yi i = 1,2,...,N (5.2)
where N is the number of subsystems in the system.
An algorithmic solution procedure is employed to adjust
the input parameters of each subsystem based on the output
parameters of the previous iteration. More specifically,
the input rate to the subsystem i ai is adjusted such that
the conservation of flow equation (5.2) is satisfied when
the system converges. When the production rate of the
subsystem i is less than the production rate of the
subsystem il, the input rate to the subsystem i is
increased to reduce the deviation between the production
rates. Similarly, when production rate of the subsystem i
is greater than the production rate of the subsystem il,
ai is decreased. This constitutes the search procedure of
the algorithm. That is a successive approximation method in
the form
Lit+1 =cit K(yit Yilt) i = 2,3,...,N (5.3)
is used to calculate the input rates to the subsystems that
give the same production rate which is the production rate
of the system. In equation (5.3), superscript t is the
discrete iteration count, K is the proportionality constant,
and N is the number of stations. Note that as
Yi + Yi1, ait+1 ) it
It is empirically observed that such procedures
converge to the approximate solution. A more rigorous
discussion of the convergence properties is presented in
later sections.
5.3 Mathematical Programming Formulation
A formal model of the decomposition approximation may
be presented in several ways. As one alternative, consider
the following mathematical programming formulation:
minimize z
subject to
z < Oi < z
Oi = Yi Yil
Yi = Fi(ai, bi+1, Pi+l)
bi = Gi(ai, bi+i, Pi+l)
Pi = Hi(ai, bi+1, Pi+l)
ai = o0
bN+1 = 0, PN+1 = c
ai 2 0
(5.4)
(5.5)
(5.6)
(5.7)
(5.8)
(5.9)
(5.10)
(5.11)
(5.12)
The decision variables are the input rates ai,
i=2,3,...,N. The conservation of flow is expressed for a
series arrangement of workstations by the objective function
and the inequalities (5.5) and (5.6). The equality (5.10)
corresponds to the assumption that the first station is
never starved. Similarly, the qualities (5.11) correspond
to the assumption that the last station is never blocked.
The functions in equations (5.7), (5.8), and (5.9), Fi, Gi,
and Hi i=l,...,N, give the output variables, the production
rate yi, the blocking probability bi, and the blocking
i=2,3,...,N
i=2,3,...,N
i=l,2,...,N
i=l,2,...,N
i=l,2,...,N
removal rate pi as functions of the input variables, the
input rate ai, the blocking probability of the downstream
bi+1, and the blocking removal rate of the downstream Pi+l.
The objective function forces the conservation of flow
condition to be satisfied, i.e., it forces the maximum
difference between the output rates of adjacent workstations
to be minimized. As z+0, yi = Yii for i=2,3,...,N.
We refer to (5.4) as the global problem and finding yi,
bi, Pi by solving each subsystem for the current values of
the decision variables and system parameters as the
subproblem.
Although this formulation provides insights and helps
visualize the decomposition approach, it should be noted
that the functions Fi, Gi, Hi, i=l,...,N, in most cases, are
not available in closed form. They may be written as
algorithms following the matrix polynomial solution
procedure.
5.4 Properties of the Subproblem
It was mentioned that the functions Fi, Gi, and Hi
given in equations (5.7), (5.8), and (5.9) are not available
in closed form. However several properties of these
functions are known. The following seven postulates are
easily justified for a wide range of production lines.
These postulates can easily be verified for the one buffer
one station subsystem by using the properties of the service
completion rate given in Section 4.7 (Also see Figures 5.6,
5.7, and 5.8).
Postulate 5.1
(i) Fi, Gi, and Hi are realvalued, continuous,
nonnegative functions.
(ii) Fi and Gi are monotonically increasing functions
of the input rate ai, i.e., for ai > ai*
Fi(ai, bi+i, Pi+) > Fi(ci*, bi+1, Pi+l) and
Gi(ai, bi+l, Pi+l) > Gi((i*, bi+l, Pi+l)
(iii)Fi and Hi are monotonically decreasing functions
of the blocking probability of the downstream
bi+l, i.e., for bi+1 > bi+l*
Fi(ai, bi+1, Pi+) < Fi(ci, bi+l*, Oi+1) and
Hi(ai, bi+i, Pi+) < Hi(ci, bi+l*, Di+l)
(iv) Gi is a monotonically increasing function of the
blocking probability of the downstream bi+1, i.e.,
for bi+1 > bi+l*
Gi(ai, bi+1, Pi+1) > Gi(ai, bi+l*, Pi+1).
(v) Fi and Hi are monotonically increasing functions
of the blocking removal rate of the downstream
Di+l, i.e., for Pi+l > Pi+l*
Fi(ai, bi+1, Pi+l) > Fi(ci, bi+1, Pi+*) and
Hi(ai, bi+l, Di+l) > Hi(ai, bi+l, Pi+l*).
(vi) Gi is a monotonically decreasing function of the
blocking removal rate of the downstream Pi+l,
i.e., for Pi+1 > Pi+l*
Gi(ai, bi+l, Pi+l) < Gi(ai, bi+l, Pi+l*)*
(vii)Fi, Gi, and Hi are bounded above and below.
0 : Fi pi
0 < Gi < 1
0 5 Hi Pi+i
where pi and Pi+l are the standalone production
rates of the subsystems i and i+1 respectively.
The Figures 5.6, 5.7, and 5.8 depict the production
rate, the blocking probability, and the blocking removal
rate of the subsystem i as functions of the input rate to
the subsystem i, the blocking probability, and the blocking
removal rate of the downstream subsystem respectively for a
specific system.
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 5.6
The production rate, the blocking probability
and the blocking removal rate as functions
of the input rate. (One stationone buffer
inspection station model: ti=0.2, Qi=0.02,
ri=0.3, Mi=4, pi=0, bi+i=0.3, Bi+i=0.9)
 gamma i
  bi
  beta i
0.2 +
I I I I
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Pi+1
Figure 5.7 The production rate, the blocking
probability and the blocking removal rate
as functions of the blocking removal rate
of the downstream station, Pi+l. (same
parameters as.given in Figure 5.6 with ai=0.2)
gamma i
. bi
.. .. beta i
..................  . . ..... ....
_7
0.3
0 2 5 . .  ~ 
0.2 ..
0.15
0.15
0.1 gamma i
 bi
0.0   beta i
0.05
0I I I I I I I I I I
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
bi+l
Figure 5.8 The production rate, the blocking
probability and the blocking removal rate
as functions of the blocking probability
of the downstream station, bi+I. (same
model given in Figure 5.6)
Since it is assumed that the last station in the
production line is never blocked, the production rate, the
blocking probability, and the blocking removal rate of the
last station are only functions of the input rate aN, i.e.,
YN = FN(N)
bN = GN((aN)
PN = HN(aN). (5.13)
Therefore, the production rate, the blocking probability,
and the blocking removal rate of the subsystem N1 are
YN = FN1(aNl, FN(aN), GN(c(N))
bN1 = GNl(aN1, FN(aN), GN(clN))
ON1 = HNI((NI, FN(ON), GN(aN)).
That is, the production rate, the blocking probability,
and the blocking removal rate of the subsystem Nl are only
functions of the input rates to the subsystems Nl and N.
Similarly, the successive substitutions indicate that
the production rates yi may be related to aj j=i,i+l,...,N.
Let fi be the function that gives yi, that is
Yi = fi(ai, Ci+l,...,aN) (5.14)
In words, the production rate of the subsystem i is
given as a function of the input rates of its downstream
subsystems and its own input rate. The blocking of the
subsystem is described in terms of the downstream input
rates. The function fi inherits and combines several
properties of the functions Fi, Gi, and Hi (see also Figures
5.12, 5.13, 5.14, 5.15).
Result 5.1: Properties of the function fi
(i) fi is a realvalued, continuous, nonnegative
function (Postulate 5.1 (i)).
(ii) fi is a monotonically increasing function of ai,
i.e., for ai>ai*
fi(ai, / ,aj, ,aN) > fi(ai* ..'aj, ... (XN)
(Postulate 5.1 (i) and (ii)).
(iii) fi is a monotonically decreasing function of aj
for j>i, i.e., for aj>aj*,
fi(air*..aj,.** N) < fi(ai,..,j*,..,aN)
(Postulate 5.1 (iii) to (vi)).
(iv) fi is bounded above and below.
0 fi 5 Pi
where pi ; standalone production rate of the
subsystem i (Postulate 5.1 (vii)).
Let (i(ai) be the function given in equation (5.6)
that gives the difference of the production rates of the
station i and i1 as a function of the input rates to the
subsystems i1 to N, i.e.,
Oi(ai) = fi(ai,ai+l,..,N) fil(ailai,.,aN)
where Xi = (ai1,ai,ai+l,..,aN). The convergence and
uniqueness properties of the decomposition method depend on
the properties of the function i((ia).
The properties of Oi(gi) follow from the previous
results. (see also Figures 6.9, 5.16, and 5.17)
