Simulating CONWIP

CONWIP is an interesting method to shift work from the shop floor to the orderbook. There are, however, numerious issues to sort out before CONWIP can be succesfully applied:

  • Which stations to include in the conwip loop?
  • What should the conwip level be?
  • What will be the throughput times on the shop floor?

Together with a few students I wrote the simulator below to obtain insight into these questions.

The simulator

First some standard imports.

import numpy as np
import scipy
from scipy.stats import poisson
import  matplotlib.pylab as plt

scipy.random.seed(30)

We consider the periodic review simulation of a queueing network consisting of three stations. External orders arrive at an orderbook which is represented as queue \(Q_0\). Orders are released from the orderbook, and sent to the station \(Q_1\). Once served at the first station the orders move on to station 2 \(Q_2\). After being processed at this last station, jobs depart from the network.

Let

  • \(a_i(t)\) be the number of arrivals at station \(i\) at the start of period \(t\)
  • \(d_i(t)\) be the number jobs that depart from station \(i\) at the end of period \(t\)
  • \(Q_i(t)\) be the amount of jobs at station \(i\) at the end of period \(t\).
  • \(C_i(t)\) be the amount of jobs that station \(i\) can serve during period \(t\).

It is clear that the following recursions must hold for the order book

\[\begin{split}d_0(t) &= \min\{Q_0(t) + a_0(t), C_0(t)\}, \\ Q_0(t) &= Q_0(t) + a_0(t) - d_0(t), \\ a_{1}(t) &= d_0(t). \\\end{split}\]

For station 1, the structure is slightly different:

\[\begin{split}d_1(t) &= \min\{Q_1(t) + a_1(t), C_1(t)\}, \\ Q_1(t) &= Q_1(t) + a_1(t) - d_1(t), \\ a_{2}(t+1) &= d_1(t). \\\end{split}\]

The last equation implements the behavior that the jobs that have been served by station \(1\) become available for transportation at the end of period \(t\). Hence, these jobs arrive at the start of the next period, i.e., period \(t+1\), at station 2. For station 2 similar recursions hold.

Under conwip the number of jobs that can be released at the start of a period from the orderbook must be such that the WIP on the shop floor does not exceed a certain maximal WIP level \(W\). In this case, the WIP at the end of period \(t\) consists of all jobs queued at station \(1\) and \(2\) plus the jobs in transit, i.e.,

\[W(t) = Q_1(t) + Q_2(t) + a_2(t+1)\]

Therefore, under conwip, the amount of jobs that can leave from the order book, i.e., the capacity, must be:

\[C_0(t+1) = W - W(t).\]

Note that \(W(t)\) is the WIP at the end of period \(t\), hence, \(W-W(t)\) orders can be released at the next period.

In code the above becomes:

def simulator(a0, C1, C2, W):
    WIP = np.zeros_like(a0)
    Q0 = np.zeros_like(a0)
    C0 = np.zeros_like(a0)
    Q1 = np.zeros_like(a0)
    Q2 = np.zeros_like(a0)
    d0 = np.zeros_like(a0)
    d1 = np.zeros_like(a0)
    d2 = np.zeros_like(a0)
    a1 = np.zeros_like(a0)
    a2 = np.zeros_like(a0)

    for t in range(1, len(a0)-1):
        d0[t] = min(Q0[t-1] + a0[t],  C0[t])
        Q0[t] = Q0[t-1] + a0[t] - d0[t]

        a1[t] = d0[t]
        d1[t] = min(Q1[t-1] + a1[t], C1[t])
        Q1[t] = Q1[t-1] + a1[t] - d1[t]

        a2[t+1] = d1[t]
        d2[t] = min(Q2[t-1] + a2[t], C2[t])
        Q2[t] = Q2[t-1] + a2[t] - d2[t]

        WIP[t] = Q1[t] + Q2[t] + a2[t+1]

        assert W >= WIP[t], "More jobs in the network than CONWIP
allows"
        C0[t+1] = W - WIP[t]

    return d2, Q0, Q1, Q2

Simulation with deterministic numbers

Lets set the conwip level \(W\) to 1 and see what happens when we choose some basic numbers.

W = 2

a0= np.ones(15)
C1 = np.ones_like(a0)
C2 = np.ones_like(a0)

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

print("The output of queue 2:")
print(d2)
The output of queue 2:
[ 0.  0.  0.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  0.]

If \(W=1\) we expect to see throughput of 0.5:

W = 1

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

print("The output of queue 2:")
print(d2)
The output of queue 2:
[ 0.  0.  0.  1.  0.  1.  0.  1.  0.  1.  0.  1.  0.  1.  0.]

And indeed, after some initialization effects, we see that the output alternates between 1 and 0.

A simulation with random data

Now try random data. Suppose that, on average, 9 jobs per period arrive at the system, but that both workstations only have capacity of 8 jobs per period.

n = 200
a0 = poisson.rvs(9, size = n)
C1 = poisson.rvs(8, size = n)
C2 = poisson.rvs(8, size = n)
W = 50

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

plt.plot(Q0, label="Q0")
plt.plot(Q1, label="Q1")
plt.plot(Q2, label="Q2")
plt.plot(Q1+Q2, label="Q1+Q2")
plt.legend(loc=2)
../../_images/conwip_figure5_1.png

It is clear that the jobs accumulate in the orderbook, i.e., in \(Q_0\), as it should. also the WIP remains limited.

The first station is a bottleneck

Suppose again that, on average, 9 jobs per period arrive at the system, but now workstation 1 has a capacity of 8 jobs per week and workstation 2 a capacity of 9.

C2 = poisson.rvs(9, size = n)

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

plt.plot(Q0, label="Q0")
plt.plot(Q1, label="Q1")
plt.plot(Q2, label="Q2")
plt.plot(Q1+Q2, label="Q1+Q2")
plt.legend(loc=2)
../../_images/conwip_figure6_1.png

Of course, most orders are queued in the orderbook, but on the floor the orders mostly queue in front of station 1.

The second station is a bottleneck

Now workstation 1 has a capacity of 9 jobs per period and workstation 2 a capacity of 8.

C1 = poisson.rvs(9, size = n)
C2 = poisson.rvs(8, size = n)

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

plt.plot(Q0, label="Q0")
plt.plot(Q1, label="Q1")
plt.plot(Q2, label="Q2")
plt.plot(Q1+Q2, label="Q1+Q2")
plt.legend(loc=2)
../../_images/conwip_figure7_1.png

And now the orders queue in front of station 2.

The second station is a bottleneck (2)

Up to now I considered a system in overload, i.e., that more jobs arrive than can be served. Like this, it is simple to see the benefits of CONWIP. What happens in the more natural situation, in which the stations have slightly higher capacity than the rate at which jobs arrive.

n = 200
a0 = poisson.rvs(8.5, size = n)
C1 = poisson.rvs(10, size = n)
C2 = poisson.rvs(9, size = n)

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

plt.plot(Q0, label="Q0")
plt.plot(Q1, label="Q1")
plt.plot(Q2, label="Q2")
plt.plot(Q1+Q2, label="Q1+Q2")
plt.legend(loc=2)
../../_images/conwip_figure8_1.png

This is interesting: in periods of high load, the queue at the orderbook builds up, and in such periods nearly all jobs on the floor build up at the queue in front of station 2. Thus, we see the behavior we like to see.

Investing in non-bottleneck capacity

One of the interesting claims in Chapter 10 of Factory Physics is that it may be advantageous to invest in increasing the capacity of non-bottleneck stations. Of course, if extending non-bottleneck capacity is just as expensive as bottleneck capacity, it is best to increase the capacity of the bottleneck. If, however, non-bottleneck capacity is much cheaper, it may pay off to invest in non-bottlenecks. Even though some people may have doubts about this, the simulations below show this point clearly.

In the first scenario, the capacity of both stations is nearly the same, the second is the bottleneck, and both perform close to the arrival rate.

W = 30

n = 20000
a0 = poisson.rvs(8.8, size = n)
C1 = poisson.rvs(10, size = n)
C2 = poisson.rvs(9, size = n)

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

plt.plot(Q0, label="Q0")
plt.plot(Q1, label="Q1")
plt.plot(Q2, label="Q2")
plt.plot(Q1+Q2, label="Q1+Q2")
plt.legend(loc=2)
../../_images/conwip_figure9_1.png

And the throughput is

print(d2.mean())
8.7528

Now we increase the production rate of, the non-bottleneck, station 1 from 10 to 20.

C1 = poisson.rvs(20, size = n)

d2, Q0, Q1, Q2 = simulator(a0, C1, C2, W)

plt.plot(Q0, label="Q0")
plt.plot(Q1, label="Q1")
plt.plot(Q2, label="Q2")
plt.plot(Q1+Q2, label="Q1+Q2")
plt.legend(loc=2)
../../_images/conwip_figure11_1.png

And the throughput is

print(d2.mean())
8.7787

The difference in queueing behavior is quite dramatic. In the first case, the queue in the order book becomes really large, while in the second case it is much smaller.

Interestingly, the throughputs are nearly the same in both cases. This is according to expectation, since the conwip level \(W\) is quite large, i.e., \(W=30\). Further experimentation shows that when \(W=20\) the network is unstable, even when the capacity of the first station is quite large, i.e., \(C = 20\).

A network of four stations

Suppose now that we have 5 stations in tandem and that station 2 is the bottleneck. I claim that stations 3, 4 and 5 need not be part of the conwip loop. Lets see whether I am right.

First I include all stations in the conwip loop.

The code is a bit different than before. I implement the queues as a matrix, such that the time runs along the columns and the queue lengths along the rows. For instance, Q[1,3] is the queue length of the first station at period 3.

W = 100

N=5    # number of workstations
L = 1000 # number of periods
C = poisson.rvs(10, size = (N+1,L))
C[2, :] = poisson.rvs(9, size = L) # the second workstation is the
bottleneck


Q = np.zeros((N+1,L), dtype = int)
d = np.zeros_like(Q)
a = np.zeros_like(Q)
WIP = np.zeros_like(Q[0,:])
a[0,:] = poisson.rvs(9.2, size = L)


for t in range(1, L-1): # run along the periods
    d[0,t] = min( Q[0,t-1] + a[0,t],  C[0,t] )
    Q[0,t] = Q[0, t-1] + a[0, t] -  d[0,t]
    a[1,t]= d[0,t]
    for n in range(1,N):
        d[n,t] = min(Q[n,t-1] + a[n,t], C[n,t])
        Q[n,t] = Q[n,t-1] + a[n,t] - d[n,t]
        a[n+1,t+1] = d[n-1,t]

    WIP[t] = np.sum(Q[1:,t]) + np.sum(a[2:,t+1])
    C[0,t+1] = W - WIP[t]


for n in range(0,N+1):
    plt.plot(Q[n,:], label="Q"+str(n))
plt.legend(loc=2)
../../_images/conwip_figure13_1.png

As expected, the orders remain mostly at the orderbook, as the arrival rate is higher than the service rate of station 2, and on the floor the jobs are queued in front of station 2.

A CONWIP loop that includes the stations up to and including the bottleneck

Next, include only stations 1 and 2 in the loop. We only have to change the line in which we compute the WIP.

for t in range(1, L-1): # run along the periods
    d[0,t] = min( Q[0,t-1] + a[0,t],  C[0,t] )
    Q[0,t] = Q[0, t-1] + a[0, t] -  d[0,t]
    a[1,t]= d[0,t]
    for n in range(1,N):
        d[n,t] = min(Q[n,t-1] + a[n,t], C[n,t])
        Q[n,t] = Q[n,t-1] + a[n,t] - d[n,t]
        a[n+1,t+1] = d[n-1,t]

    WIP[t] = np.sum(Q[[1,2],t]) + np.sum(a[2,t+1])
    C[0,t+1] = W - WIP[t]


for n in range(0,N+1):
    plt.plot(Q[n,:], label="Q"+str(n))
plt.legend(loc=2)
../../_images/conwip_figure14_1.png

Indeed, there is a stable queue in front of station 2, and there is hardly any queue at stations 3 and further down the line.

A loop without the bottleneck

Suppose we make the stupid choice to only include stations 3, 4 and 5 in the conwip loop, but not stations 1 and 2. Thus, the conwip loop does not contain the bottleneck station 2.

for t in range(1, L-1): # run along the periods
    d[0,t] = min( Q[0,t-1] + a[0,t],  C[0,t] )
    Q[0,t] = Q[0, t-1] + a[0, t] -  d[0,t]
    a[1,t]= d[0,t]
    for n in range(1,N):
        d[n,t] = min(Q[n,t-1] + a[n,t], C[n,t])
        Q[n,t] = Q[n,t-1] + a[n,t] - d[n,t]
        a[n+1,t+1] = d[n-1,t]

    WIP[t] = np.sum(Q[3:,t]) + np.sum(a[4:,t+1])
    C[0,t+1] = W - WIP[t]


for n in range(0,N+1):
    plt.plot(Q[n,:], label="Q"+str(n))
plt.legend(loc=2)
../../_images/conwip_figure15_1.png

This is also according to expectation. Jobs are released to the floor when the loop containing stations 3, 4 and 5 allow. Apparently, nearly all jobs wait in front of the station 2, i.e., on the shop floor rather than in the order book. It is clear that this is not a desirable situation.

A loop containing only the bottleneck, but not any other station

Finally, suppose we choose to only include station 2 in the loop, but not other stations. What will happen then?

for t in range(1, L-1): # run along the periods
    d[0,t] = min( Q[0,t-1] + a[0,t],  C[0,t] )
    Q[0,t] = Q[0, t-1] + a[0, t] -  d[0,t]
    a[1,t]= d[0,t]
    for n in range(1,N):
        d[n,t] = min(Q[n,t-1] + a[n,t], C[n,t])
        Q[n,t] = Q[n,t-1] + a[n,t] - d[n,t]
        a[n+1,t+1] = d[n-1,t]

    WIP[t] = np.sum(Q[2,t])  + np.sum(a[2,t+1])
    #C[0,t+1] = max(W - WIP[t], 0)
    C[0,t+1] = W - WIP[t]


for n in range(0,N+1):
    plt.plot(Q[n,:], label="Q"+str(n))
plt.legend(loc=2)
../../_images/conwip_figure16_1.png

Finally, a long simulation to see whether we missed anything. We have to set the arrival rate to 9, for otherwise the queueing process is unstable.

W = 100

N=5    # number of workstations
L = int(1e5) # number of periods
C = poisson.rvs(10, size = (N+1,L))
C[2, :] = poisson.rvs(9, size = L) # the second workstation is the
bottleneck


Q = np.zeros((N+1,L), dtype = int)
d = np.zeros_like(Q)
a = np.zeros_like(Q)
WIP = np.zeros_like(Q[0,:])
a[0,:] = poisson.rvs(9, size = L)


for t in range(1, L-1): # run along the periods
    d[0,t] = min( Q[0,t-1] + a[0,t],  C[0,t] )
    Q[0,t] = Q[0, t-1] + a[0, t] -  d[0,t]
    a[1,t]= d[0,t]
    for n in range(1,N):
        d[n,t] = min(Q[n,t-1] + a[n,t], C[n,t])
        Q[n,t] = Q[n,t-1] + a[n,t] - d[n,t]
        a[n+1,t+1] = d[n-1,t]

    WIP[t] = np.sum(Q[2,t]) + np.sum(a[2,t+1])
    C[0,t+1] = W - WIP[t]


for n in range(0,N+1):
    plt.plot(Q[n,:], label="Q"+str(n))
plt.legend(loc=2)
../../_images/conwip_figure17_1.png

CONWIP works!