1. Discrete Event Simulation

What a discrete-event simulation consists of — entities, events, the event list, the clock, and statistical counters — walked through on a single-server FCFS queue.

1.1. Stochastic Dynamic Systems

A discrete-event simulation model consists of:

1. entities — the objects in the system
2. attributes — properties of an entity
3. state — variables describing the system at time $t$
4. events — instants at which the state changes
5. event list — the scheduled future events, in time order
6. clock — the current simulated time
7. statistical counters — running totals used to compute the outputs

1.2. Single-server queueing system with FCFS service discipline

Here is an example:

flowchart LR
    A([arrival]) --> Q[waiting line] --> S((server)) --> D([departure])

Objective. Estimate the average delay of a customer, where

$$\text{delay} = \text{(time service begins)} - \text{(arrival time)}.$$

Let $D_i$ be the delay of the $i$-th customer. Then $\bar{D} = \dfrac{1}{n}\sum_{i=1}^n D_i$.

Entities: customer, server
Attributes: arrival time of a customer; service time of a customer
State: $N(t)$ = number of customers in the system at time $t$; $S(t) \in \{\text{busy}, \text{idle}\}$ = status of the server
Events: arrival, departure, end of simulation
Statistical counters: number of customers delayed, total delay $\sum D_i$, area under $N(t)$, area under $Q(t)$

1.2.1 Notation

• $X_i$: interarrival time; $a_i = a_{i-1} + X_i$ is the arrival time ($a_0 = 0$), so the 2nd customer arrives at $X_1 + X_2$.
• $Y_i$: service time.
• $d_i$: departure time (the figure writes $D_i$; it should be $d_i$).
• Delay $D_i$: time spent in the queue, excluding service.
• Waiting time $W_i$: time spent in the system, including service, so $W_i = D_i + Y_i$.

Under FCFS, customer $i$ begins service when customer $i-1$ departs (or immediately, if the server is idle):

$$D_i = \max(0,\; d_{i-1} - a_i), \qquad d_i = a_i + D_i + Y_i .$$

1.2.2. Sample path (from the board, $T = 180$)

$i$	1	2	3	4	5
$X_i$	55	32	24	40	12
$a_i$	55	87	111	151	163
$Y_i$	43	36	34	16	–
$d_i$	98	134	168	184	–
$D_i$	0	11	23	17	21

gantt
    title Single-server queue, T = 180
    dateFormat X
    axisFormat %s
    tickInterval 20second
    section Customer 1
    in service   :done, c1s, 55, 98
    section Customer 2
    in queue     :active, c2q, 87, 98
    in service   :done,   c2s, 98, 134
    section Customer 3
    in queue     :active, c3q, 111, 134
    in service   :done,   c3s, 134, 168
    section Customer 4
    in queue     :active, c4q, 151, 168
    in service   :done,   c4s, 168, 184
    section Customer 5
    in queue     :active, c5q, 163, 184

Equivalently, $N(t)$ as a step function:

interval	$[0,55)$	$[55,87)$	$[87,98)$	$[98,111)$	$[111,134)$	$[134,151)$	$[151,163)$	$[163,168)$	$[168,180]$
$N(t)$	0	1	2	1	2	1	2	3	2
$Q(t)$	0	0	1	0	1	0	1	2	1

where $Q(t) = \max(0, N(t) - 1)$ is the number in the queue.

1.2.3. Time-average quantities (M/M/1 queue)

Average number in system.

$$\hat{L} = \sum_{i} i \cdot \frac{\text{total time with } i \text{ customers in system}}{T} = \frac{1}{T}\int_0^T N(t)\,dt = \frac{193}{180}.$$

Average number in queue.

$$\hat{L}_Q = \frac{1}{T}\int_0^T Q(t)\,dt = \frac{11 + 23 + 12 + 10 + 12}{180} = \frac{68}{180}.$$

Server utilization (proportion of time the server is busy).

$$\hat{\rho} = \frac{1}{T}\int_0^T \mathbf{1}\{S(t) = \text{busy}\}\,dt = \frac{180 - 55}{180}.$$

If $X_i \sim \text{Exp}(\lambda)$ and $Y_i \sim \text{Exp}(\mu)$, this is the $M/M/1$ queue (Markovian, i.e. exponential, interarrival times; exponential service times; one server).

1.2.4. How the simulation is driven

flowchart LR
    U["random numbers<br/>U ~ U(0,1)"] --> Tr[transformation]
    In[input distributions] --> Tr
    Tr --> V["random variates<br/>X_i, Y_i"] --> P[simulation program] --> O["outputs<br/>D_1, D_2, ..., W_1, W_2, ..."]

1.3. Simulation of an Inventory System (single commodity)

Problem

flowchart RL
    C([customer]) -->|demand| W["warehouse: inventory<br/>(managed by the inventory manager)"]
    W -->|order| SUP([supplier])
    SUP -.->|delivery, after the lead time| W

• When to order?
• How much to order?

(Textbook p. 48, §1.5.)

The times between successive demands are i.i.d. $\text{Exp}(\lambda)$ with $1/\lambda = 0.1$ month (on average $\lambda = 10$ demands per month).

The sizes of demands are i.i.d. (and independent of the demand times) with

$$D = \begin{cases} 1, & \text{w.p. } 1/6, \\ 2, & \text{w.p. } 1/3, \\ 3, & \text{w.p. } 1/3, \\ 4, & \text{w.p. } 1/6. \end{cases}$$

At the beginning of each month, the inventory level is reviewed and a decision is made whether to place an order and how much to order. If an order is placed, the order lead time (delivery time) is uniformly distributed between 0.5 and 1 month, $U(0.5, 1)$.

Let $I$ = the inventory level at the beginning of a month. The following $(s, S)$ policy is employed. Let $Z$ be the order amount:

$$Z = \begin{cases} S - I, & \text{if } I < s, \\ 0, & \text{if } I \ge s. \end{cases}$$

So $s$ answers when to order (order as soon as the level at a review is below $s$), and $S$ answers how much (order up to $S$).

$I$ could be negative (demand backlog): a demand larger than the stock on hand is not lost; the excess is backlogged and filled later. When an order arrives, it first clears the backlog, and the rest goes into inventory.

One can show that for inventory models of this kind (a fixed cost per order, plus per-item ordering, holding and backlog costs) the optimal policy takes this $(s, S)$ form (Scarf, 1960). It is a theorem for this class of models, not for every inventory system. Here the $(s, S)$ form is taken as given, and simulation is used to compare different values of $s$ and $S$.

Cost structure

If an order is placed, the cost is $k + c\,Z$ (no order, no cost).

$k$ is called the setup cost. $c$ is the per-unit ordering cost: the price paid for each item ordered, so this part grows with the order size, whereas $k$ is paid once per order however many items are ordered (paperwork, transportation). With the textbook’s values $k = \$32$, $c = \$3$, ordering $Z = 20$ items costs $32 + 3 \cdot 20 = \$92$. Because of $k$, many small orders are expensive, which pushes toward ordering rarely and in large amounts.

There is a holding cost of $\$h$ per item per month held in inventory (penalizes high inventory levels). Textbook: $h = \$1$.

There is a penalty cost of $\$p$ per item in backlog per month (penalizes negative inventory levels). Textbook: $p = \$5$ (written $\pi$ there).

Simulation

States: inventory level at time $t$, $I(t)$; and the amount $Z$ of the outstanding order, if any. (The lead time is at most 1 month and reviews are monthly, so at most one order is outstanding at a time.)

Events: demand arrival, order arrival, inventory evaluation (at the beginning of each month), end of simulation.

Statistical counters: total ordering cost, area under $I^+(t)$, area under $I^-(t)$ (defined below).

Let $I^+(t) = \max\{I(t), 0\}$: items physically held in inventory at time $t$.
$I^-(t) = \max\{-I(t), 0\}$: items in backlog at time $t$ (demand already received but not yet filled).

Why separate? $I(t) = I^+(t) - I^-(t)$, and at any time at most one of the two is non-zero. They are charged at different rates: $h$ per item per month on hand, $p$ per item per month in backlog. Averaging $I(t)$ itself would let the positive and negative parts cancel (10 items in stock for half the time and 10 items short for the other half averages to $0$, which would look free), and two different prices could not be applied.

Average number of items held in inventory per month:

$$\bar I^{+} = \frac{1}{n}\int_0^n I^+(t)\,dt,$$

average number of backlogged items per month:

$$\bar I^{-} = \frac{1}{n}\int_0^n I^-(t)\,dt,$$

where $n$ is the number of months simulated. These are time averages, computed from areas just like the area under $N(t)$ in §1.2.3.

$$\text{Average cost per month} = \text{average ordering cost per month} + h\,\bar I^{+} + p\,\bar I^{-},$$

where the average ordering cost per month is $\frac{1}{n}\sum_{\text{orders}}(k + cZ)$.

The simulation is run for several candidate pairs $(s, S)$, and the pair with the smallest average cost per month is preferred. (Textbook: $n = 120$ months, starting from $I(0) = 60$ with no outstanding order; nine $(s, S)$ pairs are compared.)

Generating the inputs (derivations and a worked example: §4.7):

• interdemand time: $-0.1\ln u$ (exponential with mean $0.1$);
• demand size: discrete inverse transform with $F(1) = 1/6$, $F(2) = 1/2$, $F(3) = 5/6$, $F(4) = 1$: return $1$ if $u \le 1/6$, else $2$ if $u \le 1/2$, else $3$ if $u \le 5/6$, else $4$;
• lead time: $0.5 + 0.5\,u$, i.e. $U(0.5, 1)$.