9.3 Discrete-Event Simulation Principles and Monte Carlo Modeling

Computer simulation is a numerical technique that models the operation of complex stochastic systems over time. When industrial systems exhibit complex interactions, finite buffer constraints, non-exponential distributions, or dynamic routing logic that exceed the analytical capabilities of classical queuing theory or Markov chains, simulation becomes the primary tool for analysis and decision-making. The NCEES FE Reference Handbook outlines core simulation concepts under its Industrial and Systems Engineering section. Candidates must understand the operational architecture of Discrete-Event Simulation (DES), pseudo-random number generation, random variate generation via the Inverse Transform Method, and rigorous statistical output analysis of terminating and steady-state systems.

---

1. Discrete-Event Simulation vs. Continuous Simulation

Simulation models are broadly classified based on how system state variables evolve over time:

Continuous Simulation:                           Discrete-Event Simulation (DES):
 State Variable S(t)                              State Variable N(t)
    ▲                                                ▲
    │      ╭───╮                                     │                 ┌───────
    │     ╭╯   ╰╮                                    │         ┌───────┘
    │    ╭╯     ╰╮                                   │         │
    │   ╭╯       ╰╮                                  │ ┌───────┘
    │───╯         ╰───────                           │─┘
    └──────────────────────► Time t                  └─────────────────────────► Time t
    State changes continuously                       State changes only at discrete event times
    (Differential Equations)                         (t1, t2, t3, ...) via Event List

Continuous Simulation

In continuous simulation, state variables change continuously and smoothly over continuous time. Continuous models are typically formulated as systems of differential equations and advanced using fixed time increments $\Delta t$. Examples include modeling chemical kinetics in a distillation column, flight dynamics of aerospace vehicles, or temperature dissipation in a heat-treating furnace.

Discrete-Event Simulation (DES)

In discrete-event simulation, system state variables remain completely constant over time intervals and change instantaneously at isolated, discrete points in time. These discrete points are triggered by the occurrence of an event. Between consecutive events, no state changes occur, and the simulation clock can bypass inactive time intervals entirely.

Time-Advance Mechanisms

• Next-Event Time Advance (Event-Driven): The standard mechanism for DES. The simulation clock jumps directly from the timestamp of the current event to the timestamp of the next earliest scheduled event in the Future Event List (FEL). Inactive periods are completely skipped, maximizing computational efficiency.
• Fixed-Increment Time Advance (Time-Driven): The clock advances in uniform slices $\Delta t$. At each tick $t + \Delta t$, the engine checks whether any events occurred during that slice. While useful for continuous or real-time gaming models, it is computationally inefficient for stochastic queuing networks with sparse events.

---

2. Core Architecture and Components of DES

A discrete-event simulation model is organized around six foundational entities and data structures:

1. Entities: Dynamic objects that enter the model, navigate through processes, consume resources, wait in queues, and exit. In manufacturing, entities represent raw parts, pallets, or work orders; in service systems, they represent customers, hospital patients, or incoming calls.
2. Attributes: Local data values or variables attached to individual entities. Attributes remain with the entity as it flows through the system. Examples include part arrival timestamp (used to compute cycle time $W$), required processing recipe, physical weight, customer priority, or rework status.
3. Resources: Static or semi-static capacity-constrained system assets that provide service to entities. Resources have a defined capacity (number of units) and discrete operational states: Idle, Busy, Blocked, or Under Maintenance/Failed. An entity requests a resource; if available, the resource is seized, service begins, and upon completion, the resource is released.
4. Queues: Storage buffers where entities wait when requested resources are unavailable. Queues are governed by capacity limits and queuing disciplines (e.g., FIFO, LIFO, Earliest Due Date, Highest Priority).
5. Events: Instantaneous occurrences that alter the system state variables. Standard events in queuing models include:
   • Arrival Event: Generates a new entity, assigns its attributes, schedules the next arrival event, and checks resource availability.
   • End-of-Service / Departure Event: Releases the resource, collects entity statistics, routes the finished entity, and checks the queue for waiting entities.
