Analysis of Bernoulli Lines
Motivation: The problem of performance analysis of production systems consists of investigating their performance measures, e.g., production rate (PR), work-in-process (WIPi), blockages (BLi) and starvations (STi), as functions of machine and buffer parameters. In principle, this investigation can be carried out using computer simulations. However, the simulation approach has two drawbacks. First, it is not conducive to analysis of fundamental properties of production systems, e.g., relationships between system parameters and performance measures. Second, simulations require a relatively lengthy and costly process of developing a computer model and its multiple runs for statistical evaluation of the performance measures. In some cases, especially when numerous "what if" scenarios must be analyzed, this approach may become prohibitively expensive and slow. This problem is exacerbated by the exponential explosion of the dimensionality of the system as a function of buffer capacity. Indeed, even in the Bernoulli reliability case (i.e., when the machines are memoryless), a serial line with, say, 11 machines and buffers of capacity 9, has 1010 states, which is overwhelming for simulations. Therefore, a quick, easy and revealing method for production systems analysis, based on formulas, rather than on simulations, is of importance. The purpose of this chapter is to present such a method for serial production lines with Bernoulli machines.
Overview: The analytical approach to calculating PR, WIPi, BLi and STi is based on the mathematical models of production systems discussed in Chapter 3. Due to the complex nature of interactions among the machines, closed-form expressions for their performance measures are all but impossible to derive, except for the case of systems with two machines. Therefore, the approach, developed here, is based on a two-stage procedure: First, analytical formulas for performance analysis of two-machine lines are derived and, second, an aggregation procedure is developed, which reduces longer systems to a set of coupled two-machine lines and recursively evaluates their performance characteristics. This approach, illustrated in Figure 4.1, leads to sufficiently accurate estimates of the performance measures PR, WIPi, BLi and STi.
Figure 4.1: Block diagram of the analysis procedure
In addition, based on the analytical results obtained, this chapter investigates several system-theoretic properties, which provide qualitative insights into the behavior of serial lines.
4.1 Two-machine Lines
4.1.1 Mathematical description
System: The production system considered here is shown in Figure 4.2. The time is slotted with slot duration equal to the cycle time of the machines, and machines m1 and m2 are up during each time slot with probability p1 and p2, respectively. The buffer is of capacity N < ∞.
Figure 4.2: Two-machine Bernoulli production line
States of the system: Since the Bernoulli machines are memoryless, the states of the system coincide with the states of the buffer, i.e., the state space consists of N + 1 points: 0, 1,..., N.
Conventions: The following conventions are used to define the system at hand:
- (a) Blocked before service.
- (b) The first machine is never starved; the last machine is never blocked.
- (c) The status of the machines is determined at the beginning and the state of the buffer at the end of each time slot.
- (d) Each machine status is determined independently from the other.
- (e) Time-dependent failures.
4.1. TWO-MACHINE LINES
State transition diagram: Since the buffer occupancy can change in each time slot at most by one part, the state transition diagram is "linear", as shown in Figure 4.3.
Figure 4.3: Transition diagram of two-machine Bernoulli production line
Transition probabilities: In the expressions that follow, the events {mi is up during the time slot n+1} and {mi is down during the time slot n+1} are denoted for the sake of brevity as {mi up} and {mi down}. The buffer occupancy at slot n is denoted as h(n). In these notations, the transition probabilities are:
Using the formulas for the probability of the union of mutually exclusive events (2.1) and for the probability of the intersection of independent events (2.3), these transition probabilities can be calculated as follows:
Dynamics of the system: Since the transition probabilities (4.1) are constant, the system under consideration is a Markov chain. Let Pi(n) be the probability of state i, i = 0, 1,..., N, at time n. Then, as it is shown in Subsection 2.3.3, the evolution of Pi(n) can be described by the following constrained linear dynamical system:
Statics of the system: The steady state of the system at hand is described by the balance equations
Their solution provides a complete characterization of the system behavior. This is accomplished next.
4.1.2 Steady state probabilities
Since the states are communicating and there are "self-loops" (see Figure 4.3), the Markov chain (4.1) is ergodic. Therefore, there exists a unique stationary probability mass function. Taking into account (4.1), the balance equations can be re-written as
These equations can be solved consecutively in terms of P0. Indeed, from the first equation in (4.6) we have
It is convenient to re-write this expression as
and denote the second factor in the right hand side as
Then, from the second equation of (4.6),
and so leading to
where, for the sake of brevity, α(p1,p2) is denoted as α. Thus,
To complete the calculation, the expression for P0 must be derived. Using (4.5) and (4.8), we obtain
Thus
In the special case of identical machines, i.e., when p1 = p2 =: p,
and, therefore,
implying that
An illustration of the pmf (4.10), (4.11) is given in Figure 4.4 for p = 0.95 and p = 0.55 with N = 5. Clearly, when p is close to 1 (which is the practical case),
(a) p1 = p2 =0.95, N = 5 (b) p1 = p2 = 0.55, N =5
Figure 4.4: Stationary pmf of buffer occupancy in two-machine lines with identical Bernoulli machines
Intuitively, one would expect that the pmf of the buffer occupancy in the case p1 = p2 would be symmetric in the sense that P0 = PN . The fact that it is not is due to the blocked before service assumption, i.e., m1 is not blocked even if h = N but m2 is up. If this assumption is changed to "m1 is blocked if h = N" (i.e., irrespective of the status of m2), it is possible to show that the buffer occupancy is indeed symmetric, i.e., P0 = PN (see Problem 4.5). We follow, however, the original assumption since it is closer to the reality of the factory floor (Note that in the case of continuous time systems, discussed in Part III, both assumptions lead to the same conclusion – the pdf of buffer occupancy for the case of identical machines is symmetric.)
When p1 ≠ p2, substituting (4.7) into (4.9) we obtain:
Summing up the geometric series in the denominator, this can be re-written as
i.e.,
Substituting (4.7) into the first term in the denominator of the above expression, after simplification we finally obtain:
Thus, equations (4.8), (4.12) describe the steady state pmf of buffer occupancy in two-machine lines with non-identical Bernoulli machines. This pmf is illustrated in Figure 4.5 for the following serial lines:
Clearly, L1 and L2 are a reverse of each other, in the sense that the first (respectively, second) machine of L1 is the second (respectively, first) machine of L2, while the buffer remains the same. Similarly, L3 and L4 are also reverse of each other. The pmf's of Figure 4.5 clearly show that the buffer tends to be empty (respectively, full) if p1 − p2 < 0 (respectively, p1 − p2 > 0), and this phenomenon becomes more pronounced when | p1 − p2 | is large.
Figure 4.5: Stationary pmf of buffer occupancy in two-machine lines with nonidentical Bernoulli machines
The functions in the right hand side of (4.10) and (4.12) play an important role in the subsequent analyses. Similar functions appear in all other cases of serial lines (e.g., when the machines are exponential). Therefore, we introduce a special notation:
The properties of this function are as follows:
Lemma 4.1 Function Q(x, y, N), where 0 < x < 1, 0 < y < 1, and N ∈ {1, 2,...}, takes values on (0, 1) and is
- strictly decreasing in x,
- strictly increasing in y,
- strictly decreasing in N.
Proof: See Section 20.1.
The properties of Q, as indicated in Section 4.2, ensure the convergence of the aggregation procedure that is used for analysis of serial lines with more than two machines.
To conclude this subsection, we re-write the expression for the probability of the buffer being full (i.e., PN) in a form, which is more convenient for the performance measure formulas described in Subsection 4.1.3.
From (4.8) and (4.12) we obtain
Dividing both numerator and denominator by αN and replacing [p2(1 − p1)]/[p1(1 − p2)] with 1/α, after simplification, we obtain
Taking into account that, as it follows from (4.7),
this, finally, can be re-written as
4.1.3 Formulas for the performance measures
Using the steady state probabilities derived above and function Q defined by (4.14), the performance measures PR, WIP, BL1, and ST2 can be expressed as shown below.
Production rate: Keeping in mind conventions (a)-(e) listed at the beginning of this section and formula (2.3) for the probability of the intersection of independent events, PR can be represented as
Alternatively, PR can be expressed as
While the probability of the first event in the right hand side of (4.18) is p1, the probability of the second event is
Therefore,
where, as it follows from (4.11) and (4.16),
This implies that
and, therefore,
Thus, PR can be expressed in two equivalent ways:
These are important expressions: along with their direct value as PR of two-machine lines, they are the basis of the aggregation procedure for the analysis of M > 2-machine lines (see Section 4.2).
Work-in-process: The average value of pmf's (4.10), (4.11) and (4.8), (4.12), i.e., WIP , is given by
For p1 = p2 = p, this leads to
For p1 ≠ p2, after some algebraic manipulations, this can be reduced to
Therefore,
Blockages and starvations: As it follows from the definitions of Subsection 3.6.3,
Taking into account that the probabilities of the buffer being full, PN, and being empty, P0, are given by (4.16) and (4.14), respectively, these relationships can be expressed as
Clearly, these expressions are in agreement with formulas (4.19) for PR, which now can be understood as
4.1.4 Asymptotic properties
Expressions (4.19), (4.23) and (4.24) reveal the following properties of the performance measures as N → ∞:
Theorem 4.1 In a two-machine line with Bernoulli machines defined by conventions (a)-(e),
Proof: See Section 20.1.
The last expression in (4.26) implies that for N sufficiently large and p1 = p2, the buffer is, on the average, half full.
Theorem 4.1 is illustrated in Figure 4.6 for three serial lines defined by
กก
Figure 4.6: Performance measures of two-machine Bernoulli lines as functions of buffer capacity
Clearly, PR is monotonically increasing but with a decreasing rate, while WIP increases linearly (if p1 ≥ p2). This implies that there is nothing to be gained from having a buffer of capacity larger than 5. The issue of buffer capacity selection is explored in details in Chapter 6.
4.2 M > 2-machine Lines
4.2.1 Mathematical description and approach
System: The production system considered here is shown in Figure 4.7. The time is slotted, and each machine mi, i =1,..., M, is up during a time slot with probability pi and down with probability 1 − pi, i =1,..., M. The capacity of buffer i is Ni < ∞, i =1,..., M − 1.
Figure 4.7: M-machine Bernoulli production line
Conventions: Similar to two-machine lines, the following conventions are used:
- (a) Blocked before service.
- (b) Machine m1 is never starved for parts; machine mM is never blocked by subsequent operations.
- (c) The status of the machines is determined at the beginning and the state of the buffers at the end of each time slot.
- (d) Machines' status is determined independently from each other.
- (e) Time-dependent failures.
Throughout this book, we refer to (a)-(e) intermittently as either conventions or assumptions.
States of the system: The Bernoulli machines are memoryless, and, therefore, the states of the system coincide with the states of the buffers. Since the i-th buffer has Ni + 1 states, the system has (N1 + 1)(N2 + 1) ··· (NM−1 + 1) states. For example, if Ni = 9 for all i and M = 23, the number of states is 1022, which equals the number of molecules in a cubic centimeter of gas under normal pressure and temperature! Clearly, a direct analysis of such a large system is impossible (and, perhaps, unnecessary). Therefore, a simplification is in order. We use for this purpose an aggregation approach.
กก
Idea of the aggregation: Consider the M-machine line and aggregate the last two machines,
mM−1, and mM , into a single Bernoulli machine denoted as
, where b stands for
backward aggregation (see Figure 4.8(a)). The Bernoulli parameter,
, of this machine is assigned as the production rate
of the aggregated two-machine line, calculated using the
first expression of (4.19). Next, aggregate this machine, i.e.,
, with mM−2 and obtain another aggregated machine
. Continue this process until all the machines
are aggregated into 
, which completes the
backward phase of the aggregation procedure.
It turns out that the Bernoulli parameter of
might be quite different from the production rate of the
M-machine line under consideration. To remedy this problem, we introduce the
forward phase of the aggregation procedure defined as follows: Aggregate the
first machine, m1 with the aggregated version of the rest of the line, i.e., with
. This results in the aggregated machine, denoted as
, where f stands for forward aggregation (see Figure 4.8(b)). The Bernoulli
parameter, 
, is assigned as the production rate of the aggregated two-machine
line, calculated using the second expression of (4.19). Next, aggregate
with
, resulting in
and so on until all the machines are aggregated into
which completes the forward phase of the procedure. Again, the Bernoulli
parameter of 
may be quite different from the actual production rate of
the M-machine system.
To alleviate this discrepancy, we iterate between the backward and forward aggregations. In other words, view the above backward and forward aggregations as the
first iteration of the aggregation, i.e., as s = 1. At the second iteration,
s = 2, 
is aggregated
with mM to result in
for the second iteration, which is then aggregated with
and so on until the second
iteration of the backward aggregation is complete. Next, the second iteration of the forward aggregation is carried out, followed by the third iteration of the backward aggregation and so on, i.e.,
s =3, 4,....
We show below that the steady states of this recursive procedure lead to relatively accurate estimates of the performance measures for the M-machine line. We also indicate that the convergence to these steady states is quite rapid, typically requiring less than 10 iterations, which implies that the performance measures are evaluated within a fraction of a second.
Figure 4.8: Illustration of the aggregation procedure
4.2.2 Aggregation procedure and its properties
The mathematical representation of the recursive aggregation procedure described above and its properties are given below.
Recursive Aggregation Procedure 4.1:
with initial conditions
and boundary conditions
As before,
where
These equations are solved as follows: With i =
M − 1, using the initial condition
and the boundary condition
, solve the first equation of (4.30) to obtain
; then solve it with i =
M − 2 to obtain
, and so on, until
is obtained. Next, solve the second equation of (4.30) with
i = 2 to obtain 
; then solve it with
i = 3 to obtain 
and so on, until
is obtained. This completes the first iteration of the aggregation procedure. For the second, third, ..., iterations, this process is repeated anew
using 
, ..., respectively, in the first equation of (4.30).
กก
Example 4.1 Consider a three-machine serial line with p = [0.9, 0.9, 0.9] and N = [2, 2]. Then,
All subsequent iterations are carried out similarly.
Convergence: Clearly, (4.30) is an (M − 1)-dimensional dynamical system, which iterates pi's and Ni's and results in two sequences of numbers
defined on the interval (0,1). The properties of these sequences and their physical meaning are described below.
Based on (4.30) and Lemma 4.1, the following can be proved:
Proof: See Section 20.1.
Figure 4.9 illustrates the behavior of
, s = 0, 1, 2, ...,
for the following lines:
Obviously, L1 represents lines with identical machines, while
L2, L3 and L4 illustrate lines where
pi's are allocated according to an increasing, inverted bowl, and bowl patterns, respectively. As one can see,
indeed exhibit the properties established in Theorem 4.2. In addition, Figure 4.9 shows that convergence to the limits is quite fast: 2-4 iterations of the aggregation procedure for the uniform, ramp, and bowl allocations and about 15 iterations for the inverted bowl.
Figure 4.9: Illustration of the dynamics of
Interpretation of
: From the point of view of each buffer
bi, i =1,..., M − 1, the upstream of the line is represented by the
"virtual" Bernoulli machine
defined by the parameter
. Similarly, the downstream
is represented by the "virtual" machine
defined by
. In addition, the whole line can be represented by
. Thus, the M-machine
line can be represented as shown in Figure 4.10. Clearly, all the performance measures of the two-machine lines included in this
figure can be calculated using the formulas of Subsection 4.1.3.
Figure 4.10: Equivalent representations of Bernoulli M > 2-machine line through the aggregated machines.
4.2.3 Formulas for the performance measures
Using the equivalent representations of Figure 4.10 and the limits (4.35), estimates of the performance measures of the M > 2-machine line are introduced below.
Production rate: Based on Figure 4.10 and expression (4.19), the estimate, PR, of the production rate is defined as
Work-in-process: Using the two-machine representation of the
M > 2-machine line (Figure 4.10) and expression (4.23), the estimate,
, of the steady state occupancy of buffer
i is defined as
Obviously, the estimate of the total WIP is
Blockages and starvations: Since these probabilities must evaluate blockages and starvations of the real, rather then aggregated, machines, taking into account expressions (4.24), the estimates of these performance measures,
and
are introduced as follows:
These expressions give another interpretation of
. Indeed, using the steady states of the aggregation procedure (4.30), and expressions (4.39), (4.40), we obtain
PSE Toolbox: The recursive procedure (4.30) and performance measure estimates (4.36)-(4.40) are implemented in the Performance Analysis function of the toolbox. For a description and illustration of this tool, see Subsection 19.3.1.
4.2.4 Asymptotic properties of M > 2-machine lines
Formulas (4.36)-(4.40) can be used to investigate asymptotic properties of Bernoulli lines as N → ∞. This is carried out using the following eight serial lines:
The reasons for selecting these particular lines are as follows: Line 1 illustrates the behavior of systems with identical machines. Lines 2 and 3 represent systems with increasing and decreasing machine efficiency, respectively; clearly, L3 is the reverse of L2. Lines 4 and 5 illustrate systems with machine efficiency allocated according to a bowl and an inverted bowl patterns, respectively. Lines 6 and 7 exemplify systems with "oscillating" machine efficiency allocation. Finally, Line 8 is selected to illustrate the case of a good machine surrounded by low efficiency ones.
Figures 4.11-4.18 show the performance measures of these lines as a function of N. Based on this information, the following can be concluded: As N → ∞,
In addition, Figures 4.11 - 4.18 indicate that since
as a function of
N exhibits
saturation, while tends to infinity (except for the case of increasing machine efficiency), there is no reason for having large
N, typically, above 8 -10. A detailed discussion on selecting
N is included in Chapter 6.
Figure 4.11: Performance of Line L1 (identical machines)
Figure 4.12: Performance of Line L2 (decreasing machine efficiency)
Figure 4.13: Performance of Line L3 (increasing machine efficiency)
Figure 4.14: Performance of Line L4 (bowl machine efficiency allocation)
Figure 4.15: Performance of Line L5 (inverted bowl machine efficiency allocation)
Figure 4.16: Performance of Line L6 ("oscillating" machine efficiency allocation)
Figure 4.17: Performance of Line L7 ("oscillating" machine efficiency allocation)
Figure 4.18: Performance of Line L8 (good machine surrounded by low efficiencies ones)
4.2.5 Accuracy of the estimates
Clearly, formulas (4.36)-(4.40) are just estimates of the true values of the performance measures. Therefore, an analysis of their accuracy is of importance. In this subsection, such an analysis is carried out using both analytical and numerical tools.
Analytical investigation: Consider the
joint probability 
that buffers
i, i + 1, ..., j, contain hi,
hi+1, ..., hj parts, respectively. In general, this joint probability is not close to the product of its marginals, i.e.,
where Pi[hi] is the probability that the i-th buffer contains hi parts. However, it turns out that for certain values of hi, ..., hj, related to blockages and starvations, these probabilities are indeed close to each other. Specifically, define
and
Extensive numerical simulations show that δ is practically always small. An illustration is given in Table 4.1 for several four-machine lines with Ni = 3, i =1, 2, 3. Thus, we formulate
Numerical Fact 4.1 For Bernoulli lines defined by assumptions (a)-(e), δ << 1.
Based on this fact, the following can be proved:
Theorem 4.3 For Bernoulli lines defined by assumptions (a)-(e), the accuracy of the production rate estimate is characterized by
where δ is defined by (4.44) and O(δ) denotes a quantity of the same order of magnitude as δ.
Proof: See Section 20.1.
Numerical investigation: The accuracy of the performance measure estimates (4.36)-(4.40) has been investigated numerically using a C++ code, which simulates production lines defined by assumptions (a)-(e). (This code is included in the Simulation function of the PSE Toolbox – see Section 19.9.) Since similar numerical investigations are carried out throughout this textbook, we define them below by a standard procedure.
Simulation Procedure 4.1:
- (1) Select initial status of each machine up with probability pi and down with probability 1 − pi, i =1,..., M.
- (2) For each line under consideration, carry out 20 runs of the simulation code.
- (3) In each run, use the first 20,000 time slots as a warm-up period and the subsequent 400,000 time slots to statistically evaluate PR, WIPi, STi and BLi; this results in 95% confidence intervals of less than 0.001 for PR; 0.02 for WIPi; and 0.002 for STi and BLi.
The accuracy of the estimates has been evaluated by
The results of this numerical investigation for the set of production lines (4.43) are shown in Figures 4.19 - 4.21. From this information, the following conclusions can be derived:
- In general,
provides a relatively accurate estimate of
PR; in most cases the error is within 1% and the largest error is about 3%.
- The accuracy of
is typically lower. - The highest accuracy of all estimates is for the uniform machine efficiency pattern.
- The lowest accuracy is for the inverted bowl and "oscillating" patterns.
- Lines, which are inverse of each other, result in identical accuracy of all estimates.
In spite of this limitation in accuracy, formulas (4.36)-(4.40) provide a useful analytical tool for performance evaluation of Bernoulli lines, especially taking into account that the parameters of the machines are rarely known on the factory floor with accuracy better than 5% -10%.
Figure 4.19: Accuracy of
Figure 4.20: Accuracy of
Figure 4.21: Accuracy of
Figure 4.22: Accuracy of
4.3.1 Static laws of production systems
Many engineering systems can be characterized by their static laws. For instance, a static law of mechanical systems is:
Sum of all forces acting on a rigid body = 0. (4.49)
The statics of electric circuits are described by
Sum of currents flowing through a node = 0 (4.50)
and by
Sum of all voltage drops and rises in a closed loop = 0. (4.51)
Similarly, the static law of production systems should describe the steady state flow of parts through this system. As it has been shown above, in the case of serial lines with Bernoulli machines, this flow is characterized by the steady states of the recursive procedure (4.30), i.e., by
In the same manner as (4.49)-(4.51) characterize static behavior of mechanical and electrical systems, (4.52) characterizes static properties of the production systems under consideration. Analyzing (4.52), and taking into account the expressions for the performance measures (4.36)-(4.40), we derive system-theoretic properties of serial lines with Bernoulli machines.
4.3.2 Reversibility
Consider a serial line L, defined by assumptions (a)-(e) of Subsection 4.2.1, and its reverse Lr (see Figure 4.23).
Figure 4.23: M-machine production line and its reverse
Theorem 4.4 The performance measures of a Bernoulli serial line, L, and its reverse, Lr, are related as follows:
Proof: See Section 20.1.
The reversibility property of production lines has practical implications. Several of them are mentioned below.
- Some argue that buffers at the end of the line should be larger than those at its beginning, since more work has been put into parts as they progress further along the line and are closer to being complete. The reversibility property says that the same effect can be ensured by reversing this argument. Thus, the argument of "end of the line" or "beginning of the line" is not valid for buffer capacity assignment.
- If all machines and buffers are identical and one machine can be improved, which one should it be, so that the production rate of a serial line with an odd number of machines is maximized? Similar to the above, the reversibility property leads to a conclusion that it should be the machine in the middle of the line.
- If all machines are identical and the total buffering capacity N* must be allocated among M − 1 buffers, how should their capacity be selected so that the production rate is maximized? The reversibility property tells us that the allocation, whatever it may be, must be symmetric with respect to the middle machine (if the number of machines is odd) or with respect to the two middle machines (if the number of machines is even). So, the only question that remains is whether the buffers should be of equal capacity or not (assuming that N* is divisible by M). It is easy to show that buffers of equal capacity do not maximize PR. Therefore, due to the reversibility property and the argument of the first bullet above, optimal allocation must be of an inverted bowl shape, i.e., larger buffers are in the middle of the line. This result is quite intuitive, since the machines at the beginning and the end of the line experience less perturbations (because the first machine is not starved and the last is not blocked) and, therefore, need less protection than those in the middle of the line. In practical terms, however, the "bowl effect" is not very significant: the difference of PRs under the uniform and the optimal bowl allocation is typically within 1% (which is often below the accuracy with which parameters of the machines are known). To illustrate this fact, consider a Bernoulli 5machine line with pi =0.9. Its production rate for the uniform allocation of N* = 20 is 0.8666, while according to the best bowl allocation, it is 0.8667.
4.3.3 Monotonicity
Could PR of a serial line decrease if the machines or buffers are improved? Intuitively, it is clear that this should not happen. A formal statement to this effect follows from (4.52):
Theorem 4.5 In Bernoulli serial lines, defined by assumptions (a)-(e) of Subsection 4.2.1,
is
- strictly monotonically increasing in Ni, i =1,..., M − 1;
- strictly monotonically increasing in pi, i =1,..., M;
Proof: See Section 20.1.
:
