Monte Carlo Simulation

Monte Carlo simulation answers a question about an uncertain system by drawing random numbers and reading each one as an event: a day's demand, a lead time, a treatment, a service time. This section solves all five examples of Exercise 5 step by step: inventory control, a dentist clinic, the inverse transform for the exponential distribution, a queue with exponential times, and a manufacturing facility.

Objectives

• Turn a frequency table into two-digit random number intervals.
• Compute the mean and standard deviation of demand and lead time, then the economic order quantity (EOQ) and the reorder level (ROL).
• Build the inventory simulation table and decide when to order.
• Simulate a queue with scheduled arrivals and find the average waiting time.
• Generate exponential random variates with the inverse transform.
• Use a mixed congruential generator to simulate a manufacturing facility.

1. Example 1: inventory control

Suppose you have the demand for the past 100 days and the lead time in which the ordered quantity arrived. We need to know the demand in the coming days, the possible lead time, and when to order.

Demand for the past 100 days:

Lead time:

For the lead time we may have a standard frequency distribution rather than the given frequencies, to find the random numbers that correspond to a particular lead time.

Suppose that the inventory control has these costs:

When you have a shortage, you have to bring items from the local market at a higher price.

Step 1: mean and standard deviation

The reorder level and the other quantities come from statistical calculations on the available data. The 100 demand values add up to $\sum D_i = 4523$, and their squares add up to $\sum D_i^2 = 214077$.

The mean demand and its standard deviation are:

$$\begin{aligned} \mu_D &= \dfrac{\sum D_i}{n} = \dfrac{4523}{100} \\ &= 45.23 \text{ unit/day} \end{aligned}$$ $$\begin{aligned} \sigma_D &= \sqrt{\dfrac{\sum D_i^2 - n (\mu_D)^2}{n - 1}} \\ &= \sqrt{\dfrac{214077 - 100 (45.23)^2}{99}} \\ &= 9.8 \end{aligned}$$

The mean lead time and its standard deviation are:

$$\begin{aligned} \mu_L &= \dfrac{(28 \times 3) + (44 \times 4) + (28 \times 5)}{100} \\ &= \dfrac{400}{100} \\ &= 4 \text{ days} \end{aligned}$$ $$\begin{aligned} & 28 (3)^2 + 44 (4)^2 + 28 (5)^2 - 100 (4)^2 \\ &= 252 + 704 + 700 - 1600 \\ &= 56 \end{aligned}$$ $$\begin{aligned} \sigma_L &= \sqrt{\dfrac{56}{99}} \\ &= 0.75 \end{aligned}$$

Note on the exercise sheet. The sheet writes the mean lead time as 400/4 = 4. The division is by 100: 400 / 100 = 4. The answer of 4 days stands.

Step 2: random numbers for the demand

We find the histogram of the demand, calculate the frequency and the cumulative frequency, and assign two-digit random numbers (from 00 to 99) to each interval according to its frequency. An interval with frequency f takes f numbers, and it ends at its cumulative frequency minus 1. When we predict the demand, we take the average value of the demand in each interval.

Note on the exercise sheet. The sheet says the demand ranges from 20 to 70. The data runs from 22 to 72. The table is correct: its eleven intervals, from 20 - 24 to 70 - 74, hold all 100 values, and counting the values again gives the same frequencies.

Step 3: random numbers for the lead time

The same calculation applies to the lead time:

Step 4: the economic order quantity

$$\begin{aligned} EOQ &= \sqrt{\dfrac{2 \mu_D C_0}{C_i}} \\ &= \sqrt{\dfrac{2 \times 45.23 \times 100}{0.1}} \\ &\approx 300 \text{ units} \end{aligned}$$

Whenever the stock falls below a certain limit, the reorder level, we order this same quantity of 300 units, because this is a fixed quantity system.

Step 5: the reorder level

$$\begin{aligned} ROL &= \mu_D \mu_L \\ &\quad + K \sqrt{\mu_L \sigma_D^2 + \mu_D^2 \sigma_L^2} \end{aligned}$$

K is a constant that sets the probability of a stock shortage. For K = 1.65 the chance of a shortage is 5%, so we are 95% sure that there will not be any shortage. For K = 2.33 there is only a 1% chance of a shortage.

For K = 1.65:

$$\sqrt{4 (9.8)^2 + (45.23)^2 (0.75)^2} = 39.18$$ $$\begin{aligned} ROL &= (45.23 \times 4) + 1.65 \times 39.18 \\ &= 180.92 + 1.65 \times 39.18 \\ &\approx 245.6 \end{aligned}$$