6. Simulation Clock and Future Event List (FEL):
   • Clock ($T$): A global variable storing the current simulated time.
   • Future Event List (FEL): A dynamic priority queue data structure containing all scheduled future event notices, sorted in ascending order of scheduled event time ($t_1 \le t_2 \le t_3 \le \dots$).

                                FUTURE EVENT LIST (FEL)
                             ┌───────────────────────────────┐
                             │ [Event 1: Departure at 10:15] │ <── Head of List (Earliest)
                             │ [Event 2: Arrival at 10:22]   │
                             │ [Event 3: Breakdown at 10:45] │
                             └───────────────────────────────┘
                                             │
                                             ▼ Advance Clock to 10:15
                             ┌───────────────────────────────┐
                             │ Execute Event Routine:        │
                             │ - Update State Variables      │
                             │ - Free Resource               │
                             │ - Collect Statistics          │
                             │ - Schedule Next Event into FEL│
                             └───────────────────────────────┘

---

3. Pseudo-Random Number Generation (PRNG)

Stochastic simulations rely on pseudo-random numbers uniformly distributed on the continuous interval $(0, 1)$, denoted $U \sim U(0, 1)$. These numbers must be:

• Uniform: Equal probability of falling into any sub-interval of $(0, 1)$.
• Independent: No serial correlation between successive numbers.
• Reproducible: Given an identical initial starting integer (seed $X_0$), the algorithm must regenerate the exact same sequence for debugging and controlled experiments.
• Long Period: The sequence of generated numbers should not repeat for billions of iterations.

The Linear Congruential Generator (LCG)

The classical algorithm tested on the FE exam is the Linear Congruential Generator (LCG), governed by the recursive modular arithmetic relationship:

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

Where:

• $X_0$: Initial starting value (the seed, $X_0 \ge 0$)
• $a$: Multiplier ($a > 0$)
• $c$: Increment ($c \ge 0$); if $c = 0$, it is called a multiplicative LCG; if $c > 0$, a mixed LCG
• $m$: Modulus ($m > X_0, a, c$)
• $\pmod m$: The remainder after integer division by $m$

To map the resulting integer $X_i$ to a uniform continuous random number on $[0, 1)$:

$U_i = \frac{X_i}{m}$

Full-Period Conditions (The Hull-Dobell Theorem)

The maximum possible cycle length (period) of an LCG is $m$. By the Hull-Dobell Theorem, a mixed LCG ($c > 0$) achieves a full period of length $m$ if and only if:

1. $c$ and $m$ are relatively prime (their greatest common divisor is $1$: $\gcd(c, m) = 1$).
2. $a - 1$ is divisible by all prime factors of $m$.
3. $a - 1$ is divisible by $4$ if $m$ is divisible by $4$.

---

4. Random Variate Generation: The Inverse Transform Method

Once uniform random numbers $U \sim U(0, 1)$ are generated, they must be converted into random variates following specific physical distributions (e.g., exponential, normal, triangular, Weibull).

Theoretical Foundation: Probability Integral Transform

Let $X$ be a continuous random variable with cumulative distribution function $F(x) = P(X \le x)$. The Inverse Transform Method is based on the mathematical theorem that if $U \sim U(0, 1)$, then the random variable:

$X = F^{-1}(U)$

follows the distribution defined by $F(x)$.

       Cumulative Probability F(x)
           ▲
       1.0 │                     ╭────────────────────── F(x)
           │                   ╭╯
       U ──┼─────────────────╭─╯
           │                 │
           │                 │
       0.0 └─────────────────┼────────────────────────► Variate X
                             X = F^(-1)(U)

Derivation for the Exponential Distribution

The Exponential distribution with rate parameter $\lambda$ (and mean $\beta = 1/\lambda$) has cumulative distribution function:

$F(x) = 1 - e^{-\lambda x}, \quad x \ge 0$

Setting $F(X) = U$ and solving for $X$:

$1 - e^{-\lambda X} = U$ $e^{-\lambda X} = 1 - U$ $-\lambda X = \ln(1 - U)$ $X = -\frac{1}{\lambda} \ln(1 - U)$

