Simulation & Modelling

State variables describe the current condition of the system at any point in time:

• Number of customers in the queue — how many are waiting
• Server status — busy (1) or idle (0)
• Number of customers in the system — queue + those being served
• Time of last event — when the last arrival or departure occurred

Events in a Single-Server Model

• Arrival event — a new customer joins the system
• Departure/Service completion event — server finishes with a customer; next customer (if any) begins service
• End-of-simulation event — the stopping condition is reached

• Large solution space — there may be millions of possible combinations of input variables to explore
• Noisy responses — simulation output is stochastic; the same input can give different outputs in different runs
• Computational cost — each simulation run takes time; evaluating all options is infeasible
• Local vs global optima — algorithms may get stuck in a local optimum and miss the global best
• Warm-up period — initial transient effects bias results if not handled correctly

A queuing simulation tracks customers arriving, waiting, being served, and leaving. The key metrics to compute are:

Given: inter-arrival times and service times for each customer.

1. Arrival Time = Previous arrival time + Inter-arrival time (Customer 1 arrives at time 0)
2. Time Service Begins = MAX(Arrival Time, Previous Service End Time)
3. Wait in Queue = Time Service Begins − Arrival Time
4. Time Service Ends = Time Service Begins + Service Time
5. Time in System = Time Service Ends − Arrival Time
6. Idle Time of Server = MAX(0, Arrival Time − Previous Service End Time) [only if server was free when customer arrived]

From the 20-customer simulation table given in the exam paper, here's how to compute the three required metrics:

Step 1: Extract the totals from the table

Step 2: Calculate the metrics

When arrival rate λ and service rate μ are given (not a simulation table), use these formulae:

Let ρ = λ/μ (traffic intensity / server utilisation)

Given IAT and ST sequences, build the simulation table the same way as the bank example:

IAT: 0 1 1 6 3 7 5 2 4 1 → Customer 1 arrives at t=0, Customer 2 at t=0+1=1, etc.

Once the table is built, compute:

• Average waiting time = Sum of all queue wait times / 8 customers
• Idle time of server = Sum of all server idle time entries
• Average service time = Sum of all service times / 8 customers

Monte Carlo simulation uses random numbers to sample from probability distributions and simulate the behaviour of a system many times, then averages the results.

Steps:

1. Identify the uncertain variables and their probability distributions
2. Assign random number ranges to each outcome (cumulative distribution)
3. Use random numbers to look up simulated values
4. Calculate the outcome for each trial
5. Average the results across all trials

Convert a probability table into cumulative probabilities, then map ranges of random numbers to each value:

Train capacity = 51 trucks. Shortfall = max(0, 51 − production). Excess = max(0, production − 51).

Random numbers: 05, 34, 78, 56, 45, 90, 04, 58, 92, 39

Avg shortfall: 17/10 = 1.7 trucks/day  |  Avg excess: 6/10 = 0.6 trucks/day

Formula: Profit = (Selling Price − Variable Cost) × Demand − Investment

Each trial: use random numbers to look up a selling price, a variable cost, and a demand. Calculate profit for that trial. Average across all 10 trials.

• Uniformly distributed on [0,1]
• Statistically independent — knowing one number tells you nothing about the next
• Long period — the sequence should not repeat for a very long time
• Reproducible — same seed gives same sequence (useful for debugging)
• Fast to generate computationally

Where: X₀ = seed, a = multiplier, c = increment, m = modulus

To get U ∈ [0,1]: divide by m: U = X / m

When c = 0, it's the Multiplicative Congruential Generator:

X₁ = (19 × 63) mod 100 = 1197 mod 100 = 97
X₂ = (19 × 97) mod 100 = 1843 mod 100 = 43
X₃ = (19 × 43) mod 100 = 817 mod 100  = 17
X₄ = (19 × 17) mod 100 = 323 mod 100  = 23
X₅ = (19 × 23) mod 100 = 437 mod 100  = 37

Skip-Ahead Formula

This lets you jump n steps ahead without computing all intermediate values.