Note on the exercise sheet. The sheet writes ROL = 245. The exact value is about 245.6, so 245 is rounded down. We keep 245, because the simulation table uses it. With 246 the orders fall on the same days, since the lowest quantity that gives no order is 248.

Step 6: starting the simulation

We start the simulation with EOQ = 300 and ROL = 245 units. We assume the initial stock is 300 units on Friday evening, so the simulation runs from Monday onwards.

The random numbers from 00 to 99 can be generated in Excel or by using other random number generators.

Note on the exercise sheet. The sheet says to put 100 into the rand() function. In Excel, RAND() takes no argument, and it gives a number below 1. Multiply by 100 instead: INT(RAND()*100) gives a number from 00 to 99.

The sheet takes its random numbers from a printed table of two-digit numbers:

• For demand, from the left of the first row: 81, 78, 89, 59, 40, 33, 60, 67, 19, 21.
• For lead time, from the right of the last row: 65, 67, 22. We draw one only when we place an order.

Note on the exercise sheet. The text names the demand numbers as 81, 78, 59, which skips one. The first row and the simulation table both have 81, 78, 89, 59. Wednesday uses 89, so its demand is 57.

Step 7: the simulation table

The table covers two working weeks (10 working days), with initial stock = 300 units, EOQ = 300 units and ROL = 245 units:

The rules behind each row:

1. The initial stock is the previous final stock.
2. The demand comes from the random number through the table in Step 2: 81 falls in 69 - 81, so the demand is 52.
3. Final stock = initial stock - daily demand + quantity received.
4. Quantity in hand + order = final stock + the order on the way.
5. If that quantity falls below ROL, we order EOQ and draw a lead time from the table in Step 3.

Now follow the first days:

• First Monday: 300 - 52 = 248 units. This is above 245, so there is no order.
• Tuesday: the initial stock is 248 and the final stock is 248 - 52 = 196. This is below 245, so we have to reorder. The random number 65 falls in 28 - 71, so the lead time is 4 days.
• Wednesday to Friday: the final stock is below 245, but the 300 units on the way count: 139 + 300 = 439, 92 + 300 = 392 and 50 + 300 = 350. There is no new order.
• Second Wednesday: the stock has fallen to 214 with nothing on the way, so we order again. The random number 67 gives a lead time of 4 days.

Note on the exercise sheet. The sheet defines quantity in hand + order as final stock + quantity ordered. On the day of an order that gives 196 + 300 = 496 and 214 + 300 = 514, but the table shows 196 and 214. In the table, an order joins this column from the next day. We keep the table, and the order days are the same either way.

Step 8: when the orders arrive, and the costs

When the lead time is 4 working days from Tuesday, the quantity ordered is received on next Monday, and the final stock is:

$$(50 - 42) + 300 = 308 \text{ units}$$

The second ordered quantity is received on next Tuesday, after the 10 days of the table.

We can also calculate parameters like the cost in a certain period:

• Cost of holding inventory = holding cost per unit × final stock.
• Cost of ordering = number of orders × order cost.

2. Example 2: a dentist clinic

At a dentist clinic, apply Monte Carlo simulation to find the average waiting time of customers. The simulation starts at 8 a.m., and the arrival of customers is scheduled (in other cases it may be random, depending on some probability distribution). The first customer arrives at 8 a.m., and one customer comes every 30 minutes for one of five different treatments. The five treatment services are filling, crown, cleaning, extraction and checkup, with different probabilities and service times.

Step 1: the random number intervals

We do not know the treatment service a customer is coming for. That guess comes from random numbers generated by any generator, through this table:

Step 2: which service?

The random numbers 40, 82, 11, 34, 25, 66, 17, 79 predict the service for 8 patients:

Note that 40 is a crown, not a filling: the filling interval ends at 39.

Step 3: the simulation table

A service starts at the later of the arrival and the previous service end. The waiting time is the service start minus the arrival.

For example, patient 2 arrives at 8:30 while the dentist is busy until 9:00, so the wait is 30 minutes. That service ends at 9:00 + 15 = 9:15, so patient 3, who arrives at 9:00, waits 15 minutes.

The average waiting time is:

$$\begin{aligned} & 0 + 30 + 15 + 25 + 35 \\ & \quad + 45 + 30 + 40 = 220 \end{aligned}$$ $$\dfrac{220}{8} = 27.5 \text{ minutes}$$

Note on the exercise sheet. The sheet asks for the average waiting time but does not state it. Its table also runs off the page: it prints 6 of the 8 waits and 7 of the 8 service ends. The last values (the end 12:55 and the waits 30 and 40) and the average of 27.5 minutes come from the same rule.