Because $U$ is uniformly distributed on $(0, 1)$, its complement $1 - U$ is also identically distributed as $U(0, 1)$. Therefore, the standard computational formula is:

$X = -\frac{1}{\lambda} \ln(U) = -\beta \ln(U)$

Where $\beta = 1/\lambda$ is the mean of the exponential distribution.

Derivation for Continuous Uniform Distribution $U(a, b)$

The CDF of a uniform variable on the interval $[a, b]$ is $F(x) = \frac{x - a}{b - a}$. Setting $F(X) = U$ yields:

$X = a + (b - a)U$

Discrete Variate Generation via Table Lookup

For a discrete random variable with probability mass function $p(x_i) = P(X = x_i)$, the cumulative distribution $F(x_k) = \sum_{i=1}^k p(x_i)$ partitions the interval $[0, 1)$ into sub-intervals. If a generated uniform number $U$ falls within $[F(x_{k-1}), F(x_k))$, the variate $X = x_k$ is returned.

---

5. Monte Carlo Simulation Methodology

While discrete-event simulation explicitly tracks state dynamics across time, Monte Carlo simulation uses repeated stochastic sampling to solve static or deterministic mathematical models where time advance is not the governing factor.

Key Applications in Industrial Engineering

1. Engineering Economics Risk Analysis: Evaluating Net Present Value (NPV) or Internal Rate of Return (IRR) where future cash flows, interest rates, and asset lifespans follow subjective probability distributions rather than fixed single-point estimates.
2. Mechanical Tolerance Stack-Up: Simulating dimensional tolerance accumulation across multiple mated parts in complex assemblies to determine expected assembly defect rates.
3. Reliability Network Analysis: Estimating system reliability for complex topologies (bridge networks, non-series/parallel systems) where analytical derivations are intractable.

---

6. Statistical Output Analysis and Experimental Design

Because simulation inputs are random variates, simulation outputs are also random variables. A single simulation run represents merely one sample point ($n = 1$) from an underlying output population. Rigorous statistical analysis is required to draw valid engineering conclusions.

Terminating vs. Non-Terminating (Steady-State) Systems

Initialization Bias and the Warm-Up Period

When simulating a non-terminating steady-state system, the simulation typically begins with the model empty and idle (zero entities in queue and all servers idle). This creates initialization bias, because early observations systematically underestimate steady-state queue lengths, WIP, and server utilization.

• Remedy: Warm-Up Period Truncation: Allow the simulation to run for a warm-up period $T_w$ until transient startup effects dissipate and steady-state conditions are reached. All data collected during $[0, T_w]$ is discarded (truncated), and statistics are accumulated strictly over $[T_w, T_{\text{final}}]$.
• Welch's Method: A standard graphical procedure that plots moving averages of output metrics across multiple replications to visually identify the time point $T_w$ where the metric stabilizes.

  WIP Level L(t)
      ▲
      │                        Steady-State Phase (Data Recorded)
      │               ╭─╮     ╭───╮   ╭─╮     ╭─╮   ╭───╮
      │              ╭╯ ╰─────╯   ╰─╮╭╯ ╰─────╯ ╰───╯   ╰─╮
      │             ╭╯              ╰╯
      │            ╭╯
      │  Transient │◄────────────── Data Retained ──────────────►│
      │  (Startup) │
      │─╭──────────╯
      └─┴────────────────────────────────────────────────────────► Time t
        0          Tw (Warm-up Truncation Point)
        │◄────────►│ Data Discarded

Independent Replications and Confidence Intervals

To obtain statistically valid estimators for a performance measure $\mu_Y$:

1. Execute $R$ independent simulation runs (replications), each using different random number seeds and starting from identical initial conditions.
2. Let $Y_r$ be the sample mean output from replication $r$ ($r = 1, 2, \dots, R$). The point estimator of the overall mean is: $\bar{Y} = \frac{1}{R} \sum_{r=1}^R Y_r$
3. The sample variance across the $R$ independent replications is: $S^2 = \frac{1}{R - 1} \sum_{r=1}^R (Y_r - \bar{Y})^2$
4. A $(1 - \alpha)100%$ confidence interval for the true mean $\mu_Y$ is computed using the Student's $t$-distribution with $\nu = R - 1$ degrees of freedom: $\bar{Y} \pm t_{\alpha/2, R-1} \frac{S}{\sqrt{R}}$