Combines three multiplicative generators for a very long period (~8×10¹²):

Parameters: a₁=157, m₁=32363  |  a₂=146, m₂=31727  |  a₃=142, m₃=31657

X₁,₁ = (157 × 100) mod 32363 = 15700 mod 32363 = 15700
X₂,₁ = (146 × 300) mod 31727 = 43800 mod 31727 = 12073
X₃,₁ = (142 × 500) mod 31657 = 71000 mod 31657 = 7686

Y₁ = (15700 - 12073 + 7686) mod 32362 = 11313 mod 32362 = 11313
U₁ = 11313 / 32363 ≈ 0.3496

Proposed by John von Neumann. To generate n-digit numbers:

1. Start with an n-digit seed
2. Square it to get a 2n-digit number
3. Extract the middle n digits as the next number
4. Repeat

Fibonacci Generator

# R implementation
fibonacci_rng <- function(n, x0) {
  x <- numeric(n + 2)
  x[1] <- x0[1]; x[2] <- x0[2]  # need 2 seeds
  for (i in 3:(n+2)) x[i] <- (x[i-1] + x[i-2]) %% 1
  return(x[3:(n+2)])
}

Mitchell-Moore Generator

# R implementation
mitchell_moore <- function(n, x0) {  # x0 must be length 55
  x <- c(x0, numeric(n))
  for (i in 56:(55+n)) x[i] <- (x[i-24] + x[i-55]) %% 1
  return(x[56:(55+n)])
}

To transform U ~ Uniform[0,1] to Uniform[a, b]:

Examples:

• Uniform[−11, 17]: X = −11 + (17 − (−11)) × U = −11 + 28U
• Uniform[−148, 1128]: X = −148 + 1276U
• Uniform[−160, 200]: X = −160 + 360U

The most common technique. If F(x) is the CDF, set F(x) = U (uniform random number) and solve for x:

Works perfectly when F⁻¹ has a closed form.

Exponential Distribution (mean = β)

Generate 5 values using U: 0.12, 0.45, 0.78, 0.33, 0.91 (example)
X = -3.7 × ln(U)
X₁ = -3.7 × ln(0.12) = -3.7 × (-2.12) = 7.85 days
X₂ = -3.7 × ln(0.45) = -3.7 × (-0.80) = 2.95 days  ...

Uniform Distribution on [a, b]

Normal Distribution (mean=μ, variance=σ²)

Use Box-Muller or standard normal table. In R: rnorm(n, mean=μ, sd=σ)

Weibull Distribution (shape α, scale β)

For the fluorescent bulb: P(T > t) = e^(−(t/1000)⁵). Set U = e^(−(T/1000)⁵), solve for T:

Discrete Distribution (Empirical)

Build a CDF table. For each U, find the x value where F(x) first exceeds U.

Given: F(x) = x(x+1)(2x+1) / [n(n+1)(2n+1)] for x = 1, 2, …, n. With n=4:

Find x such that F(x−1) < U ≤ F(x):

• R₁ = 0.83 → F(3)=0.467 < 0.83 ≤ F(4)=1.000 → x = 4
• R₂ = 0.24 → F(2)=0.167 < 0.24 ≤ F(3)=0.467 → x = 3
• R₃ = 0.57 → F(3)=0.467 < 0.57 ≤ F(4)=1.000 → x = 4

Range (a, b), mode at c. True mean = (a + b + c) / 3

For range (1,10), mode at 4: True mean = (1+10+4)/3 = 5

Beta Distribution

In R: rbeta(n, shape1=β₁, shape2=β₂) generates values on [0,1]. To transform to [L, H]: X = L + (H−L) × U_beta

# Exponential (mean = 3.7)
rexp(5, rate = 1/3.7)

# Normal (mean=33, variance=4, so sd=2)
rnorm(5, mean=33, sd=2)

# Gamma (shape=3, scale=1)
rgamma(10000, shape=3, scale=1)

# Beta
rbeta(10, shape1=1.47, shape2=2.16)

# Weibull (shape=5, scale=1000)
rweibull(n, shape=5, scale=1000)

# Chi-squared (k degrees of freedom)
rchisq(n, df=k)  # or: sum of k squared standard normals

# Binomial P(X = 25000) where X ~ Bin(50000, 0.5)
# Method 1:
dbinom(25000, size=50000, prob=0.5)

# Method 2 (using sum, log, exp only):
n <- 50000; k <- 25000; p <- 0.5
log_prob <- sum(log(1:n)) - sum(log(1:k)) - sum(log(1:(n-k))) + k*log(p) + (n-k)*log(1-p)
exp(log_prob)

The simulation table tracks: Clock time, Server status, Queue length, Future events.

Given: IAT = 1, 5, 6, 3, 8 and ST = 3, 5, 4, 1, 5

Two probability tables: arrivals per day and unloading rate per day. Build RN ranges for each, then simulate day by day tracking:

• Barges arriving each day
• Barges unloaded each day (limited by unloading rate and barges available)
• Barges delayed (queue carrying over to next day)

Metrics: Avg barges delayed = Total delayed / Days; Avg arrivals/day; Avg offloaded/day.

When arrival rates change over time (e.g., rush hours), the arrival rate λ(t) is estimated from data:

Example: 8:00–10:00 has counts 22, 24, 20, 28 across 4 days → avg = 94/4 = 23.5 arrivals → λ = 23.5/2 = 11.75 arrivals/hour

When data shows autocorrelation (successive observations are related), use time-series models.

AR(1) — Autoregressive Model

Where φ is the autocorrelation coefficient and ε_n ~ Normal(0, σ²(1−φ²)). Works well for continuous, normally distributed data.

EAR(1) — Exponential Autoregressive Model

Better for count data (integers, non-negative). Uses a mixture approach with a Bernoulli random variable.

To choose between them: Look at a histogram of the data. If bell-shaped/symmetric → AR(1). If skewed/count-like → EAR(1).

Hotel patrons data (20, 14, 21, 19, 14, 18, 21, 25, 27, 26, 22, 18, 13, 18, 18, 18, 25, 23, 20, 21) appears roughly normal → AR(1) likely fits better.

# R code for AR(1)
x <- c(20,14,21,19,14,18,21,25,27,26,22,18,13,18,18,18,25,23,20,21)
fit <- arima(x, order=c(1,0,0))
print(fit)
hist(x, main="Hotel Patrons")

Test if simulation output is consistent with real system behaviour.

Data: 7 replications: 18.9, 22.0, 19.4, 22.1, 19.8, 21.9, 20.2. Real-world average: μ₀ = 22.5

H₀: μ = 22.5 (model is valid) vs H₁: μ ≠ 22.5

n = 7
x̄ = (18.9+22.0+19.4+22.1+19.8+21.9+20.2)/7 = 144.3/7 = 20.614
s² = Σ(xᵢ − x̄)² / (n−1)

Deviations: -1.714, 1.386, -1.214, 1.486, -0.814, 1.286, -0.414
Squared:     2.939,  1.921,  1.474,  2.209,  0.663,  1.654,  0.171
Sum = 11.031  →  s² = 11.031/6 = 1.839  →  s = 1.356

t = (x̄ − μ₀) / (s/√n) = (20.614 − 22.5) / (1.356/√7)
  = −1.886 / 0.5124 = −3.68

Critical value: t₀.₀₂₅,₆ = 2.447 (two-tailed, α=0.05)

|t| = 3.68 > 2.447 → Reject H₀: the simulation output is NOT consistent with the real system. The model needs revision.

Data: 0.740, 1.28, 1.46, 2.36, 0.354, 0.750, 0.912, 4.44, 0.114, 3.08, 3.24, 1.10, 1.59, 1.47, 1.17, 1.27, 9.12, 11.5, 2.42, 1.77

1. Compute summary statistics: mean, variance, min, max
2. Plot histogram to observe shape (skewed right → try Exponential or Gamma)
3. Estimate parameters: For Exponential, mean ≈ 1/λ so λ̂ = 1/x̄
4. Goodness-of-fit test: Chi-square or K-S test to confirm fit