3. Example 3: the inverse transform for the exponential distribution

Let $X$ have the exponential distribution with mean $\beta$. The distribution function is:

$$F(x) = \begin{cases} 1 - e^{-x/\beta} & \text{if } x \ge 0 \\ 0 & \text{otherwise} \end{cases}$$

To find $F^{-1}$, we set $u = F(x)$ and solve for $x$:

$$\begin{aligned} u &= 1 - e^{-x/\beta} \\ e^{-x/\beta} &= 1 - u \\ F^{-1}(u) = x &= -\beta \ln (1 - u) \end{aligned}$$

To generate the desired random variate, we first generate $U \sim U(0, 1)$ and then let:

$$X = -\beta \ln U$$

It is possible in this case to use $U$ instead of $1 - U$, since $1 - U$ and $U$ have the same $U(0, 1)$ distribution. This saves a subtraction.

4. Example 4: a queue with exponential times

Suppose that in a queueing system the inter-arrival time is exponentially distributed with mean $\beta_1 = 2$, and the service time is also exponentially distributed with mean $\beta_2 = 1.33$. From the simulation table for the arrival, service and departure times of 8 customers, calculate the average waiting time for a customer.

Step 1: the random variates

The inter-arrival time $X_i$ and the service time $Y_i$ come from the inverse transform equations:

$$\begin{aligned} X_i &= -\beta_1 \ln U_1 \\ Y_i &= -\beta_2 \ln U_2 \end{aligned}$$

Here $U_1$ and $U_2$ lie between 0 and 1. Given a stream of two-digit random numbers $u_1, u_2$ from 00 to 99, we get $U_1, U_2$ by dividing each number by 100. The streams are:

• u1 = 20, 23, 86, 09, 92, 35, 38, 01
• u2 = 21, 44, 27, 70, 73, 36, 59, 42

For the first customer:

$$\begin{aligned} U_1 &= 20 / 100 = 0.20 \\ X_1 &= -2 \ln (0.20) = 3.22 \\ U_2 &= 21 / 100 = 0.21 \\ Y_1 &= -1.33 \ln (0.21) = 2.08 \end{aligned}$$

Step 2: the simulation table

Each variate is rounded to two decimals and then added, as the sheet does. The arrival time is the previous arrival plus $X_i$. A service starts at the later of the arrival and the previous departure, and the customer departs after $Y_i$. The waiting time is the service start minus the arrival.

Only two customers wait:

• Customer 3 arrives at 6.46, and customer 2 leaves at 7.25: the wait is 7.25 - 6.46 = 0.79, and the departure is 7.25 + 1.74 = 8.99.
• Customer 5 arrives at 11.45, and customer 4 leaves at 11.75: the wait is 11.75 - 11.45 = 0.30, and the departure is 11.75 + 0.42 = 12.17.

The average waiting time is:

$$\dfrac{0.79 + 0.30}{8} = \dfrac{1.09}{8} \approx 0.14$$

We can calculate other measures of performance from the same table.

Note on the exercise sheet. For u2 = 27, 70, 42 the sheet has service times 1.75, 0.48, 1.16. With the stated $\beta_2 = 1.33$ they are 1.74, 0.47, 1.15, and the table above uses these. The sheet's whole column matches $\beta_2 = 4/3$ instead. So the sheet's departures of 9.00, 11.76, 12.18, 25.86 become 8.99, 11.75, 12.17, 25.85, and its second wait of 0.31 becomes 0.30. The sheet's average, (0.79 + 0.31) / 8 = 0.14, rounds to the same value.

5. Example 5: manufacturing

Given the units produced in a manufacturing facility and the number of vehicles available for transportation at a particular station, with these probabilities:

Step 1: the mixed congruential generator

A mixed congruential generator with the three parameters a, c and m generates two-digit random numbers between 00 and 99:

$$X_i = (a X_{i-1} + c) \bmod m$$

• For $u_1$: a = 21, c = 3, m = 100, seed $X_0 = 20$.
• For $u_2$: a = 41, c = 7, m = 100, seed $X_0 = 1$.

Note on the exercise sheet. The sheet writes $X_i = a X_{i-1} + c \bmod (m)$, without the outer brackets. Its steps compute the bracketed form, for example 423 mod 100 = 23, and every value of both streams follows that form.

For $u_1$:

$$\begin{aligned} X_1 &= (21 \times 20 + 3) \bmod 100 \\ &= 423 \bmod 100 = 23 \\ X_2 &= (21 \times 23 + 3) \bmod 100 \\ &= 486 \bmod 100 = 86 \end{aligned}$$

For $u_2$:

$$\begin{aligned} X_1 &= (41 \times 1 + 7) \bmod 100 \\ &= 48 \bmod 100 = 48 \\ X_2 &= (41 \times 48 + 7) \bmod 100 \\ &= 1975 \bmod 100 = 75 \end{aligned}$$

And so on to $X_9$:

Step 2: the simulation table for 10 days

The ten days use $X_0$ to $X_9$, so day 1 uses the seed itself. Each number is read through the interval table above.

Average number of units produced per day:

$$\dfrac{6000}{10} = 600 \text{ units/day}$$

Average number of available vehicles per day:

$$\dfrac{64}{10} = 6.4 \text{ vehicle/day}$$

6. Common mistakes

1. Using the limits of an interval as the demand. The demand drawn from an interval is its average: 81 falls in 69 - 81, and that row gives 52.
2. Missing a boundary. 40 is a crown because the filling interval ends at 39, and 79 is an extraction because its interval is 70 - 79.
3. Testing the final stock alone against ROL. On the first Wednesday the final stock is 139, but the 300 units on the way give 439, so there is no new order.
4. Counting lead time in calendar days. The count is in working days: 4 working days from Tuesday is next Monday.
5. Starting every service on arrival. A service starts at the later of the arrival and the previous end. Waiting time is measured from the arrival.
6. Forgetting to divide by 100, or dropping the minus sign. $U = u / 100$ and $X = -\beta \ln U$. Since $\ln U$ is negative, the minus sign makes $X$ positive.

7. Practice

Solve each problem, then check the answers at the end of this page.

1. With the demand and lead time tables of Example 1, find the demand for the random numbers 94, 51, 37 and the lead time for 20, 71, 72.
2. Use the cost formulas of Example 1 to find the cost of holding inventory and the cost of ordering for the 10 days of the simulation table.
3. At the dentist clinic of Example 2, 4 patients arrive every 30 minutes from 8:00, with the random numbers 39, 55, 80, 70. Find the waiting times and their average.
4. With $\beta_1 = 2$ and $\beta_2 = 1.33$ as in Example 4, simulate 3 customers with u1 = 50, 90, 40 and u2 = 30, 60, 80. Find the average waiting time.
5. Use the $u_1$ generator of Example 5 with the seed $X_0 = 5$ to get 4 random numbers (day 1 uses the seed), then find the units produced on each day and their average.

Key takeaways

• An interval takes as many two-digit numbers as its frequency, and it ends at its cumulative frequency minus 1.
• EOQ sets how much to order, and ROL sets when.
• An order counts from the day after it is placed, and it enters the stock on the day it arrives.
• In a queue, the waiting time is the service start minus the arrival.
• The inverse transform gives the exponential variate $X = -\beta \ln U$.
• A mixed congruential generator computes $X_i = (a X_{i-1} + c) \bmod m$.

Answers

Practice 1

94 falls in 94 - 95, so the demand is 62. 51 falls in 26 - 51, so the demand is 42. 37 also falls in 26 - 51, so the demand is 42. For the lead time, 20 falls in 00 - 27 and gives 3 days, 71 falls in 28 - 71 and gives 4 days, and 72 falls in 72 - 99 and gives 5 days.

Practice 2

The final stocks add up to 248 + 196 + 139 + 92 + 50 + 308 + 261 + 214 + 177 + 140 = 1825 unit-days. The cost of holding inventory is 0.1 × 1825 = 182.5 pounds. The table places 2 orders, so the cost of ordering is 2 × 100 = 200 pounds. There is no shortage, so there is no shortage cost.

Practice 3

39 is a filling (40 min), 55 is a cleaning (15 min), 80 is a checkup (15 min) and 70 is an extraction (45 min).

The average waiting time is 10 / 4 = 2.5 minutes.

Practice 4

The average waiting time is (1.39 + 0.24) / 3 = 0.54.

Practice 5

$$\begin{aligned} X_1 &= (21 \times 5 + 3) \bmod 100 \\ &= 108 \bmod 100 = 8 \\ X_2 &= (21 \times 8 + 3) \bmod 100 \\ &= 171 \bmod 100 = 71 \\ X_3 &= (21 \times 71 + 3) \bmod 100 \\ &= 1494 \bmod 100 = 94 \end{aligned}$$

The numbers 05, 08, 71, 94 give 500, 550, 650 and 700 units. The total is 2400, so the average is 2400 / 4 = 600 units/day.