Determining Required Number of Replications ($R^*$)

If an engineer requires an estimate with a specified maximum allowable half-width error $\epsilon$ at confidence level $(1 - \alpha)$, the required total number of replications is approximated by:

$R^* \ge \left( \frac{t_{\alpha/2, R-1} \cdot S}{\epsilon} \right)^2 \approx \left( \frac{z_{\alpha/2} \cdot S}{\epsilon} \right)^2$

Variance Reduction Techniques (VRT)

Variance reduction techniques improve the statistical precision (narrowing the confidence interval) of simulation estimators without increasing sample size or simulation run length:

• Common Random Numbers (CRN): Used when comparing alternative system configurations (e.g., System A vs. System B). By synchronizing the exact same random number streams across both models (e.g., identical part arrival times and identical machine cycle times), the variance of the difference is reduced: $\text{Var}(Y_A - Y_B) = \text{Var}(Y_A) + \text{Var}(Y_B) - 2 \text{Cov}(Y_A, Y_B)$ Positive correlation ($\text{Cov}(Y_A, Y_B) > 0$) significantly reduces $\text{Var}(Y_A - Y_B)$, ensuring differences reflect design changes rather than stochastic noise.
• Antithetic Variates (AV): Used for estimating the performance of a single system. For every run driven by uniform numbers $U_1, U_2, \dots$, a complementary paired run is driven by $1 - U_1, 1 - U_2, \dots$, inducing negative covariance ($\text{Cov}(Y_1, Y_2) < 0$) and reducing variance of the pooled average.

---

7. Step-by-Step Worked Engineering Examples

Worked Example 9.3.1: LCG and Exponential Variate Generation

Problem:

1. A mixed Linear Congruential Generator is defined by: $X_{i+1} = (7 X_i + 3) \pmod{16}$ Using a starting seed $X_0 = 5$, compute the next three integer values ($X_1, X_2, X_3$) and their associated uniform random variates ($U_1, U_2, U_3$).
2. Use the Inverse Transform Method to generate an interarrival time following an Exponential distribution with mean arrival rate $\lambda = 4$ parts per hour, using the generated uniform number $U_1$.

Solution:

1. Compute LCG Sequence:
   • Step 1 ($i = 0 \to 1$): $X_1 = (7 X_0 + 3) \pmod{16} = (7 \times 5 + 3) \pmod{16} = (35 + 3) \pmod{16} = 38 \pmod{16}$ $38 / 16 = 2 \text{ with a remainder of } 6 \implies X_1 = 6$ $U_1 = \frac{X_1}{m} = \frac{6}{16} = 0.3750$
   • Step 2 ($i = 1 \to 2$): $X_2 = (7 X_1 + 3) \pmod{16} = (7 \times 6 + 3) \pmod{16} = (42 + 3) \pmod{16} = 45 \pmod{16}$ $45 / 16 = 2 \text{ with a remainder of } 13 \implies X_2 = 13$ $U_2 = \frac{X_2}{m} = \frac{13}{16} = 0.8125$
   • Step 3 ($i = 2 \to 3$): $X_3 = (7 X_2 + 3) \pmod{16} = (7 \times 13 + 3) \pmod{16} = (91 + 3) \pmod{16} = 94 \pmod{16}$ $94 / 16 = 5 \text{ with a remainder of } 14 \implies X_3 = 14$ $U_3 = \frac{X_3}{m} = \frac{14}{16} = 0.8750$
2. Generate Exponential Variate via Inverse Transform:
   • Arrival rate: $\lambda = 4$ parts/hour
   • Mean interarrival time: $\beta = 1/\lambda = 1/4 = 0.25$ hours $= 15$ minutes
   • Apply inverse transform formula with $U_1 = 0.3750$: $X = -\frac{1}{\lambda} \ln(U_1) = -0.25 \times \ln(0.3750)$ $\ln(0.3750) \approx -0.98083$ $X = -0.25 \times (-0.98083) = 0.2452 \text{ hours}$
   • In minutes: $0.2452 \times 60 = 14.71$ minutes.

