9.3 Discrete-Event Simulation Principles and Monte Carlo Modeling

FE Industrial and Systems Study Guide Exam

Track your progress as you read

Finished sections get a check mark.

1 Chapter 1: Exam Overview, CBT Strategy & Handbook Navigation

1.1 FE Industrial Exam Structure, CBT Logistics, and Pacing 1.2 Using the NCEES FE Reference Handbook and Approved Calculators

2 Chapter 2: Mathematics for Industrial Engineers

2.1 Analytic Geometry, Coordinate Systems, and Trigonometry 2.2 Calculus: Differential, Integral, and Series Approximations 2.3 Linear Algebra: Matrix Operations, Systems of Equations, and Vectors

3 Chapter 3: Engineering Sciences Core

3.1 Thermodynamics Fundamentals and Fluid Flow Principles 3.2 Statics, Dynamics, and Mechanics of Materials Core 3.3 Electrical Circuits, DC/AC Analysis, and Power in Industrial Systems

4 Chapter 4: Ethics and Professional Practice

4.1 NCEES Model Rules of Professional Conduct and Technical Society Ethics 4.2 Engineering Licensure, Public Welfare, and Professional Liability 4.3 Engineering Contracts, Intellectual Property, and Conflicts of Interest

5 Chapter 5: Engineering Economics

5.1 Time Value of Money, Discount Factors, and Equivalence 5.2 Evaluating Alternatives: Present Worth, Equivalent Annual Cost, and IRR 5.3 Incremental Rate of Return, Replacement Analysis, and Break-Even 5.4 Depreciation Methods, MACRS, Income Taxes, and Inflation

6 Chapter 6: Probability Foundations & Distributions

6.1 Basic Probability Concepts, Conditional Probability, and Bayes' Theorem 6.2 Discrete Probability Distributions 6.3 Continuous Probability Distributions

7 Chapter 7: Statistical Inference, Regression & Design of Experiments

7.1 Sampling Distributions, Central Limit Theorem, and Confidence Intervals 7.2 Hypothesis Testing: Means, Proportions, and Variances 7.3 Linear Regression, Correlation Analysis, and Design of Experiments (ANOVA)

8 Chapter 8: Data, Logic, Modeling & Linear Programming

8.1 Optimization Modeling: Decision Variables, Objectives, and Constraints 8.2 Linear Programming: Graphical Solutions and Simplex Fundamentals 8.3 Duality, Sensitivity Analysis, and Integer Linear Programming 8.4 Data, Logic Development, and Analytics

9 Chapter 9: Stochastic Models & Simulation

9.1 Queuing Theory: Little's Law, M/M/1, and Multi-Server Systems 9.2 Discrete-Time Markov Chains: Transition Matrices and Steady-State Probabilities 9.3 Discrete-Event Simulation Principles and Monte Carlo Modeling

10 Chapter 10: Engineering Management, Project Planning & Performance Measurement

10.1 Organizational Structure, Leadership, and Worker Motivation Theories 10.2 Project Management: WBS, Critical Path Method (CPM), PERT, and Crashing 10.3 Earned Value Management (EVM), Decision Trees, and Risk Preference 10.4 Performance Measurement: Productivity, KPIs, Wage Systems, and the Balanced Scorecard

11 Chapter 11: Manufacturing & Production Systems

11.1 Manufacturing Processes, Material Processing, and Tolerancing 11.2 Assembly Line Balancing, Takt Time, Cycle Time, and Throughput 11.3 Sustainable Manufacturing, Waste Minimization, and Life-Cycle Assessment

12 Chapter 12: Production Planning, Inventory Control & Sequencing

12.1 Demand Forecasting Methods and Forecast Error Tracking 12.2 Aggregate Planning, Master Production Scheduling (MPS), and MRP 12.3 Deterministic Inventory Models: EOQ, EPQ, and Safety Stock 12.4 Job Sequencing, Dispatching Rules, and Flow Shop Scheduling

13 Chapter 13: Lean Systems & Theory of Constraints

13.1 Lean Foundations: Value Stream Mapping, Seven Wastes, and 5S 13.2 Pull Systems, Kanban, Just-In-Time (JIT), and SMED Setup Reduction 13.3 Theory of Constraints: Bottleneck Identification and Drum-Buffer-Rope

14 Chapter 14: Facilities Planning & Material Handling

14.1 Facility Location Analysis: Qualitative Factor Rating and Centroid Method 14.2 Facility Layout Design and Flow Analysis 14.3 Material Handling Equipment, Warehousing, and AS/RS

15 Chapter 15: Supply Chain Design & Logistics

15.1 Supply Chain Network Architecture and Transportation Modes 15.2 Distribution Strategies, Cross-Docking, and Vehicle Routing 15.3 The Bullwhip Effect, Risk Pooling, and Supply Chain Coordination

16 Chapter 16: Human Factors & Cognitive Ergonomics

16.1 Cognitive Engineering, Human-Machine Interface, Displays, and Controls 16.2 Anthropometric Principles, Percentile Design, and Biomechanical Posture 16.3 Environmental Ergonomics: Industrial Noise, Lighting, and Heat Stress

17 Chapter 17: Physical Ergonomics & Workplace Safety

17.1 Physical Ergonomics: The NIOSH Manual Lifting Equation 17.2 Cumulative Trauma Disorders and Postural Assessment 17.3 Workplace Safety: OSHA Standards, Hazard Analysis, and Control Hierarchy

18 Chapter 18: Work Design & Measurement

18.1 Methods Engineering, Process Charting, and Principles of Motion Economy 18.2 Direct Stopwatch Time Study: Observed Time, Performance Rating, and Allowances 18.3 Predetermined Time Systems, Work Sampling, and Learning Curves

19 Chapter 19: Quality Management, Six Sigma & Analytical Tools

19.1 Quality Philosophies, Total Quality Management, and Cost of Quality 19.2 Six Sigma DMAIC Methodology and Quality Function Deployment 19.3 The Seven Basic Quality Tools, Root Cause Analysis, and Taguchi Loss Function

20 Chapter 20: Statistical Process Control

20.1 Statistical Process Control Principles, Variation, and Run Rules 20.2 Variable Control Charts: X-bar & R, X-bar & S, and Individual-Moving Range 20.3 Attribute Control Charts: Defectives (p, np) and Defects (c, u)

21 Chapter 21: Process Capability & Acceptance Sampling

21.1 Process Capability Indices and Defect Levels 21.2 Acceptance Sampling Concepts and Operating Characteristic Curves 21.3 Sampling Risks, AQL, LTPD, and Standard Sampling Plans

22 Chapter 22: Systems Engineering & Life-Cycle Management

22.1 The Systems Engineering V-Model, Stakeholder Requirements, and Functional Architecture 22.2 Configuration Management, Interface Control, and Verification & Validation 22.3 Life-Cycle Costing, Maintainability, and System Availability

23 Chapter 23: Reliability Engineering & Risk Analysis

23.1 Reliability Mathematics, Hazard Rates, and the Bathtub Curve 23.2 System Reliability Configurations: Series, Parallel, and k-out-of-n Systems 23.3 Failure Modes and Effects Analysis and Fault Tree Analysis

9.3 Discrete-Event Simulation Principles and Monte Carlo Modeling

Key Takeaways Explain This Topic AI

Discrete-Event Simulation (DES) models dynamic systems where state variables change instantaneously at distinct event occurrences, utilizing a next-event time-advance mechanism governed by a Future Event List (FEL), contrasting with continuous simulation driven by differential equations.

The foundational components of a DES architecture comprise entities (dynamic flowing objects with local attributes), resources (capacity-constrained service assets with states like idle, busy, blocked), queues, and events.

Linear Congruential Generators (LCG: X_{i+1} = (aX_i + c) mod m) produce pseudo-random integers that are normalized to U(0, 1) random numbers; full period m requires meeting the Hull-Dobell criteria.

The Inverse Transform Method generates non-uniform random variates by applying the inverse cumulative distribution function to uniform random numbers (X = F^{-1}(U)), specifically generating exponential variates via X = -(1/λ)ln(U).

Rigorous statistical output analysis requires identifying and truncating the initial warm-up period to eliminate startup initialization bias, executing independent replications with different seeds, and calculating confidence intervals via Student's t-distribution.

Last updated: September 2026

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 Δ t \Delta t Δ 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 Δ t \Delta t Δ t . At each tick t + Δ t t + \Delta t t + Δ 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:

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.

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 W W ), required processing recipe, physical weight, customer priority, or rework status.

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.

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

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.

Simulation Clock and Future Event List (FEL) :

Clock ( T T 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 ≤ t 2 ≤ t 3 ≤ … t_1 \le t_2 \le t_3 \le \dots t 1 ​ ≤ t 2 ​ ≤ t 3 ​ ≤ … ). 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 ) (0, 1) ( 0 , 1 ) , denoted U ∼ U ( 0 , 1 ) U \sim U(0, 1) U ∼ U ( 0 , 1 ) . These numbers must be:

Uniform : Equal probability of falling into any sub-interval of ( 0 , 1 ) (0, 1) ( 0 , 1 ) .

Independent : No serial correlation between successive numbers.

Reproducible : Given an identical initial starting integer ( seed X 0 X_0 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 ) ( m o d m ) X_{i+1} = (a X_i + c) \pmod m X i + 1 ​ = ( a X i ​ + c ) ( mod m )

Where:

• X 0 X_0 X 0 ​ : Initial starting value (the seed , X 0 ≥ 0 X_0 \ge 0 X 0 ​ ≥ 0 )
• a a a : Multiplier ( a > 0 a > 0 a > 0 )

c c c : Increment ( c ≥ 0 c \ge 0 c ≥ 0 ); if c = 0 c = 0 c = 0 , it is called a multiplicative LCG; if c > 0 c > 0 c > 0 , a mixed LCG

• m m m : Modulus ( m > X 0 , a , c m > X_0, a, c m > X 0 ​ , a , c )
• ( m o d m ) \pmod m ( mod m ) : The remainder after integer division by m m m

To map the resulting integer X i X_i X i ​ to a uniform continuous random number on [ 0 , 1 ) [0, 1) [ 0 , 1 ) : U i = X i m U_i = \frac{X_i}{m} U i ​ = m X i ​ ​

Full-Period Conditions (The Hull-Dobell Theorem)

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

c c c and m m m are relatively prime (their greatest common divisor is 1 1 1 : gcd ⁡ ( c , m ) = 1 \gcd(c, m) = 1 g cd ( c , m ) = 1 ).

a − 1 a - 1 a − 1 is divisible by all prime factors of m m m .

a − 1 a - 1 a − 1 is divisible by 4 4 4 if m m m is divisible by 4 4 4 .

4. Random Variate Generation: The Inverse Transform Method

Once uniform random numbers U ∼ U ( 0 , 1 ) U \sim U(0, 1) U ∼ 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 X X be a continuous random variable with cumulative distribution function F ( x ) = P ( X ≤ x ) F(x) = P(X \le x) F ( x ) = P ( X ≤ x ) . The Inverse Transform Method is based on the mathematical theorem that if U ∼ U ( 0 , 1 ) U \sim U(0, 1) U ∼ U ( 0 , 1 ) , then the random variable: X = F − 1 ( U ) X = F^{-1}(U) X = F − 1 ( U )

• follows the distribution defined by F ( x ) F(x) 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 β = 1 / λ \beta = 1/\lambda β = 1/ λ ) has cumulative distribution function: F ( x ) = 1 − e − λ x , x ≥ 0 F(x) = 1 - e^{-\lambda x}, \quad x \ge 0 F ( x ) = 1 − e − λ x , x ≥ 0

Setting F ( X ) = U F(X) = U F ( X ) = U and solving for X X X : 1 − e − λ X = U 1 - e^{-\lambda X} = U 1 − e − λ X = U ≠ − λ X = 1 − U \ne^{-\lambda X} = 1 - U  = − λ X = 1 − U − λ X = ln ⁡ ( 1 − U ) -\lambda X = \ln(1 - U) − λ X = ln ( 1 − U ) X = − 1 λ ln ⁡ ( 1 − U ) X = -\frac{1}{\lambda} \ln(1 - U) X = − λ 1 ​ ln ( 1 − U )

Because U U U is uniformly distributed on ( 0 , 1 ) (0, 1) ( 0 , 1 ) , its complement 1 − U 1 - U 1 − U is also identically distributed as U ( 0 , 1 ) U(0, 1) U ( 0 , 1 ) . Therefore, the standard computational formula is: X = − 1 λ ln ⁡ ( U ) = − β ln ⁡ ( U ) X = -\frac{1}{\lambda} \ln(U) = -\beta \ln(U) X = − λ 1 ​ ln ( U ) = − β ln ( U )

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

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

The CDF of a uniform variable on the interval [ a , b ] [a, b] [ a , b ] is F ( x ) = x − a b − a F(x) = \frac{x - a}{b - a} F ( x ) = b − a x − a ​ . Setting F ( X ) = U F(X) = U F ( X ) = U yields: X = a + ( b − a ) U X = a + (b - a)U 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 ) p(x_i) = P(X = x_i) p ( x i ​ ) = P ( X = x i ​ ) , the cumulative distribution F ( x k ) = ∑ i = 1 k p ( x i ) F(x_k) = \sum_{i=1}^k p(x_i) F ( x k ​ ) = ∑ i = 1 k ​ p ( x i ​ ) partitions the interval [ 0 , 1 ) [0, 1) [ 0 , 1 ) into sub-intervals. If a generated uniform number U U U falls within [ F ( x k − 1 ) , F ( x k ) ) [F(x_{k-1}), F(x_k)) [ F ( x k − 1 ​ ) , F ( x k ​ )) , the variate X = x k X = x_k 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

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.

Mechanical Tolerance Stack-Up : Simulating dimensional tolerance accumulation across multiple mated parts in complex assemblies to determine expected assembly defect rates.

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 n = 1 n = 1 ) from an underlying output population. Rigorous statistical analysis is required to draw valid engineering conclusions.

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

Classification Physical Characteristics Starting / Stopping Conditions Statistical Analysis Approach

Terminating System System has a natural, fixed starting condition and a distinct terminating event Opens empty at fixed time t 0 t_0 t 0 ​ , closes at defined time T E T_E T E ​ or upon batch completion Execute R R R independent replications from time 0 0 0 to T E T_E T E ​ ; analyze across replications

Non-Terminating (Steady-State) System operates continuously without a predefined endpoint (e.g., 24/7 semiconductor fab) Runs indefinitely; interested in long-run stationary equilibrium behavior Discard warm-up period to remove initialization bias; execute long runs or multiple replications

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 T_w T w ​ until transient startup effects dissipate and steady-state conditions are reached. All data collected during [ 0 , T w ] [0, T_w] [ 0 , T w ​ ] is discarded (truncated), and statistics are accumulated strictly over [ T w , T final ] [T_w, T_{\text{final}}] [ T w ​ , T 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 T_w 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 μ Y \mu_Y μ Y ​ :

Execute R R R independent simulation runs ( replications ), each using different random number seeds and starting from identical initial conditions.

Let Y r Y_r Y r ​ be the sample mean output from replication r r r ( r = 1 , 2 , … , R r = 1, 2, \dots, R r = 1 , 2 , … , R ). The point estimator of the overall mean is: Y ˉ = 1 R ∑ r = 1 R Y r \bar{Y} = \frac{1}{R} \sum_{r=1}^R Y_r Y ˉ = R 1 ​ r = 1 ∑ R ​ Y r ​

The sample variance across the R R R independent replications is: S 2 = 1 R − 1 ∑ r = 1 R ( Y r − Y ˉ ) 2 S^2 = \frac{1}{R - 1} \sum_{r=1}^R (Y_r - \bar{Y})^2 S 2 = R − 1 1 ​ r = 1 ∑ R ​ ( Y r ​ − Y ˉ ) 2

A ( 1 − α ) 100 % (1 - \alpha)100\% ( 1 − α ) 100% confidence interval for the true mean μ Y \mu_Y μ Y ​ is computed using the Student's t t t -distribution with ν = R − 1 \nu = R - 1 ν = R − 1 degrees of freedom: Y ˉ ± t α / 2 , R − 1 S R \bar{Y} \pm t_{\alpha/2, R-1} \frac{S}{\sqrt{R}} Y ˉ ± t α /2 , R − 1 ​ R ​ S ​

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

If an engineer requires an estimate with a specified maximum allowable half-width error ϵ \epsilon ϵ at confidence level ( 1 − α ) (1 - \alpha) ( 1 − α ) , the required total number of replications is approximated by: R ∗ ≥ ( t α / 2 , R − 1 ⋅ S ϵ ) 2 ≈ ( z α / 2 ⋅ S ϵ ) 2 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 R ∗ ≥ ( ϵ t α /2 , R − 1 ​ ⋅ S ​ ) 2 ≈ ( ϵ z α /2 ​ ⋅ S ​ ) 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: Var ( Y A − Y B ) = Var ( Y A ) + Var ( Y B ) − 2 Cov ( Y A , Y B ) \text{Var}(Y_A - Y_B) = \text{Var}(Y_A) + \text{Var}(Y_B) - 2 \text{Cov}(Y_A, Y_B) Var ( Y A ​ − Y B ​ ) = Var ( Y A ​ ) + Var ( Y B ​ ) − 2 Cov ( Y A ​ , Y B ​ ) Positive correlation ( Cov ( Y A , Y B ) > 0 \text{Cov}(Y_A, Y_B) > 0 Cov ( Y A ​ , Y B ​ ) > 0 ) significantly reduces Var ( Y A − Y B ) \text{Var}(Y_A - Y_B) 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 , … U_1, U_2, \dots U 1 ​ , U 2 ​ , … , a complementary paired run is driven by 1 − U 1 , 1 − U 2 , … 1 - U_1, 1 - U_2, \dots 1 − U 1 ​ , 1 − U 2 ​ , … , inducing negative covariance ( Cov ( Y 1 , Y 2 ) < 0 \text{Cov}(Y_1, Y_2) < 0 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 :

A mixed Linear Congruential Generator is defined by: X i + 1 = ( 7 X i + 3 ) ( m o d 16 ) X_{i+1} = (7 X_i + 3) \pmod{16} X i + 1 ​ = ( 7 X i ​ + 3 ) ( mod 16 ) Using a starting seed X 0 = 5 X_0 = 5 X 0 ​ = 5 , compute the next three integer values ( X 1 , X 2 , X 3 X_1, X_2, X_3 X 1 ​ , X 2 ​ , X 3 ​ ) and their associated uniform random variates ( U 1 , U 2 , U 3 U_1, U_2, U_3 U 1 ​ , U 2 ​ , U 3 ​ ).

Use the Inverse Transform Method to generate an interarrival time following an Exponential distribution with mean arrival rate λ = 4 \lambda = 4 λ = 4 parts per hour, using the generated uniform number U 1 U_1 U 1 ​ .

Solution :

Compute LCG Sequence :

Step 1 ( i = 0 → 1 i = 0 \to 1 i = 0 → 1 ) : X 1 = ( 7 X 0 + 3 ) ( m o d 16 ) = ( 7 × 5 + 3 ) ( m o d 16 ) = ( 35 + 3 ) ( m o d 16 ) = 38 ( m o d 16 ) X_1 = (7 X_0 + 3) \pmod{16} = (7 \times 5 + 3) \pmod{16} = (35 + 3) \pmod{16} = 38 \pmod{16} X 1 ​ = ( 7 X 0 ​ + 3 ) ( mod 16 ) = ( 7 × 5 + 3 ) ( mod 16 ) = ( 35 + 3 ) ( mod 16 ) = 38 ( mod 16 ) 38 / 16 = 2 with a remainder of 6 ⟹ X 1 = 6 38 / 16 = 2 \text{ with a remainder of } 6 \implies X_1 = 6 38/16 = 2 with a remainder of 6 ⟹ X 1 ​ = 6 U 1 = X 1 m = 6 16 = 0.3750 U_1 = \frac{X_1}{m} = \frac{6}{16} = 0.3750 U 1 ​ = m X 1 ​ ​ = 16 6 ​ = 0.3750

Step 2 ( i = 1 → 2 i = 1 \to 2 i = 1 → 2 ) : X 2 = ( 7 X 1 + 3 ) ( m o d 16 ) = ( 7 × 6 + 3 ) ( m o d 16 ) = ( 42 + 3 ) ( m o d 16 ) = 45 ( m o d 16 ) X_2 = (7 X_1 + 3) \pmod{16} = (7 \times 6 + 3) \pmod{16} = (42 + 3) \pmod{16} = 45 \pmod{16} X 2 ​ = ( 7 X 1 ​ + 3 ) ( mod 16 ) = ( 7 × 6 + 3 ) ( mod 16 ) = ( 42 + 3 ) ( mod 16 ) = 45 ( mod 16 ) 45 / 16 = 2 with a remainder of 13 ⟹ X 2 = 13 45 / 16 = 2 \text{ with a remainder of } 13 \implies X_2 = 13 45/16 = 2 with a remainder of 13 ⟹ X 2 ​ = 13 U 2 = X 2 m = 13 16 = 0.8125 U_2 = \frac{X_2}{m} = \frac{13}{16} = 0.8125 U 2 ​ = m X 2 ​ ​ = 16 13 ​ = 0.8125

Step 3 ( i = 2 → 3 i = 2 \to 3 i = 2 → 3 ) : X 3 = ( 7 X 2 + 3 ) ( m o d 16 ) = ( 7 × 13 + 3 ) ( m o d 16 ) = ( 91 + 3 ) ( m o d 16 ) = 94 ( m o d 16 ) X_3 = (7 X_2 + 3) \pmod{16} = (7 \times 13 + 3) \pmod{16} = (91 + 3) \pmod{16} = 94 \pmod{16} X 3 ​ = ( 7 X 2 ​ + 3 ) ( mod 16 ) = ( 7 × 13 + 3 ) ( mod 16 ) = ( 91 + 3 ) ( mod 16 ) = 94 ( mod 16 ) 94 / 16 = 5 with a remainder of 14 ⟹ X 3 = 14 94 / 16 = 5 \text{ with a remainder of } 14 \implies X_3 = 14 94/16 = 5 with a remainder of 14 ⟹ X 3 ​ = 14 U 3 = X 3 m = 14 16 = 0.8750 U_3 = \frac{X_3}{m} = \frac{14}{16} = 0.8750 U 3 ​ = m X 3 ​ ​ = 16 14 ​ = 0.8750

Generate Exponential Variate via Inverse Transform :

Arrival rate: λ = 4 \lambda = 4 λ = 4 parts/hour

Mean interarrival time: β = 1 / λ = 1 / 4 = 0.25 \beta = 1/\lambda = 1/4 = 0.25 β = 1/ λ = 1/4 = 0.25 hours = 15 = 15 = 15 minutes

Apply inverse transform formula with U 1 = 0.3750 U_1 = 0.3750 U 1 ​ = 0.3750 : X = − 1 λ ln ⁡ ( U 1 ) = − 0.25 × ln ⁡ ( 0.3750 ) X = -\frac{1}{\lambda} \ln(U_1) = -0.25 \times \ln(0.3750) X = − λ 1 ​ ln ( U 1 ​ ) = − 0.25 × ln ( 0.3750 ) ln ⁡ ( 0.3750 ) ≈ − 0.98083 \ln(0.3750) \approx -0.98083 ln ( 0.3750 ) ≈ − 0.98083 X = − 0.25 × ( − 0.98083 ) = 0.2452 hours X = -0.25 \times (-0.98083) = 0.2452 \text{ hours} X = − 0.25 × ( − 0.98083 ) = 0.2452 hours

In minutes: 0.2452 × 60 = 14.71 0.2452 \times 60 = 14.71 0.2452 × 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 R = 5 R = 5 independent replications. The observed average hourly throughputs (totes/hour) across the 5 replications are: Y = [ 124 , 130 , 122 , 136 , 128 ] Y = [124, 130, 122, 136, 128] Y = [ 124 , 130 , 122 , 136 , 128 ]

Calculate the point estimate Y ˉ \bar{Y} Y ˉ and sample standard deviation S S S .

Construct a 95 % 95\% 95% confidence interval for the true mean hourly throughput (for ν = 4 \nu = 4 ν = 4 degrees of freedom, t 0.025 , 4 = 2.776 t_{0.025, 4} = 2.776 t 0.025 , 4 ​ = 2.776 ).

If management requires the estimate to have a maximum half-width error of ϵ = 2.0 \epsilon = 2.0 ϵ = 2.0 totes/hour with 95 % 95\% 95% confidence ( z 0.025 = 1.96 z_{0.025} = 1.96 z 0.025 ​ = 1.96 ), how many total replications R ∗ R^* R ∗ are required?

Solution :

Point Estimator and Sample Variance :

Sample mean: Y ˉ = 124 + 130 + 122 + 136 + 128 5 = 640 5 = 128.0 totes/hour \bar{Y} = \frac{124 + 130 + 122 + 136 + 128}{5} = \frac{640}{5} = 128.0 \text{ totes/hour} Y ˉ = 5 124 + 130 + 122 + 136 + 128 ​ = 5 640 ​ = 128.0 totes/hour

Deviations about the mean: ( 124 − 128 ) 2 = ( − 4 ) 2 = 16 (124 - 128)^2 = (-4)^2 = 16 ( 124 − 128 ) 2 = ( − 4 ) 2 = 16 ( 130 − 128 ) 2 = ( 2 ) 2 = 4 (130 - 128)^2 = (2)^2 = 4 ( 130 − 128 ) 2 = ( 2 ) 2 = 4 ( 122 − 128 ) 2 = ( − 6 ) 2 = 36 (122 - 128)^2 = (-6)^2 = 36 ( 122 − 128 ) 2 = ( − 6 ) 2 = 36 ( 136 − 128 ) 2 = ( 8 ) 2 = 64 (136 - 128)^2 = (8)^2 = 64 ( 136 − 128 ) 2 = ( 8 ) 2 = 64 ( 128 − 128 ) 2 = ( 0 ) 2 = 0 (128 - 128)^2 = (0)^2 = 0 ( 128 − 128 ) 2 = ( 0 ) 2 = 0 ∑ ( Y r − Y ˉ ) 2 = 16 + 4 + 36 + 64 + 0 = 120 \sum (Y_r - \bar{Y})^2 = 16 + 4 + 36 + 64 + 0 = 120 ∑ ( Y r ​ − Y ˉ ) 2 = 16 + 4 + 36 + 64 + 0 = 120

Sample variance and standard deviation: S 2 = 120 5 − 1 = 120 4 = 30.0 ⟹ S = 30.0 ≈ 5.477 totes/hour S^2 = \frac{120}{5 - 1} = \frac{120}{4} = 30.0 \implies S = \sqrt{30.0} \approx 5.477 \text{ totes/hour} S 2 = 5 − 1 120 ​ = 4 120 ​ = 30.0 ⟹ S = 30.0 ​ ≈ 5.477 totes/hour

Construct 95 % 95\% 95% Confidence Interval :

Standard error of the mean: SE = S R = 5.477 5 = 5.477 2.2361 ≈ 2.449 totes/hour \text{SE} = \frac{S}{\sqrt{R}} = \frac{5.477}{\sqrt{5}} = \frac{5.477}{2.2361} \approx 2.449 \text{ totes/hour} SE = R ​ S ​ = 5 ​ 5.477 ​ = 2.2361 5.477 ​ ≈ 2.449 totes/hour

Half-width error: H = t 0.025 , 4 × SE = 2.776 × 2.449 ≈ 6.80 totes/hour H = t_{0.025, 4} \times \text{SE} = 2.776 \times 2.449 \approx 6.80 \text{ totes/hour} H = t 0.025 , 4 ​ × SE = 2.776 × 2.449 ≈ 6.80 totes/hour

95 % 95\% 95% Confidence Interval: 128.0 ± 6.80 ⟹ [ 121.20 , 134.80 ] totes/hour 128.0 \pm 6.80 \implies [121.20, \; 134.80] \text{ totes/hour} 128.0 ± 6.80 ⟹ [ 121.20 , 134.80 ] totes/hour

Determine Required Replications ( R ∗ R^* R ∗ ) :

Target half-width ϵ = 2.0 \epsilon = 2.0 ϵ = 2.0 totes/hour: R ∗ ≥ ( z 0.025 ⋅ S ϵ ) 2 = ( 1.96 × 5.477 2.0 ) 2 = ( 10.735 2.0 ) 2 = ( 5.368 ) 2 ≈ 28.81 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 R ∗ ≥ ( ϵ z 0.025 ​ ⋅ S ​ ) 2 = ( 2.0 1.96 × 5.477 ​ ) 2 = ( 2.0 10.735 ​ ) 2 = ( 5.368 ) 2 ≈ 28.81

Rounding up to the next integer yields R ∗ = 29 R^* = 29 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 / λ ) ln ⁡ ( U ) = − β ln ⁡ ( U ) X = -(1/\lambda)\ln(U) = -\beta \ln(U) X = − ( 1/ λ ) ln ( U ) = − β ln ( U ) . If given a mean service time of β = 10 \beta = 10 β = 10 minutes, do not write − 1 / 10 ln ⁡ ( U ) -1/10 \ln(U) − 1/10 ln ( U ) . The parameter λ \lambda λ is the rate ( 1 / 10 1/10 1/10 ), so 1 / λ = 10 1/\lambda = 10 1/ λ = 10 . Always check whether the problem specifies the rate λ \lambda λ or the mean β \beta β .

The Negative Sign in Exponential Generation : Because U ∈ ( 0 , 1 ) U \in (0, 1) U ∈ ( 0 , 1 ) , ln ⁡ ( U ) \ln(U) ln ( U ) is always a negative number . The leading minus sign in X = − ( 1 / λ ) ln ⁡ ( U ) X = -(1/\lambda)\ln(U) X = − ( 1/ λ ) ln ( U ) is mandatory to ensure the generated physical time X X 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 ( Y ˉ A − Y ˉ B ) (\bar{Y}_A - \bar{Y}_B) ( Y ˉ A ​ − 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 X_0 X 0 ​ will produce five identical runs, yielding zero sample variance and invalidating statistical conclusions.