# R code for input modelling
data <- c(0.740,1.28,1.46,2.36,0.354,0.750,0.912,4.44,0.114,
           3.08,3.24,1.10,1.59,1.47,1.17,1.27,9.12,11.5,2.42,1.77)
mean(data); var(data); sd(data)
hist(data, probability=TRUE)
# If variance ≈ mean², try exponential
# If right-skewed, try gamma or lognormal
library(MASS)
fitdistr(data, "exponential")
fitdistr(data, "gamma")

Where Oᵢ = observed frequency, Eᵢ = expected frequency under the hypothesised distribution.

Degrees of freedom = k − 1 − (number of estimated parameters)

• For Poisson with unknown mean: df = k − 2 (estimated λ from data)
• For Poisson with specified mean: df = k − 1

Reject H₀ if χ² > χ²_{α, df} from table.

Given: IAT = 0,1,1,6,3,7,5,2 and ST = 4,2,5,4,1,5,4,1 for 8 customers (Customer 1 at t=0)

Total wait: 0+3+4+3+4+0+0+2 = 16 | Avg wait: 16/8 = 2.0 min

Total idle: 2 min | Sim ends at: t=28

P(wait): 5/8 = 0.625 | P(idle): 2/28 ≈ 0.071

Avg svc time: (4+2+5+4+1+5+4+1)/8 = 26/8 = 3.25 min

Formula: X = −11 + 28U

Given U ~ Uniform[0,1], generate X ~ Uniform[−11, 17]:

• U = 0.25 → X = −11 + 28(0.25) = −11 + 7 = −4
• U = 0.80 → X = −11 + 28(0.80) = −11 + 22.4 = 11.4

Set U = F(x) = x⁴/16, solve for x:

True mean: E[X] = ∫₀² x × f(x)dx where f(x) = x³/4

Generate sample values using U₁, U₂, ... → Xᵢ = 2Uᵢ^(0.25). Sample mean should ≈ 1.6.

simulate_die <- function() {
  rolls <- 0
  while (TRUE) {
    rolls <- rolls + 1
    if (sample(1:6, 1) == 6) return(rolls)
  }
}

# Run it
simulate_die()

# Many runs to get average
mean(replicate(10000, simulate_die()))  # Should be ≈ 6

Data: injuries per month over 100 months. Test if Poisson distributed.

data <- c(rep(0,35), rep(1,40), rep(2,13), rep(3,6), rep(4,4), rep(5,1), rep(6,1))
lambda_hat <- mean(data)  # = (0×35+1×40+2×13+3×6+4×4+5×1+6×1)/100 = 100/100 = 1.0

# Expected frequencies (n=100, Poisson with λ=1)
E <- 100 * dpois(0:6, lambda=lambda_hat)
# Observed
O <- c(35, 40, 13, 6, 4, 1, 1)

# Merge last cells if E < 5
# chi-square
chi_sq <- sum((O - E)^2 / E)
df <- length(O) - 1 - 1  # minus 1 for estimated lambda
qchisq(0.95, df)  # critical value

Queuing Metrics

Random Number Generation

Random Variate Generation

Statistical Tests

1. Q1 (Compulsory): 2 benefits + 1 limitation of simulation; 4 real-world applications; steps in simulation; queuing table (compute avg wait, P_wait, P_idle)
2. Monte Carlo Production Simulation: Build RN ranges from probability table; simulate 10 days; compute shortfalls/excesses
3. Random Number Generation: LCG computation (X₀, a, c, m given); L'Ecuyer combined generator; transforming U[0,1] to U[a,b]
4. Random Variate Generation: CDF inversion; exponential; discrete CDF table
5. R Code: Simulation of distributions (gamma, chi-square, exponential, die roll)
6. AR(1)/EAR(1): Fitting to hotel patron data; deciding which fits better from histogram
7. Validation t-test: job shop replication data vs. μ₀ = 22.5