• TABLE OF CONTENTS
HIDE
 Title Page
 Dedication
 Acknowledgement
 Table of Contents
 Abstract
 Introduction
 Past work
 Matrix polynomial formulation of...
 The one station-one buffer-inspection...
 An approximation method for multistation...
 Appendix
 Bibliography
 Biographical sketch
 Copyright














Title: decomposition method for multistation production systems
CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00090911/00001
 Material Information
Title: decomposition method for multistation production systems
Series Title: decomposition method for multistation production systems
Physical Description: Book
Creator: Tan, Baris,
 Record Information
Bibliographic ID: UF00090911
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 001981553
oclc - 31920320

Table of Contents
    Title Page
        Page i
    Dedication
        Page ii
    Acknowledgement
        Page iii
    Table of Contents
        Page iv
        Page v
    Abstract
        Page vi
        Page vii
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
    Past work
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
    Matrix polynomial formulation of production lines with limited buffer capacity
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
    The one station-one buffer-inspection station model
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
    An approximation method for multistation production lines with limited buffer capacity
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
    Appendix
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
    Bibliography
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
    Biographical sketch
        Page 173
        Page 174
        Page 175
    Copyright
        Copyright
Full Text





















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 STATION-ONE BUFFER-INSPECTION 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 Closed-Form 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 Station-One 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 STEADY-STATE
PROBABILITY VECTOR......................... 143

3 EXTENSION TO PHASE-TYPE 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

steady-state probabilities of a single-buffer production

line model. Perhaps the most significant result is that the

computational effort in finding the steady-state

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 station-one 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 on-line 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 two-pronged approach:

1. Develop an expedient and flexible methodology to

solve the subproblems consisting of a single-buffer

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 steady-state 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

single-buffer 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

work-in-process 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 work-in-

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 Quasi-Birth-Death Processes (QBDP). More specifically,

we develop a matrix polynomial formulation of the state-

homogeneous finite-state QBDPs. The properties of QBDPs are

then used to assist in the solution procedure for computing

the steady-state 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 station-one

buffer-inspection 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.

Closed-form 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 station-one











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 station-one

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 station-one 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 steady-state

probability vector from the matrix polynomial is presented.

The model presented in Chapter 4 is extended to phase-type

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 steady-state solutions,

computing the steady-state probabilities suffices. Most

approaches lead to Markovian models. Because of

computational complexities, the majority of production line

models are limited to two- and three-stage 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 phase-type 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 three-station 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), Avi-Itzak and Yadin (1965),

and Rao (1975) presented models of two- and three-station

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 three-station 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 three-station

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 one-step 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 Quasi-Birth-Death 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 steady-state 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 steady-state probabilities of their models.

These are the recursive solution procedure, the matrix

geometric procedure, and the scalar geometric procedure. A

general state non-homogeneous discrete flow model is given by

Yeralan and Muth. They present a very general view of the

two-station 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 three-station 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 inter-station 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 two-station 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 two-station

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 steady-state 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 zero-buffer

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 two-station

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 two-station

models is required. The method developed in Chapter 3 is

proposed as a general approach to evaluate a wide range of

single-buffer 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 steady-state

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 steady-state 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 steady-state probabilities. The model describes an

entire production system when there is only one buffer. It

is shown that for single-buffer 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 sum-of-products 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 one-step 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-

Birth-Death 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 state-homogeneous finite-state

QBDPs. The properties of QBDPs are then used to enhance a

solution procedure for computing the steady-state

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 steady-state probabilities. Most notably, once

the eigenvalues and eigenvectors of the associated matrix

polynomial are obtained, the steady-state 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 two-station

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 single-buffer

subsystem is required. The method developed in this chapter

is proposed as a general approach to evaluate a wide range

of single-buffer 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 steady-state

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 steady-state probability vector is unique

if R is irreducible. Any vector satisfying (3.1) but not

(3.2) is called a nonnormalized steady-state 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)

gi-iC + giB +gi+1A = 0 i = 1,2,...,M-1 (3.5)

gM-1C + 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 steady-state probabilities, our

approach exploits the fact that the interior equations are

repeated M-l 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,

ki-leC + 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 eigenvalue-eigenvector

pairs solve equation (3.5). From equations (3.9) and (3.10)

we get

eL(X) = 0 (3.11)

Any eigenvalue-eigenvector 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 2n-N-Z, 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 j-th 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 J-N+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(X-l), 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(k-1) corresponding to

X'=0. Equation (3.12) can be expressed in matrix form as



gi = WJFiXF i=0,1,2,...,M-1 (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
eigenvalue-eigenvector pair. Similarly, Wo is a row vector

of size 2n-Z-N. 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 i-l. 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,...,M-1 (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 k-l, i.e.,

gk-1Cu = gkAu (3.20)











Again, using property 1 and equation (3.5) for i=k, we

obtain

gk-iC1 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 eigenvalue-eigenvector 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

m-1 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)
1-L 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(C-A). 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(C-A) = 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) = n-1, SQ is an

n-1 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(C-A)u = g u = 0, that is g I u. Since

SQ I u and g I u, g is in the (n-1) 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

right-hand 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(Cu-Au) 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(C-A)u (3.32)

