# 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)


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)


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)


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)


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)


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)


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)


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)


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)


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)


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)


CONWIP works!