# R
simulate_hpp <- function(lambda, T_max) {
t <- 0
times <- c()
while (TRUE) {
t <- t + rexp(1, lambda)
if (t > T_max) break
times <- c(times, t)
}
times
}31 Workshop Activities
In this workshop, we will practice simulating
- Homogeneous Poisson processes (HPP)
- Non-homogeneous Poisson processes (NHPP)
- M/M/\(c\) queues
Exercise 1: Simulating HPP
Consider a homogeneous Poisson process with rate \(\lambda = 5\) events per hour.
- Simulate the process over a time interval of 8 hours (i.e., from \(t=0\) to \(t=8\)) using the inter-arrival time method. Set seed using set.seed(1234) before simulating the process.
- Plot the number of events \(N(t)\) against time \(t\) for the simulated process.
- Calculate the average number of events per hour from your simulation and compare it to the theoretical rate \(\lambda\).
- Repeat the simulation 1000 times and plot the distribution of the total number of events in 8 hours. Does it match the expected Poisson distribution with mean \(\lambda \times T\)?
- Calculate the empirical probability of observing exactly 40 events in 8 hours and compare it to the theoretical probability from the Poisson distribution
- Calculate the empirical probability of observing more than 50 events in 8 hours and compare it to the theoretical probability from the Poisson distribution.
Hint: HPP can be simulated by generating inter-arrival times from an exponential distribution. You may use the following function to simulate a HPP:
# Python
import numpy as np
def simulate_hpp(lmbda, T_max):
t = 0.0
times = []
while True:
t += np.random.exponential(scale=1 / lmbda)
if t > T_max:
break
times.append(t)
return np.array(times)Exercise 2: Simulating a NHPP with a piecewise intensity function
A non-homogeneous Poisson process has an intensity function that varies over time. Consider the following piecewise intensity function:
\[ \lambda_t = \begin{cases} 1 & \text{if } 0 \leq t < 3 \\ 4 & \text{if } 3 \leq t < 7 \\ 2 & \text{if } 7 \leq t \leq 10 \end{cases} \]
- Plot the intensity function \(\lambda(t)\) over the time interval \([0, 10]\).
- Use the thinning algorithm to simulate a NHPP with this intensity function over the time interval \([0, 10]\). Set seed using
set.seed(1234)before simulating the process. - Report the number of events that occurred in the simulated NHPP.
- Plot the cumulative count of events over time, i.e., plot \(N(t)\) against \(t\).
- Do you think the thinning algorithm is efficient for simulating this NHPP? Why or why not?
Hint: NHPP can be simulated using the thinning algorithm, which involves simulating a HPP with a rate equal to the maximum intensity of the NHPP and then accepting or rejecting events based on the ratio of the NHPP intensity to the maximum intensity. You may use the following function to simulate a NHPP using thinning:
# R
simulate_nhpp_thinning <- function(lambda_fun, lambda_max, T_max) {
t <- 0
times <- c()
while (TRUE) {
t <- t + rexp(1, lambda_max)
if (t > T_max) break
if (runif(1) <= lambda_fun(t) / lambda_max) {
times <- c(times, t)
}
}
times
}# Python
import numpy as np
def simulate_nhpp_thinning(lambda_fun, lambda_max, T_max):
t = 0.0
times = []
while True:
# Propose next event time from HPP(lambda_max)
t += np.random.exponential(scale=1 / lambda_max)
if t > T_max:
break
# Accept–reject step
if np.random.rand() <= lambda_fun(t) / lambda_max:
times.append(t)
return np.array(times)Exercise 3: Simulating M/M/\(c\) (Extended)
Consider the M/M/2 queueing system of the EV charging station dicussed in 🚗 EV Charging Bays Simulation with the following parameters:
- initial state (number of customers) \(X_0 = 0\).
- arrival rate \(\lambda = 3\)
- service rate \(\mu = 2\) per server
- a total simulation time of 3000 hours (250 working days for 12 hours per day).
- Conduct a comparative analysis if EV adoption increases, i.e., re-run the simulation with \(\lambda = 3.5\) cars/hour. How do the average waiting time (\(W_q\)) and utilisation (\(\rho\)) change?
- What if charging is slower in winter? Increase the mean charging time to 37.5 minutes (i.e., \(\mu = 1.6\)). How does this affect utilisation and waiting?
- The university is planning to add more chargers to the station as the analysis (from part (a) and (b)) shows that the current system will be heavily utilised if the arrival rate increases or if the service rate decreases only by just a little. Suppose that in the next couple years the arrival rate is expected to increase to 4 cars/hour and the mean charging time during winter is 37.5 minutes (i.e., \(\mu = 1.6\)), how many chargers should the university plan to install?
Hint: M/M/c Simulation Function in R
simulate_mmc <- function(lambda = 2, mu = 3, c = 2,
sim_time = 10, seed = 1234) {
if (!is.null(seed)) set.seed(seed)
# System state
T <- 0
X <- 0
# Output dataframe
df <- data.frame(
time = 0,
state = 0,
event = "Start",
server = NA,
start = NA,
end = NA
)
# Next arrival
T_a <- rexp(1, lambda)
# Active services: one row per busy server
services <- data.frame(
server = integer(0),
start = numeric(0),
end = numeric(0)
)
while (T < sim_time) {
T_s <- if (nrow(services) > 0) min(services$end) else Inf
if (min(T_a, T_s) > sim_time) break
if (T_a < T_s) {
# Arrival
T <- T_a
X <- X + 1
df <- rbind(
df,
data.frame(
time = T, state = X, event = "Arrival",
server = NA, start = NA, end = NA
)
)
T_a <- T + rexp(1, lambda)
# Start service immediately if a server is free
if (X <= c) {
busy <- services$server
idle <- setdiff(seq_len(c), busy)
srv <- sample(idle, 1) # assign to a random free server
dur <- rexp(1, mu)
services <- rbind(
services,
data.frame(
server = srv,
start = T,
end = T + dur
)
)
}
} else {
# Departure
idx <- which.min(services$end)
srv <- services$server[idx]
s0 <- services$start[idx]
s1 <- services$end[idx]
T <- s1
X <- X - 1
df <- rbind(
df,
data.frame(
time = T, state = X, event = "Departure",
server = srv, start = s0, end = s1
)
)
# Remove completed service
services <- services[-idx, ]
# If queue non‑empty, same server starts next service immediately
if (X >= c) {
dur <- rexp(1, mu)
services <- rbind(
services,
data.frame(
server = srv,
start = T,
end = T + dur
)
)
}
}
}
rownames(df) <- NULL
df
}Hint: M/M/c Simulation Function in Python
import numpy as np
import pandas as pd
def simulate_mmc(lambda_=2, mu=3, c=2, sim_time=10, seed=1234):
if seed is not None:
np.random.seed(seed)
# System state
T = 0.0 # current time
X = 0 # number in system
# Output data
records = [{
"time": 0.0,
"state": 0,
"event": "Start",
"server": np.nan,
"start": np.nan,
"end": np.nan
}]
# Next arrival
T_a = np.random.exponential(scale=1 / lambda_)
# Active services: one row per busy server
services = pd.DataFrame(columns=["server", "start", "end"])
while T < sim_time:
T_s = services["end"].min() if not services.empty else np.inf
if min(T_a, T_s) > sim_time:
break
if T_a < T_s:
# ----------------
# Arrival event
# ----------------
T = T_a
X += 1
records.append({
"time": T,
"state": X,
"event": "Arrival",
"server": np.nan,