2.3 Discrete-Event Simulation, Random Variates & Model Validation

2.3 Discrete-Event Simulation, Random Variates & Model Validation

Discrete-Event Simulation (DES) is an indispensable computational modeling paradigm used by industrial engineers to analyze complex, dynamic, stochastic systems that defy closed-form analytical solutions. While queueing theory models provide exact steady-state equations for idealized Markovian structures, real-world manufacturing plants, supply chains, and hospital operating suites feature non-exponential distributions, finite buffer blocking, sequence-dependent machine setups, routing recirculations, and dynamic labor scheduling that necessitate simulation.

---

Discrete-Event Simulation (DES) Architecture & System State Mechanics

Unlike continuous simulation (where state variables change continuously according to differential equations) or Monte Carlo simulation (which evaluates static stochastic relationships without a time dimension), Discrete-Event Simulation models systems where state changes occur instantaneously at discrete, separated points in time.

Core DES Elements

1. Entities: Dynamic objects that flow through the simulation model (e.g., parts, pallets, customers, work orders). Entities are created, travel through processes, seize resources, wait in queues, and are ultimately disposed.
2. Attributes: Local data values attached to individual entities that travel with them (e.g., entity creation timestamp, part geometry type, priority level, rework count).
3. Resources: Static operational units that provide service to entities (e.g., CNC milling machines, automated storage cranes, forklifts, certified technicians). A resource possesses a finite capacity and distinct states: IDLE, BUSY, FAILED, or BLOCKED.
4. Queues: Ordered lists of entities waiting to seize a resource when its capacity is fully utilized. Queues operate under specific ranking rules (FIFO, LIFO, Lowest Processing Time First, Earliest Due Date).
5. Events: Instantaneous occurrences that change the state of the system at a specific simulation clock instant. Common events include:
   • ARRIVAL_EVENT: A new entity enters the system.
   • DEPARTURE_EVENT / END_SERVICE_EVENT: A resource completes processing an entity.
   • BREAKDOWN_EVENT: An operational resource experiences an unscheduled failure.
6. Simulation Clock: A global variable that holds the current point in simulated time. It does not progress uniformly; rather, it leaps forward instantaneously to the timestamp of the next scheduled event.
7. State Variables: The set of mathematical variables required to completely describe the status of the system at clock time $t$ (e.g., $N_q(t)$: number of entities in queue at time $t$; $B(t)$: number of busy servers at time $t$).
8. Future Event List (FEL) / Event Calendar: A priority queue containing all scheduled future events, sorted chronologically by scheduled occurrence time ($t_{e1} \le t_{e2} \le \dots \le t_{ek}$).

The Next-Event Time-Advance Mechanism

The engine of any discrete-event simulator operates on the Next-Event Time-Advance algorithm:

1. Initialization: Set simulation clock $CLK = 0$. Initialize state variables (e.g., $N_q = 0, B = 0$). Schedule the initial arrival event in the FEL.
2. Identify Imminent Event: Retrieve the event from the head of the FEL with the smallest timestamp $t^* = \min \{t_e \in \text{FEL}\}$.
3. Advance Clock: Jump the simulation clock directly to the event time: $CLK = t^*$.
4. Execute Event Routine: Update system state variables corresponding to the event type. For example, if an END_SERVICE_EVENT occurs:
   • Decrement the count of busy resources.
   • If $N_q > 0$, remove the first entity from the queue, decrement $N_q$, assign it to the newly freed resource, and generate a new END_SERVICE_EVENT at $CLK + T_{service}$.
   • If $N_q = 0$, set resource state to IDLE.
5. Statistical Accumulation: Update time-weighted area accumulators for performance metrics (e.g., $\int N_q(t) dt$).
6. Check Stopping Criteria: If the target simulation run length or processed entity count is reached, terminate; otherwise, return to Step 2.

Loading diagram...

---

Pseudo-Random Numbers and Random Variate Generation

Stochastic simulation models require a stream of continuous uniform random numbers $U \sim U(0, 1)$ that are subsequently transformed into random variates representing operational phenomena (interarrival times, process durations, machine times-to-failure).

Linear Congruential Generators (LCG)

Pseudo-random numbers are generated deterministically using mathematical recurrence relations. The most widely implemented algorithm is the Linear Congruential Generator (LCG):

$Z_i = (a Z_{i-1} + c) \pmod m$ $U_i = \frac{Z_i}{m}$

Where:

• $m$: Modulus (a large prime number or power of 2, typically $2^{31}-1$ or $2^{64}$).
• $a$: Multiplier integer.
• $c$: Increment integer.
• $Z_0$: The initial seed value.

Under the Hull-Dobell Theorem, an LCG achieves its maximum possible cycle period of $m$ if and only if $c$ and $m$ are relatively prime, $(a - 1)$ is divisible by all prime factors of $m$, and $(a - 1)$ is divisible by 4 if $m$ is divisible by 4.

The Inverse Transform Method

The Inverse Transform Method is the primary mathematical technique for converting standard uniform random numbers $U \sim U(0, 1)$ into non-uniform random variates $X$ following an arbitrary continuous probability distribution with Cumulative Distribution Function (CDF) $F(x)$.

Theoretical Foundation (Probability Integral Transform): Let $X$ be a continuous random variable with strictly increasing CDF $F(x) = P(X \le x)$. The variable $U = F(X)$ is uniformly distributed on the interval $(0, 1)$. Therefore, setting $U = F(X)$ and solving for $X$ yields:

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

Step-by-Step Derivation for Exponential Distribution

1. Define CDF: $F(x) = 1 - e^{-\lambda x}$ for $x \ge 0$.
2. Set $F(X) = U$: $1 - e^{-\lambda X} = U$
3. Solve algebraically for $X$: $e^{-\lambda X} = 1 - U \implies -\lambda X = \ln(1 - U) \implies X = -\frac{1}{\lambda} \ln(1 - U)$
4. Since $(1 - U)$ is also uniformly distributed on $(0, 1)$, the formulation is commonly written as: $X = -\frac{1}{\lambda} \ln(U)$

Common Inverse Transform Formulations

Other Random Variate Generation Methods

• Acceptance-Rejection Method: Used when the CDF cannot be analytically inverted (e.g., Gamma, Beta distributions). Variates are sampled from an easy-to-generate proposal distribution and accepted with probability proportional to the ratio of the target density to the proposal density.
• Box-Muller Transformation (for Normal Distributions): Generates two independent standard normal variates $Z_1, Z_2 \sim N(0, 1)$ from two independent uniform numbers $U_1, U_2 \sim U(0, 1)$: $Z_1 = \sqrt{-2 \ln U_1} \cos(2\pi U_2), \quad Z_2 = \sqrt{-2 \ln U_1} \sin(2\pi U_2)$ To obtain $X \sim N(\mu, \sigma^2)$, scale by $X = \mu + \sigma Z_1$.

---

Simulation Experiment Design: Terminating vs. Steady-State Systems

Simulation models are classified by their operating time horizon:

1. Terminating (Transient) Simulations

A terminating simulation possesses a natural, physical event that defines the start and end of the system's operational cycle:

• Examples: A bank branch open from 9:00 AM to 5:00 PM; a manufacturing facility processing a discrete work order of 500 units; an emergency room during a 24-hour hurricane event.
• The initial conditions (e.g., empty facility at opening) are a realistic part of the operational system.
• Output analysis focuses on transient performance measures across $R$ independent replications, each beginning at identical initial conditions.

2. Non-Terminating (Steady-State) Simulations

A non-terminating simulation models a system that runs continuously without a defined stopping point:

• Examples: A 24/7 semiconductor fabrication plant; an automated distribution hub; an internet telecommunication router.
• The objective is to estimate long-term steady-state performance parameters: $\theta = \lim_{t \to \infty} E[Y(t)]$

The Initialization Bias & Welch's Graphical Procedure

When a steady-state simulation starts from an empty-and-idle state, early observations are biased downward (shorter queues, zero initial waiting times). Including these initial data points introduces initialization bias.

To eliminate initialization bias, the model must run through a warm-up period ($T_0$), during which all generated statistical observations are deleted:

Welch's Graphical Procedure for Determining Warm-Up Period:

1. Conduct $R$ independent replications of the simulation, each of length $m$ observations ($R \ge 5$).
2. Compute the cross-replication average for each discrete observation index $i = 1, 2, \dots, m$: $\bar{Y}_i = \frac{1}{R} \sum_{r=1}^R Y_{ri}$
3. Smooth the cross-replication averages using a symmetric moving average with window half-width $w$: $\bar{Y}_i(w) = \begin{cases} \frac{1}{2w + 1} \sum_{s=-w}^w \bar{Y}_{i+s} & \text{for } i = w+1, \dots, m-w \\ \frac{1}{2i - 1} \sum_{s=-(i-1)}^{i-1} \bar{Y}_{i+s} & \text{for } i = 1, \dots, w \end{cases}$
4. Plot $\bar{Y}_i(w)$ against observation index $i$. The warm-up truncation point $d$ is identified visually as the point beyond which the moving average curve stabilizes into horizontal steady-state oscillations. Delete all data for $i \le d$.

---

Statistical Output Analysis & Replication Framework

Because simulation outputs are stochastic random variables, a single simulation run represents only a single sample point ($N = 1$). Valid engineering conclusions require statistical confidence intervals.

Method of Independent Replications

To obtain statistically independent, identically distributed (i.i.d.) observations:

1. Execute $R$ independent replications of the simulation model, resetting the system state but initializing each replication with an independent pseudo-random number seed.
2. For each replication $r$, calculate the mean performance measure $\bar{X}_r$ after deleting the warm-up period: $\bar{X}_r = \frac{1}{n - d} \sum_{i=d+1}^n X_{ri}$
3. Compute the grand sample mean across all $R$ replications: $\bar{\bar{X}} = \frac{1}{R} \sum_{r=1}^R \bar{X}_r$
4. Compute the sample variance between replication means: $S^2 = \frac{1}{R - 1} \sum_{r=1}^R (\bar{X}_r - \bar{\bar{X}})^2$
5. Construct the $(1 - \alpha)100\%$ Confidence Interval for the true mean $\mu_X$ using the Student's $t$-distribution with $R - 1$ degrees of freedom:

$\bar{\bar{X}} \pm t_{\alpha/2, R-1} \cdot \frac{S}{\sqrt{R}}$

Where the term $HW = t_{\alpha/2, R-1} \cdot \frac{S}{\sqrt{R}}$ is the half-width of the confidence interval.

Sample Size Determination

If an engineering study requires the confidence interval half-width to be no larger than an acceptable error margin $\epsilon$, the required number of replications $R^*$ is estimated by:

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

---

Model Credibility: Verification, Validation & Testing Protocols

A simulation model is useless unless stakeholders trust its decisions. Credibility requires rigorous separation of verification and validation.

Verification vs. Validation Matrix

Verification Techniques

• Structured Walkthroughs: Peer review of simulation model logic by an independent engineering team.
• Trace Analysis: Step-by-step printing of every state transition, FEL change, and attribute update over a small time window.
• Animation Inspection: Visual review of dynamic 2D/3D graphics to detect bottlenecks, buffer overflows, or unexpected resource deadlocks.
• Conservation of Flow Audits: Ensuring that over a closed run, total entities entering equals total entities exiting plus entities remaining in queues and service.
• Extreme Condition / Stress Testing: Setting arrival rate $\lambda \to 0$ (system should be completely idle) or $\lambda > s\mu$ (queues must explode systematically).

Validation Protocols

• Face Validity: Subject Matter Experts (SMEs, shop floor supervisors, line operators) review model input parameters and observe animations to verify reasonable operational behavior.
• Historical Data Validation: Drive the simulation model using actual historical input traces (e.g., recorded arrival times from last month) and test whether model performance measures match historical outputs.
• Statistical Goodness-of-Fit Tests: Validating input probability distributions using:
  • Chi-Square Goodness-of-Fit Test ($\chi^2$): For discrete or continuous data with large sample sizes ($n \ge 50$).
  • Kolmogorov-Smirnov Test ($K$-$S$): For continuous empirical data, measuring the maximum vertical deviation $D = \max |F_n(x) - F_0(x)|$ between empirical and theoretical CDFs.
• The Turing Test: Present operational output reports from both the real physical facility and the simulation model to experienced managers in an identical format. If the managers cannot distinguish between real and simulated data, the model has demonstrated high operational validity.