The proof elAu = elCu and elAu = e2(C-A)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(C-A)u (3.35)
Thus if elAu = e2(Cu-Au) then w2 need not be zero. But

Theorem 3.3 assures that if elAu = elCu and elAu = e2(Cu-Au)

there are three or more eigenvalues with value one.

Therefore if there are only two eigenvalues with value one

then elAu # e2(Cu-Au), 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

single-buffer 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 double-precision

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 eigenvalue-eigenvector 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 empty-buffer states, or to full-buffer 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-

of-products 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,

broken-down, idle, and blocked stations respectively. Then

the steady-state probability vector, xi, is



7 i = ["WiW, WiD, 7DiW, lDiD o

The relationship between the steady-state 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 steady-state 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 n-rank(A) = 2

of the eigenvalues will be infinitely large. The remaining

2n-2-2 = 4 eigenvalues will be finite and nonzero.











Moreover, one of these eigenvalues will be one, leaving

three eigenvalue-eigenvector 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

stand-alone 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(C-A)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

one-buffer 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 two-stage 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


Two-stage 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 lower-right submatrix of P into the standard form given

in equation (2.1). The operation constitutes expressing the

a
steady-state probability of state (BMW) as times the
b

steady-state 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, n-rank(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 stand-alone 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 steady-state 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 steady-state

probability vector. The new production-rate 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 stand-alone 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 steady-state

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 STATION-ONE BUFFER-INSPECTION 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 two-station 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 station-one 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 station-one buffer-inspection station model is given.

The closed-form expressions of the performance measures are

derived in terms of the Jordan pair of the associated matrix

polynomial. Next, its relation to the sum-of-products

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 buffer-one station-inspection 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 lost-arrivals 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 phase-type distribution.

In this study, only the first-order 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 station-one buffer-inspection
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 steady-state

probability of state (0,I) as a function of the steady-state

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 broken-down 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

steady-state probabilities as functions of the Jordan pair.

In the next section, the closed-form expressions for the

production rate, and the expected buffer level are derived.



4.5 Closed-Form Expressions for the Performance Measures of
the Workstation


Let the nonnormalized steady-state 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 steady-state

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 (I-JF)-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 steady-state

probabilities for the case the eigenvalue one does not have

multiplicity, can be written in closed-form 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 steady-state probability vector be



7i = [(Wi'r Bi' "Di] = Sb-1 9i
Let the normalized steady-state probability vector for

the steady-state 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 closed-form 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 (X-l)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 buffer-one 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 closed-form 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

station-one buffer with inspection model are derived. The

relationship between the sum-of-products solution and the

matrix polynomial solution is discussed. The state-

transition diagram for the one buffer-one station-inspection

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 steady-state 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 buffer-inspection 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) = (p-w(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 steady-state probabilities. Gershwin and

Berman (1981) proposed a sum-of-products solution for a

similar set of difference equations for a two-station

unreliable production line with random processing times. In

the spirit of their approach, for our model, the following

solution for the interior steady-state 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.22-24) 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.22-4.24) into the equation

(4.17) yields,



_t pbhYw
YB = (4.25)
X(a + p)- a


Furthermore, substituting the equations (4.22-4.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 buffer-one station-inspection 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 station-one 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

phase-type 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)
s-0 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 stand-alone 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
mn-k 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 =
b-0 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


]i-i s = 0o lbr I-Is = 1
11 +


iS =
8b

ajLS
5J


s Ls
ap

8LS -
Fi -
~W


(13i + pb( P)' + (PP)b


> 0











Balanced Workstation. One station-one 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 Station-One 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 phase-type distribution whose

Laplace transform is given in equation (4.28). Therefore

one station-one buffer model is equivalent to M/G/1/N queue

with switch-off 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 buffer-one station model, the switch-off 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)
1-pP(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(X-1). 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/(z-zj), 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(1-l).











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 buffer-one

station-inspection 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 buffer-one station-inspection station

model is modified to accommodate phase-type 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 buffer-one station-inspection

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 State-transition diagram of the first
station


The steady-state probabilities for the states that the

station is working, blocked, or down, KW, xB, or KD are

calculated from the state-transition 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 first-order 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 multiple-inputs and

multiple-outputs. 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 i-1. 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










i-1. This mechanism is analogous to state-feedback

mechanism in control-systems 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 i-l, 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 i-l,

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 Yi-lt) 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 -+ Yi-1, 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 Yi-l

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 = Yi-i 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 real-valued, 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 stand-alone 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 station-one 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 N-1 are











YN- = FN-1(aN-l, FN(aN), GN(c(N))
bN-1 = GN-l(aN-1, FN(aN), GN(clN))

ON-1 = HN-I((N-I, FN(ON), GN(aN)).


That is, the production rate, the blocking probability,

and the blocking removal rate of the subsystem N-l are only

functions of the input rates to the subsystems N-l 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 real-valued, 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(a-i,..,j*,..,aN)
(Postulate 5.1 (iii) to (vi)).

(iv) fi is bounded above and below.

0 fi 5 Pi

where pi ; stand-alone 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 i-1 as a function of the input rates to the

subsystems i-1 to N, i.e.,



Oi(ai) = fi(ai,ai+l,.-.,N) fi-l(ai-lai,-.,aN)


where Xi = (ai-1,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)




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs