Simulation And Modeling (ME-804) - Unit 3 Short Notes

UNIT 3: SIMULATION IMPLEMENTATION, VALIDATION, AND ANALYSIS

3.1. Random Number and Random Variate Generation

3.1.1. Properties of Good Random Number Generators (RNGs)

A high-quality RNG must satisfy:

• Uniformity: Numbers are equally likely over the interval.

• Independence: No correlation between successive numbers.

• Long Period: Sequence does not repeat for a very large number of calls.

• Repeatability: Same seed produces identical sequence (for debugging).

• Speed: Efficient computation.

• Portability: Consistent behavior across different hardware/software.

3.1.2. Pseudo-Random Number Generators (PRNGs)

Linear Congruential Generator (LCG):

The most common PRNG. Recurrence relation:

$$X_{n+1} = (a X_n + c) \mod m$$

where:

• $m$ = modulus (defines range $0$ to $m-1$)

• $a$ = multiplier

• $c$ = increment

• $$\displaystyle X_0 $$ = seed

Types:

• Mixed LCG: $c \neq 0$ (full period possible if parameters chosen correctly).

• Multiplicative LCG: $$\displaystyle c = 0 $$ (period at most $m/2$).

Combined Multiple Recursive Generators (CMRG):

Combine outputs of two or more MRGs to achieve longer periods and better statistical properties. Example: Marsaglia-Multicarry or L'Ecuyer's CMRG.

Seeding & Stream Management:

• Seeding: Initial value $$\displaystyle X_0 $$.

• Independent Streams: Use different seeds/parameters for separate simulation components (e.g., arrival times, service times) to avoid correlation.

• Non-overlapping Streams: Ensure sequences from different streams do not overlap (using leapfrog or sequence splitting).

[!TIP] Exam Focus: LCG formula and parameters are frequently tested. Know how to compute next number manually.

3.1.3. Testing Random Number Generators

Theoretical Tests:

• Spectral Test: For LCGs, examines distance between hyperplanes in multi-dimensional space. Detects lattice structure.

Empirical Tests (on generated sequences):

• Chi-Square Test: Groups numbers into $k$ equal intervals. Test statistic:

$$\chi^2 = \sum_{i=1}^{k} \frac{(O_i - E_i)^2}{E_i}$$

where $$\displaystyle O_i $$ = observed frequency, $$\displaystyle E_i = n/k $$ (expected). Compare to $$\displaystyle \chi^2_{k-1, \alpha} $$.

• Kolmogorov-Smirnov (K-S) Test: For uniformity. Compares empirical CDF $$\displaystyle S_n(x) $$ to theoretical $$\displaystyle F(x)=x $$. Statistic:

$$D = \sup_x |S_n(x) - F(x)|$$

Compare to critical value $$\displaystyle D_\alpha $$.

• Runs Test: Checks independence by counting runs (consecutive above/below median). Test statistic based on normal approximation.

3.1.4. Random Variate Generation from Non-Uniform Distributions

Inverse Transform Technique:

1. Generate $U \sim \text{Uniform}(0,1)$.

2. Compute $$\displaystyle X = F^{-1}(U) $$, where $F$ is target CDF.

3. $X$ follows desired distribution.

Applications:

• Exponential ($\lambda$): $$\displaystyle X = -\frac{1}{\lambda} \ln(1-U) $$

• Uniform $(a,b)$: $$\displaystyle X = a + (b-a)U $$

• Weibull ($\beta, \eta$): $$\displaystyle X = \eta (-\ln(1-U))^{1/\beta} $$

[!TIP] Common Pitfall: For exponential, use $-\ln(U)$ since $1-U$ is also uniform.

Acceptance-Rejection Technique:

1. Find PDF $f(x)$ and bounding function $c g(x)$ where $c g(x) \ge f(x)$ for all $x$, and $g(x)$ is easy to sample (e.g., uniform).

2. Generate $Y \sim g(y)$, $U \sim \text{Uniform}(0,1)$.

3. If $$\displaystyle U \le \frac{f(Y)}{c g(Y)} $$, accept $$\displaystyle X=Y $$; else reject and repeat.

4. Efficiency $\approx 1/c$.

Convolution Method:

Sum of independent random variables. Example:

• Erlang ($k$ exponential with rate $\lambda$): Sum $k$ i.i.d. exponentials.

• Binomial ($n,p$): Sum $n$ i.i.d. Bernoulli trials.

Specialized Techniques:

• Normal: Box-Muller transform: $$\displaystyle Z_1 = \sqrt{-2\ln U_1} \cos(2\pi U_2) $$, $$\displaystyle Z_2 = \sqrt{-2\ln U_1} \sin(2\pi U_2) $$.

• Lognormal: Generate normal $Y$, then $$\displaystyle X = e^Y $$.

• Empirical: Use inverse transform on empirical CDF from data.

3.2. Input Data Analysis

3.2.1. Identifying Input Processes & Data Collection

• Stationary Process: Statistical properties (mean, variance) constant over time.

• Non-Stationary Process: Properties change over time (e.g., time-dependent arrival rates). Requires careful modeling (e.g., time-varying distributions).

• Data Collection: Ensure data represents the system under study. Avoid bias (e.g., collect during typical operations, not only peak hours).

3.2.2. Fitting Probability Distributions to Data

Qualitative Identification:

• Histogram: Visual shape (symmetry, tail behavior).

• Q-Q Plot: Compare quantiles of data to theoretical distribution. Linear fit indicates good fit.

Quantitative Goodness-of-Fit Tests:

• Chi-Squared Test (for grouped data):

$$\chi^2 = \sum_{i=1}^{k} \frac{(O_i - E_i)^2}{E_i}$$

$$\displaystyle E_i = n \cdot p_i $$ (expected count from fitted distribution). Requires $$\displaystyle E_i \ge 5 $$ typically. Reject if $$\displaystyle \chi^2 > \chi^2_{k-1-c, \alpha} $$, where $c$ = number of estimated parameters.

• Kolmogorov-Smirnov (K-S) Test (for continuous distributions):

$$D = \sup_x |F_n(x) - F(x)|$$

Compare to critical value $$\displaystyle D_\alpha $$ (depends on $n$). More powerful for fully specified distributions (no estimated parameters).

• Anderson-Darling Test: More sensitive to tail deviations. Statistic:

$$A^2 = -n - \frac{1}{n} \sum_{i=1}^{n} (2i-1) \left[ \ln F(x_i) + \ln(1-F(x_{n+1-i})) \right]$$

Parameter Estimation:

• Method of Moments (MOM): Equate sample moments to theoretical moments. Solve for parameters.

• Maximum Likelihood Estimation (MLE): Maximize likelihood function $$\displaystyle L(\theta) = \prod f(x_i|\theta) $$. Often requires numerical optimization.

3.2.3. Selecting the "Best" Distribution

• Use p-value from goodness-of-fit tests: high p-value (>0.05) suggests cannot reject fit.

• Consider process knowledge (e.g., interarrival times often exponential).

• Examine tail behavior (critical for queueing systems).

• Use software (e.g., ExpertFit, @RISK) that automates fitting and compares multiple distributions.

[!TIP] Exam Trap: K-S test is for continuous distributions only. Chi-squared can be used for discrete too.

3.3. Model Verification and Validation (V&V)

3.3.1. Verification: "Are we building the model right?"

• Debugging: Modular code, trace debugging (print state at key events), output tracing.

• Logic Checks: Ensure model implements specifications correctly (e.g., entity flow, resource scheduling).

3.3.2. Validation: "Are we building the right model?"

• Face Validation: Domain experts review model logic and outputs for reasonableness.

• Input-Output Validation: Perform sensitivity analysis; check if outputs respond plausibly to input changes.

• Model Structure Validation: Compare with simpler, analytically solvable model or validated model.

• Historical Data Validation: Compare model outputs to real system historical data (if available).

• Traces & Animation: Step-by-step comparison of model and actual system behavior.

3.3.3. The V&V Process

• Iterative: Performed throughout model development, not just at end.

• Documentation: Record all V&V activities, results, and changes made.

[!TIP] Key Distinction: Verification = correct implementation; Validation = correct representation of reality.

3.4. Output Analysis for a Single Model Configuration

3.4.1. Types of Simulation Outputs

• Terminating Simulation: Fixed time horizon or event count (e.g., simulate one day, one project). Initial conditions matter.

• Steady-State (Non-Terminating) Simulation: Long-run behavior (e.g., factory throughput over years). Initial conditions should not affect results after warm-up.

3.4.2. Statistical Analysis Challenges

• Correlation/Autocorrelation: Successive observations are dependent (especially in steady-state).

• Non-Stationarity: Early observations may reflect initialization bias (warm-up period).

3.4.3. Analysis of Terminating Runs

Method of Independent Replications:

1. Perform $n$ independent replications (different RNG streams).

2. For each replication $j$, compute estimator $$\displaystyle \hat{\theta}_j $$ (e.g., mean throughput).

3. Point estimator: $$\displaystyle \bar{\theta} = \frac{1}{n} \sum_{j=1}^{n} \hat{\theta}_j $$

4. Sample variance across replications: $$\displaystyle S^2 = \frac{1}{n-1} \sum_{j=1}^{n} (\hat{\theta}_j - \bar{\theta})^2 $$

5. Confidence Interval (CI) for true mean $\mu$:

$$\bar{\theta} \pm t_{\alpha/2, n-1} \frac{S}{\sqrt{n}}$$

where $t$ is t-distribution critical value.

Determining $n$:

• Start with $n \approx 10-15$.

• Check half-width of CI relative to $\bar{\theta}$ (desired precision, e.g., 5% relative error).

• Increase $n$ if CI too wide.

3.4.4. Analysis of Steady-State Runs

Warm-up Period (Initial Transient):

• Purpose: Discard initial biased observations.

• Methods:

  • Welch's Method: Plot moving averages of cumulative means; discard until plot stabilizes.

  • Relative Precision: Run long simulation, compute CI for mean over increasing time; stop when CI width stabilizes.

  • Autocorrelation: Check autocorrelation function; discard until autocorrelation becomes negligible.

Method of Batch Means:

1. After warm-up, run simulation for total time $T$.

2. Divide remaining observations into $k$ large, non-overlapping batches of equal size $m$ (so $$\displaystyle T = k \cdot m $$).

3. Compute batch means: $$\displaystyle \bar{Y}_i = \frac{1}{m} \sum_{t \in \text{batch } i} Y_t $$

4. Treat $$\displaystyle \bar{Y}_1, ..., \bar{Y}_k $$ as approximately independent (if $m$ large enough to exceed autocorrelation).

5. CI for steady-state mean $\mu$:

$$\bar{Y} \pm t_{\alpha/2, k-1} \frac{S_{\bar{Y}}}{\sqrt{k}}$$

where $$\displaystyle \bar{Y} = \frac{1}{k} \sum \bar{Y}_i $$, $$\displaystyle S_{\bar{Y}}^2 = \frac{1}{k-1} \sum (\bar{Y}_i - \bar{Y})^2 $$.

Replication-Deletion Approach:

• Perform $n$ replications, each with warm-up deletion.

• Treat each replication's post-warm-up mean as independent observation.

• Compute CI as in terminating case.

3.4.5. Comparing Scenarios (System Configurations)

Paired-t Comparison (for correlated outputs):

• Use Common Random Numbers (CRN): same RNG streams across scenarios.

• For replication $j$, compute difference $$\displaystyle D_j = \hat{\theta}_{j,A} - \hat{\theta}_{j,B} $$.

• CI for difference $$\displaystyle \mu_A - \mu_B $$:

$$\bar{D} \pm t_{\alpha/2, n-1} \frac{S_D}{\sqrt{n}}$$

where $$\displaystyle \bar{D} = \frac{1}{n} \sum D_j $$, $$\displaystyle S_D^2 = \frac{1}{n-1} \sum (D_j - \bar{D})^2 $$.

• Most powerful when CRN induces positive correlation.

Independent-t Comparison (for independent replications):

• Scenario A: $$\displaystyle n_A $$ replications, mean $$\displaystyle \bar{\theta}_A $$, variance $$\displaystyle S_A^2 $$.

• Scenario B: $$\displaystyle n_B $$ replications, mean $$\displaystyle \bar{\theta}_B $$, variance $$\displaystyle S_B^2 $$.

• Pooled variance (if equal variances assumed):

$$S_p^2 = \frac{(n_A-1)S_A^2 + (n_B-1)S_B^2}{n_A + n_B - 2}$$

• CI:

$$(\bar{\theta}_A - \bar{\theta}_B) \pm t_{\alpha/2, n_A+n_B-2} \cdot S_p \sqrt{\frac{1}{n_A} + \frac{1}{n_B}}$$

Multiple Comparison Procedures:

• When comparing $$\displaystyle >2 $$ systems, use Bonferroni adjustment: test at $$\displaystyle \alpha' = \alpha / \text{number of comparisons} $$ to control family-wise error rate.

[!TIP] Critical Rule: Use paired-t when same RNG streams used (CRN); independent-t when different streams or independent replications.

3.5. Advanced Topics in Output Analysis

3.5.1. Variance Reduction Techniques (VRTs)

Goal: Reduce variance of estimator for same computational effort.

• Common Random Numbers (CRN):

  • Use identical RNG streams for all competing scenarios.

  • Induces positive correlation between corresponding outputs.

  • Reduces variance of difference estimator (paired-t).

  • Must ensure synchronization of random events across scenarios.

• Antithetic Variates:

  • For each run with $$\displaystyle U_1,...,U_n $$, also run with $$\displaystyle 1-U_1,...,1-U_n $$.

  • Creates negatively correlated pairs.

  • Estimator: $$\displaystyle \hat{\theta} = \frac{1}{2n} \sum_{i=1}^{n} [\phi(U_i) + \phi(1-U_i)] $$

  • Effective when $\phi(u)$ is monotonic.

• Control Variates:

  • Use output $Y$ from a correlated control variable $X$ with known mean $E[X]$.

  • Adjusted estimator: $$\displaystyle \hat{\theta}_{cv} = \bar{Y} + b (\bar{X} - E[X]) $$

  • Choose $$\displaystyle b = \frac{\text{Cov}(Y,X)}{\text{Var}(X)} $$ (optimal).

  • Requires $E[X]$ known and $X$ correlated with $Y$.

• Importance Sampling:

  • Change underlying probability measure to sample more "important" regions (e.g., rare events).

  • Weight observations by likelihood ratio $$\displaystyle w(x) = \frac{f(x)}{g(x)} $$, where $f$=original, $g$=new distribution.

  • Estimator: $$\displaystyle \hat{\theta} = \frac{1}{n} \sum w(x_i) h(x_i) $$.

3.5.2. Optimization via Simulation

Challenges:

• Stochastic objective function (noisy).

• No gradient information.

• Expensive function evaluations.

Methods Overview:

• Ordinary Search: Bisection, Fibonacci search (for 1-D).

• Response Surface Methodology (RSM):

  • Fit polynomial (usually quadratic) to response from designed experiments.

  • Find optimum via gradient/ascent on fitted surface.

• Stochastic Approximation (e.g., Robbins-Monro):

  • Iterative: $$\displaystyle \theta_{n+1} = \theta_n + a_n (Y_n - \text{target}) $$, where $$\displaystyle a_n $$ step size.

  • Estimates gradient from noisy observations.

• Heuristics:

  • Simulated Annealing: Probabilistic search allowing uphill moves.

  • Genetic Algorithms: Population-based evolution (selection, crossover, mutation).

[!TIP] Exam Focus: Know when to use paired-t vs independent-t. Understand CRN as both VRT and comparison tool.