Simulation frameworks are important tools for the analysis and design of communication networks and protocols, but they can result extremely costly and/or complex (for the case of very specialized tools), or too naive and lacking proper features and support (for the case of ad-hoc tools). In this paper, we present an analysis of three 5G scenarios using
simmer
, a recent R package for discrete-event simulation that sits between the above two paradigms. As our results show, it provides a simple yet very powerful syntax, supporting the efficient simulation of relatively complex scenarios at a low implementation cost.
This vignette contains the code associated to the article Design and Analysis of 5G Scenarios with simmer
: An R Package for Fast DES Prototyping (see the draft version on arXiv), published in the IEEE Communications Magazine (see citation("simmer")
). Refer to the article for a full description and analysis of each scenario.
This scenario is motivated by the Cloud Radio Access Network (C-RAN) paradigm, where the mobile base station functionality is split into simple Remote Radio Heads (RRH), spread across the deployment and connected by fiber to centralized (and possibly virtualized) Base-Band Units (BBU), at the operators’ premises.
In this C-RAN paradigm, fronthaul (FH) traffic from the RRH has stringent delay requirements, while backhaul (BH) traffic from the BBU has mild delay requirements. In a general topology, such as the one illustrated in the figure, packet switches will forward both types of traffic. We use simmer
to simulate the scenario and decide whether introducing service differentiation might improve the ability to fulfil the delivery guarantees of FH traffic.
These are the configuration parameters considered:
<- 40e9 # link capacity [bps]
C <- 0.75 # utilization
rho
<- (80 + 46) * 8 / C # avg. FH interarrival time [s]
FH_EX <- c(40, 576, 1500) # BH traffic distribution
BH_bytes <- c(7, 4, 1) / 12 #
Weights <- sum(Weights * (8 * BH_bytes / C)) # avg. BH interarrival time [s]
BH_EX
<- 0.5 # FH traffic ratio
FH_w <- 1 - FH_w # BH traffic ratio
BH_w <- FH_w * rho / FH_EX # FH avg. rate [pkts/s]
FH_lambda <- BH_w * rho / BH_EX # BH avg. rate [pkts/s]
BH_lambda
<- 1e5 / (FH_lambda + BH_lambda) # simulation time Tsim
The model is defined and encapsulated into the simulate
function. Then, several scenarios with different parameters are defined in cases
, which can be run in parallel:
library(simmer)
set.seed(1234)
# parametrized simulation function
# - param[["n_switch"]] is the number of switches
# - param[["fh_prio"]] is the priority level for FH traffic
# - param[["preemptive"]] is the preemptiveness flag for switches
<- function(param) {
simulate
# single FH traffic trajectory traversing n_switch elements sequentially
<- lapply(1:param[["n_switch"]], function(i) {
fh_traffic trajectory() %>%
seize(paste0("switch", i)) %>%
timeout(FH_EX) %>%
release(paste0("switch", i))
%>% join()
})
# list of n_switch trajectories, one per element, modeling the interfering BH traffic
<- lapply(1:param[["n_switch"]], function(i) {
bh_traffic trajectory() %>%
seize(paste0("switch", i)) %>%
timeout(function() sample(BH_bytes*8/C, size=1, prob=Weights)) %>%
release(paste0("switch", i))
})
# simulation environment
<- simmer() %>%
env # generator of FH traffic
add_generator("FH_0_", fh_traffic, function() rexp(100, FH_lambda), priority=param[["fh_prio"]])
for (i in 1:param[["n_switch"]])
%>%
env # n_switch resources, one per switch
add_resource(paste0("switch", i), 1, Inf, mon=FALSE, preemptive=param[["preemptive"]]) %>%
# n_switch generators of BH traffic
add_generator(paste0("BH_", i, "_"), bh_traffic[[i]], function() rexp(100, BH_lambda))
%>%
env run(until=Tsim) %>%
wrap()
}
# grid of scenarios
<- expand.grid(n_switch = c(1, 2, 5),
cases fh_prio = c(0, 1),
preemptive = c(TRUE, FALSE))
# parallel simulation
system.time({
<- parallel::mclapply(split(cases, 1:nrow(cases)), simulate,
envs mc.cores=nrow(cases), mc.set.seed=FALSE)
})
Finally, the information is extracted, summarised and represented in a few lines of code:
library(tidyverse)
<- function(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95)) {
bp.vals <- quantile(x, probs=probs, na.rm=TRUE)
r names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
r
}
<- envs %>%
arrivals get_mon_arrivals() %>%
left_join(rowid_to_column(cases, "replication")) %>%
filter(!(fh_prio==0 & preemptive==TRUE)) %>%
separate(name, c("Flow", "index", "n")) %>%
mutate(total_time = end_time - start_time,
waiting_time = total_time - activity_time,
Queue = forcats::fct_recode(interaction(fh_prio, preemptive),
"without SP" = "0.FALSE",
"with SP" = "1.FALSE",
"with SP & preemption" = "1.TRUE"))
%>%
arrivals filter(n_switch == 1) %>%
# plot
ggplot(aes(Queue, waiting_time*1e9, color=Flow)) + theme_bw() +
stat_summary(fun.data=bp.vals, geom="boxplot", position="dodge") +
theme(legend.justification=c(0, 1), legend.position=c(.02, .98)) +
labs(y = "Queueing Delay [ns]", x = element_blank())
%>%
arrivals filter(Flow == "FH") %>%
# plot
ggplot(aes(factor(n_switch), waiting_time*1e9, color=Queue)) + theme_bw() +
stat_summary(fun.data=bp.vals, geom="boxplot", position="dodge") +
theme(legend.justification=c(0, 1), legend.position=c(.02, .98)) +
labs(y = "Queueing Delay [ns]", x = "No. of switches")
We next consider the case of a residencial area with a Fiber-To-The-Premises (FTTx) infrastructure, that is, an Optical Distribution Network (ODN), composed of the Optical Line Terminal (OLT), splitters, and the Optical Network Unit (ONU) at the users’ premises. As the figure illustrates, we assume that an operator is planning to deploy an antenna, carrying the mobile traffic over the ODN, and is considering two implementation options:
In both cases, we analyze the upstream channel of a Time-Division Multiplexed Passive Optical Network (TDM-PON) providing broadband access to the residential users and the mobile users.
These are the configuration parameters considered:
<- 1.25e9 # link capacity [bps]
C <- 1e-6 # guard time [s]
Tg
<- c(40, 576, 1500) # traffic distribution
bytes <- c(7, 4, 1) / 12 #
Weights <- sum(Weights * (8 * bytes / C)) # avg. interarrival time [s]
EX <- 20 # burst length [pkts]
burst
<- 20e6 / C / EX # avg. rate per ONU [pkts/s]
ONU_lambda <- 150e6 / C / EX # avg. rate for the SC [pkts/s]
SC_lambda <- 8 * 6000 / C # RRH on period [s]
RRH_on <- RRH_on * C / 720e6 # RRH total period [s]
RRH_period <- RRH_period - RRH_on # RRH off period [s] RRH_off
The model is defined and encapsulated into the simulate
function. Then, several scenarios with different parameters are defined in cases
, which can be run in parallel:
# helper function: round-robin-based selection of ONUs
<- function(n) {
cyclic_counter <- as.numeric(n)
n <- 1
i function(op) {
if (!missing(op)) {
<<- i + op
i if (i > n) i <<- 1
else if (i < 1) i <<- n
}
i
}
}
# helper function: generator of packet sizes of n packets
<- function(n) sample(bytes*8/C, size=n, prob=Weights, replace=TRUE)
gen_pkts
# helper function: double-exponential generator of interarrival times
<- function(lambda, burst) function()
gen_exp unlist(lapply(rexp(100/burst, lambda/burst), function(i) c(i, rep(0, rpois(1, burst)))))
library(simmer)
set.seed(1234)
# parametrized simulation function
# - param[["scenario"]] is the scenario identifier (A=small cell, B=RRH)
# - param[["ONU_n"]] is the number of ONUs
# - param[["limit"]] is the max. transmission window
<- function(param) {
simulate
# global variables
<- rep(ONU_lambda, param[["ONU_n"]])
lambda <- rep(param[["limit"]], param[["ONU_n"]]) * 8 / C
limit if (param[["scenario"]] == "A") {
"ONU_n"]] <- param[["ONU_n"]] + 1
param[[<- c(lambda, SC_lambda)
lambda <- c(limit, Inf)
limit <- RRH_on + Tg + Tg * 0:param[["ONU_n"]]
tx_t else {
} <- RRH_off %/% Tg - 1
Tg_n <- rep(1:Tg_n, length.out=param[["ONU_n"]]+1)
Tg_i <- rep(0:param[["ONU_n"]], each=Tg_n, length.out=param[["ONU_n"]]+1)
period_i <- RRH_on + RRH_period * period_i + Tg * Tg_i
tx_t
}<- rep(list(NULL), param[["ONU_n"]])
eq_pkts <- rep(list(NULL), param[["ONU_n"]])
tx_pkts <- tail(tx_t, 1)
t_head <- cyclic_counter(param[["ONU_n"]])
onu_i <- Inf
remaining
# DBA logic
<- function() {
set_next_window # time until the next RRH_on
if (param[["scenario"]] == "B")
<- (t_head %/% RRH_period + 1) * RRH_period - t_head
remaining # generate new pkt lengths
<- get_queue_count(env, paste0("ONU", onu_i()))
pkts onu_i()]] <<- c(eq_pkts[[onu_i()]], gen_pkts(pkts - length(eq_pkts[[onu_i()]])))
eq_pkts[[# reserve the next transmission window
<- cumsum(eq_pkts[[onu_i()]])
eq_cumsum <- sum(eq_cumsum <= min(limit[[onu_i()]], remaining - Tg))
n if ((pkts && !n) || (remaining <= Tg)) {
# go past the next RRH_on
<<- t_head + remaining + RRH_on + Tg
t_head <- sum(eq_cumsum <= min(limit[[onu_i()]], RRH_off - 2*Tg))
n
}onu_i()]] <<- head(eq_pkts[[onu_i()]], n)
tx_pkts[[onu_i()]] <<- tail(eq_pkts[[onu_i()]], -n)
eq_pkts[[onu_i()]] <<- t_head
tx_t[[<<- t_head + sum(tx_pkts[[onu_i()]]) + Tg
t_head <<- 0
index onu_i(1)]] - now(env)
tx_t[[
}
# list of ONU_n trajectories, one per ONU
<- lapply(1:param[["ONU_n"]], function(i) {
ONUs trajectory() %>%
seize(paste0("ONU", i)) %>%
set_capacity(paste0("ONU", i), 0) %>%
seize("link") %>%
timeout(function() {
<<- index + 1
index
tx_pkts[[i]][index]%>%
}) release("link") %>%
release(paste0("ONU", i))
})
# OLT logic
<- trajectory() %>%
OLT ::select(function() paste0("ONU", onu_i())) %>%
simmerset_attribute("ONU", function() onu_i()) %>%
set_capacity_selected(function() length(tx_pkts[[onu_i()]])) %>%
timeout(function() sum(tx_pkts[[onu_i()]])) %>%
timeout(set_next_window) %>%
rollback(amount=6, times=Inf)
# RRH logic
<- trajectory() %>%
RRH seize("link") %>%
timeout(RRH_on) %>%
release("link") %>%
timeout(RRH_off) %>%
rollback(amount=4, times=Inf)
# simulation environment
<- simmer() %>%
env # the fiber link as a resource
add_resource("link", 1, Inf)
if (param[["scenario"]] == "B")
# RRH worker
%>% add_generator("RRH_", RRH, at(0))
env
for (i in 1:param[["ONU_n"]])
%>%
env # ONU_n resources, one per ONU
add_resource(paste0("ONU", i), 0, Inf) %>%
# ONU_n traffic generators, one per ONU
add_generator(paste0("ONU_", i, "_"), ONUs[[i]], gen_exp(lambda[i], burst))
%>%
env # OLT worker
add_generator("token_", OLT, at(RRH_on + Tg), mon=2) %>%
run(until=1e5/sum(lambda)) %>%
wrap()
}
# grid of scenarios
<- data.frame(scenario = c(rep("A", 4), rep("B", 1)),
cases ONU_n = c(rep(31, 4), rep(7, 1)),
limit = c(Inf, 1500, 3000, 6000, Inf))
# parallel simulation
system.time({
<- parallel::mclapply(split(cases, 1:nrow(cases)), simulate,
envs mc.cores=nrow(cases), mc.set.seed=FALSE)
})
Finally, the information is extracted, summarised and represented in a few lines of code:
library(tidyverse)
<- function(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95)) {
bp.vals <- quantile(x, probs=probs, na.rm=TRUE)
r names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
r
}
%>%
envs get_mon_arrivals() %>%
left_join(rowid_to_column(cases, "replication")) %>%
mutate(scenario = forcats::fct_recode(
`SmallCell + 31 ONUs`="A", `RRH + 7 ONUs`="B")) %>%
scenario, separate(name, c("flow", "index", "n")) %>%
mutate(flow = if_else(index == ONU_n+1, "SmallCell", flow)) %>%
mutate(total_time = end_time - start_time,
waiting_time = total_time - activity_time) %>%
# plot
ggplot(aes(forcats::fct_rev(flow), waiting_time*1e6, color=factor(limit))) + theme_bw() +
facet_grid(~scenario, scales="free", space="free") +
stat_summary(fun.data=bp.vals, geom="boxplot", position="dodge") +
theme(legend.justification=c(0, 1), legend.position=c(.02, .98)) +
labs(y=expression(paste("Queueing Delay [", mu, "s]")), x=element_blank(), color="Limit [bytes]")
Finally, we consider the case of a massive Internet-of-Things (mIoT) scenario, a use case for Long Term Evolution (LTE) and next-generation 5G networks, as defined by the Third Generation Partnership Project (3GPP). As the figure (left) illustrates, we consider a single LTE macrocell in a dense urban area. The buildings in the cell area are populated with \(N\) smart meters (for electricity, gas and water), and each meter operates independently as a Narrowband IoT (NB-IoT) device. The devices’ behaviour is modeled following the diagram depicted in the figure (right).
These are the configuration parameters considered:
<- 9 # number of RA retries
m <- 0.03 # sleep power consumption [mW]
Ps <- 10 # idle power consumption [mW]
Pi <- 100 # rx power consumption [mW]
Prx <- 32.18 # tx power consumption per RB pair [mW]
Prbp <- 32.18 # preamble tx power consumption [mW]
Ppre
<- 0.0025 # time before preamble transmission [s]
Tbpre <- 0.0014 # time for preamble transmission [s]
Tpre <- 0.01 # rx time to perform RA [s]
Trarx <- 0.016 # sum of processing delays to establish a connection [s]
Tcrcp <- 0.054 # S1 processing and transfer delay [s]
Twait
<- 36 # bytes per RB pair (QPSK modulation)
Brbp <- 7 # RRC request message size (bytes)
Breq <- 129 # RRC setup complete + NAS control plane SR + data for CP (bytes)
Bcompcp
<- Tbpre + Trarx + Tpre
Tra <- (Tbpre*Pi + Trarx*Prx + Tpre*Ppre) / Tra
Pra <- (ceiling(Breq/Brbp) + ceiling(Bcompcp/Brbp)) * 1e-3
rbs <- Tcrcp + rbs + Twait
Ttx <- (Tcrcp*Prx + rbs*Prbp + Twait*Prx) / Ttx
Ptx
<- 3600 # time between transmissions (seconds)
tx_period <- 24*3600 # simulation time (seconds) Tsim
The model is defined and encapsulated into the simulate
function. Then, several scenarios with different parameters are defined in cases
, which can be run in parallel:
library(simmer)
set.seed(1234)
# parametrized simulation function
# - param[["meter_n"]] is the number of IoT devices in the cell
# - param[["backoff"]] is the backoff time
<- function(param) {
simulate
# identifiers for RA preambles
<- paste0("preamble_", 1:54)
preambles
# IoT device logic
<- trajectory() %>%
meter trap("reading") %>%
# sleep
set_attribute("P", 0) %>%
wait() %>%
timeout(function() round(runif(1, 0, param[["backoff"]]), 3)) %>%
# ra start
::select(preambles, policy="random") %>%
simmerseize_selected(
continue=c(TRUE, TRUE),
# ra & tx
post.seize=trajectory() %>%
set_attribute("P", Pra) %>%
timeout(Tra) %>%
release_selected() %>%
set_attribute("P", Ptx) %>%
timeout(Ttx),
# ra & backoff & retry
reject=trajectory() %>%
set_attribute("P", Pra) %>%
timeout(Tra) %>%
set_attribute("P", Pi) %>%
timeout(function() sample(1:20, 1) * 1e-3) %>%
rollback(6, times=m)
%>%
) rollback(5, times=Inf)
# trigger a reading for all the meters every tx_period
<- trajectory() %>%
trigger timeout(tx_period) %>%
send("reading") %>%
rollback(2, times=Inf)
# simulation environment
<- simmer() %>%
env # IoT device workers
add_generator("meter_", meter, at(rep(0, param[["meter_n"]])), mon=2) %>%
# trigger worker
add_generator("trigger_", trigger, at(0), mon=0)
for (i in preambles)
# one resource per preamble
%>% add_resource(i, 1, 0, mon=FALSE)
env
%>%
env run(until=Tsim) %>%
wrap()
}
# grid of scenarios
<- expand.grid(meter_n = c(5e3, 1e4, 3e4), backoff = c(5, 10, 30, 60))
cases
# parallel simulation
system.time({
<- parallel::mclapply(split(cases, 1:nrow(cases)), simulate,
envs mc.cores=nrow(cases), mc.set.seed=FALSE)
})
Finally, the information is extracted, summarised and represented in a few lines of code:
library(tidyverse)
<- function(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95)) {
bp.vals <- quantile(x, probs=probs, na.rm=TRUE)
r names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
r
}
%>%
envs get_mon_attributes() %>%
group_by(replication, name) %>%
summarise(dE = sum(c(0, head(value, -1) * diff(time)))) %>%
left_join(rowid_to_column(cases, "replication")) %>%
# plot
ggplot(aes(factor(meter_n), dE*tx_period/Tsim, color=factor(backoff))) +
stat_summary(fun.data=bp.vals, geom="boxplot", position="dodge") +
theme_bw() + theme(legend.justification=c(0, 1), legend.position=c(.02, .98)) +
labs(y="Energy per Tx [mJ]", x="No. of devices", color="Backoff window [s]")