Worked Example 9.3.2: Replications, Confidence Interval, and Sample Sizing

Problem: A distribution center simulates an automated tote-picking cell across $R = 5$ independent replications. The observed average hourly throughputs (totes/hour) across the 5 replications are:

$Y = [124, 130, 122, 136, 128]$

1. Calculate the point estimate $\bar{Y}$ and sample standard deviation $S$.
2. Construct a $95%$ confidence interval for the true mean hourly throughput (for $\nu = 4$ degrees of freedom, $t_{0.025, 4} = 2.776$).
3. If management requires the estimate to have a maximum half-width error of $\epsilon = 2.0$ totes/hour with $95%$ confidence ($z_{0.025} = 1.96$), how many total replications $R^*$ are required?

Solution:

1. Point Estimator and Sample Variance:
   • Sample mean: $\bar{Y} = \frac{124 + 130 + 122 + 136 + 128}{5} = \frac{640}{5} = 128.0 \text{ totes/hour}$
   • Deviations about the mean: $(124 - 128)^2 = (-4)^2 = 16$ $(130 - 128)^2 = (2)^2 = 4$ $(122 - 128)^2 = (-6)^2 = 36$ $(136 - 128)^2 = (8)^2 = 64$ $(128 - 128)^2 = (0)^2 = 0$ $\sum (Y_r - \bar{Y})^2 = 16 + 4 + 36 + 64 + 0 = 120$
   • Sample variance and standard deviation: $S^2 = \frac{120}{5 - 1} = \frac{120}{4} = 30.0 \implies S = \sqrt{30.0} \approx 5.477 \text{ totes/hour}$
2. Construct $95%$ Confidence Interval:
   • Standard error of the mean: $\text{SE} = \frac{S}{\sqrt{R}} = \frac{5.477}{\sqrt{5}} = \frac{5.477}{2.2361} \approx 2.449 \text{ totes/hour}$
   • Half-width error: $H = t_{0.025, 4} \times \text{SE} = 2.776 \times 2.449 \approx 6.80 \text{ totes/hour}$
   • $95%$ Confidence Interval: $128.0 \pm 6.80 \implies [121.20, \; 134.80] \text{ totes/hour}$
3. Determine Required Replications ($R^*$):
   • Target half-width $\epsilon = 2.0$ totes/hour: $R^* \ge \left( \frac{z_{0.025} \cdot S}{\epsilon} \right)^2 = \left( \frac{1.96 \times 5.477}{2.0} \right)^2 = \left( \frac{10.735}{2.0} \right)^2 = (5.368)^2 \approx 28.81$
   • Rounding up to the next integer yields $R^ = 29$ replications*.

---

8. NCEES Reference Handbook Tips & Realistic Exam Traps

• Rate vs. Mean in Inverse Transform: For an exponential variate, the formula is $X = -(1/\lambda)\ln(U) = -\beta \ln(U)$. If given a mean service time of $\beta = 10$ minutes, do not write $-1/10 \ln(U)$. The parameter $\lambda$ is the rate ($1/10$), so $1/\lambda = 10$. Always check whether the problem specifies the rate $\lambda$ or the mean $\beta$.
• The Negative Sign in Exponential Generation: Because $U \in (0, 1)$, $\ln(U)$ is always a negative number. The leading minus sign in $X = -(1/\lambda)\ln(U)$ is mandatory to ensure the generated physical time $X$ is positive.
• Common Random Numbers Scope Trap: A frequent FE exam conceptual distractor claims that Common Random Numbers (CRN) reduces the variance of an absolute performance metric in a single system. This is false! CRN is strictly used when comparing two or more alternative designs to induce positive covariance and reduce the variance of the difference $(\bar{Y}_A - \bar{Y}_B)$.
• Seed Management: To achieve statistically independent replications, each simulation run must use a distinct, uncorrelated seed. Running five replications with the exact same starting seed $X_0$ will produce five identical runs, yielding zero sample variance and invalidating statistical conclusions.