# Design and Analysis of 5G Scenarios

#### Iñaki Ucar, José Alberto Hernández, Pablo Serrano, Arturo Azcorra

#### 2020-06-06

## Abstract

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

## Crosshauling of FH and BH traffic

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:

    C <- 40e9                                   # link capacity [bps]
    rho <- 0.75                                 # utilization

    FH_EX <- (80 + 46) * 8 / C                  # avg. FH interarrival time [s]
    BH_bytes <- c(40, 576, 1500)                # BH traffic distribution
    Weights <- c(7, 4, 1) / 12                  #
    BH_EX <- sum(Weights * (8 * BH_bytes / C))  # avg. BH interarrival time [s]

    FH_w <- 0.5                                 # FH traffic ratio
    BH_w <- 1 - FH_w                            # BH traffic ratio
    FH_lambda <- FH_w * rho / FH_EX             # FH avg. rate [pkts/s]
    BH_lambda <- BH_w * rho / BH_EX             # BH avg. rate [pkts/s]

    Tsim <- 1e5 / (FH_lambda + BH_lambda)       # simulation time

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
    simulate <- function(param) {

      # single FH traffic trajectory traversing n_switch elements sequentially
      fh_traffic <- lapply(1:param[["n_switch"]], function(i) {
        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
      bh_traffic <- lapply(1:param[["n_switch"]], function(i) {
        trajectory() %>%
          seize(paste0("switch", i)) %>%
          timeout(function() sample(BH_bytes*8/C, size=1, prob=Weights)) %>%
          release(paste0("switch", i))
      })

      # simulation environment
      env <- simmer() %>%
        # 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
    cases <- expand.grid(n_switch = c(1, 2, 5),
                         fh_prio = c(0, 1),
                         preemptive = c(TRUE, FALSE))

    # parallel simulation
    system.time({
      envs <- parallel::mclapply(split(cases, 1:nrow(cases)), simulate,
                                 mc.cores=nrow(cases), mc.set.seed=FALSE)
    })

Finally, the information is extracted, summarised and represented in a few lines of code:

    library(tidyverse)

    bp.vals <- function(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95)) {
      r <- quantile(x, probs=probs, na.rm=TRUE)
      names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
      r
    }

    arrivals <- envs %>%
      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")

## Mobile traffic backhauling with FTTx

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:

- Deployment of a Small Cell, reducing the amount and requirements of the generated traffic.
- Deployment of an RRH, following the C-RAN paradigm discussed above, which would therefore generate time-sensitive FH traffic.

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:

    C <- 1.25e9                           # link capacity [bps]
    Tg <- 1e-6                            # guard time [s]

    bytes <- c(40, 576, 1500)             # traffic distribution
    Weights <- c(7, 4, 1) / 12            #
    EX <- sum(Weights * (8 * bytes / C))  # avg. interarrival time [s]
    burst <- 20                           # burst length [pkts]

    ONU_lambda <- 20e6 / C / EX           # avg. rate per ONU [pkts/s]
    SC_lambda <- 150e6 / C / EX           # avg. rate for the SC [pkts/s]
    RRH_on <- 8 * 6000 / C                # RRH on period [s]
    RRH_period <- RRH_on * C / 720e6      # RRH total period [s]
    RRH_off <- RRH_period - RRH_on        # RRH off period [s]

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
    cyclic_counter <- function(n) {
      n <- as.numeric(n)
      i <- 1
      function(op) {
        if (!missing(op)) {
          i <<- i + op
          if (i > n) i <<- 1
          else if (i < 1) i <<- n
        }
        i
      }
    }

    # helper function: generator of packet sizes of n packets
    gen_pkts <- function(n) sample(bytes*8/C, size=n, prob=Weights, replace=TRUE)

    # helper function: double-exponential generator of interarrival times
    gen_exp <- function(lambda, burst) function()
      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
    simulate <- function(param) {

      # global variables
      lambda <- rep(ONU_lambda, param[["ONU_n"]])
      limit <- rep(param[["limit"]], param[["ONU_n"]]) * 8 / C
      if (param[["scenario"]] == "A") {
        param[["ONU_n"]] <- param[["ONU_n"]] + 1
        lambda <- c(lambda, SC_lambda)
        limit <- c(limit, Inf)
        tx_t <- RRH_on + Tg + Tg * 0:param[["ONU_n"]]
      } else {
        Tg_n <- RRH_off %/% Tg - 1
        Tg_i <- rep(1:Tg_n, length.out=param[["ONU_n"]]+1)
        period_i <- rep(0:param[["ONU_n"]], each=Tg_n, length.out=param[["ONU_n"]]+1)
        tx_t <- RRH_on + RRH_period * period_i + Tg * Tg_i
      }
      eq_pkts <- rep(list(NULL), param[["ONU_n"]])
      tx_pkts <- rep(list(NULL), param[["ONU_n"]])
      t_head <- tail(tx_t, 1)
      onu_i <- cyclic_counter(param[["ONU_n"]])
      remaining <- Inf

      # DBA logic
      set_next_window <- function() {
        # time until the next RRH_on
        if (param[["scenario"]] == "B")
          remaining <- (t_head %/% RRH_period + 1) * RRH_period - t_head
        # generate new pkt lengths
        pkts <- get_queue_count(env, paste0("ONU", onu_i()))
        eq_pkts[[onu_i()]] <<- c(eq_pkts[[onu_i()]], gen_pkts(pkts - length(eq_pkts[[onu_i()]])))
        # reserve the next transmission window
        eq_cumsum <- cumsum(eq_pkts[[onu_i()]])
        n <- sum(eq_cumsum <= min(limit[[onu_i()]], remaining - Tg))
        if ((pkts && !n) || (remaining <= Tg)) {
          # go past the next RRH_on
          t_head <<- t_head + remaining + RRH_on + Tg
          n <- sum(eq_cumsum <= min(limit[[onu_i()]], RRH_off - 2*Tg))
        }
        tx_pkts[[onu_i()]] <<- head(eq_pkts[[onu_i()]], n)
        eq_pkts[[onu_i()]] <<- tail(eq_pkts[[onu_i()]], -n)
        tx_t[[onu_i()]] <<- t_head
        t_head <<- t_head + sum(tx_pkts[[onu_i()]]) + Tg
        index <<- 0
        tx_t[[onu_i(1)]] - now(env)
      }

      # list of ONU_n trajectories, one per ONU
      ONUs <- lapply(1:param[["ONU_n"]], function(i) {
        trajectory() %>%
          seize(paste0("ONU", i)) %>%
          set_capacity(paste0("ONU", i), 0) %>%
          seize("link") %>%
          timeout(function() {
            index <<- index + 1
            tx_pkts[[i]][index]
          }) %>%
          release("link") %>%
          release(paste0("ONU", i))
      })

      # OLT logic
      OLT <- trajectory() %>%
        simmer::select(function() paste0("ONU", onu_i())) %>%
        set_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
      RRH <- trajectory() %>%
        seize("link") %>%
        timeout(RRH_on) %>%
        release("link") %>%
        timeout(RRH_off) %>%
        rollback(amount=4, times=Inf)

      # simulation environment
      env <- simmer() %>%
        # the fiber link as a resource
        add_resource("link", 1, Inf)

      if (param[["scenario"]] == "B")
        # RRH worker
        env %>% add_generator("RRH_", RRH, at(0))

      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
    cases <- data.frame(scenario = c(rep("A", 4), rep("B", 1)),
                        ONU_n = c(rep(31, 4), rep(7, 1)),
                        limit = c(Inf, 1500, 3000, 6000, Inf))

    # parallel simulation
    system.time({
      envs <- parallel::mclapply(split(cases, 1:nrow(cases)), simulate,
                                 mc.cores=nrow(cases), mc.set.seed=FALSE)
    })

Finally, the information is extracted, summarised and represented in a few lines of code:

    library(tidyverse)

    bp.vals <- function(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95)) {
      r <- quantile(x, probs=probs, na.rm=TRUE)
      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(
        scenario, `SmallCell + 31 ONUs`="A", `RRH + 7 ONUs`="B")) %>%
      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]")

## Energy efficiency for massive IoT

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:

    m <- 9            # number of RA retries
    Ps <- 0.03        # sleep power consumption [mW]
    Pi <- 10          # idle power consumption [mW]
    Prx <- 100        # rx power consumption [mW]
    Prbp <- 32.18     # tx power consumption per RB pair [mW]
    Ppre <- 32.18     # preamble tx power consumption [mW]

    Tbpre <- 0.0025   # time before preamble transmission [s]
    Tpre <- 0.0014    # time for preamble transmission [s]
    Trarx <- 0.01     # rx time to perform RA [s]
    Tcrcp <- 0.016    # sum of processing delays to establish a connection [s]
    Twait <- 0.054    # S1 processing and transfer delay [s]

    Brbp <- 36        # bytes per RB pair (QPSK modulation)
    Breq <- 7         # RRC request message size (bytes)
    Bcompcp <- 129    # RRC setup complete + NAS control plane SR + data for CP (bytes)

    Tra <- Tbpre + Trarx + Tpre
    Pra <- (Tbpre*Pi + Trarx*Prx + Tpre*Ppre) / Tra
    rbs <- (ceiling(Breq/Brbp) + ceiling(Bcompcp/Brbp)) * 1e-3
    Ttx <- Tcrcp + rbs + Twait
    Ptx <- (Tcrcp*Prx + rbs*Prbp + Twait*Prx) / Ttx

    tx_period <- 3600 # time between transmissions (seconds)
    Tsim <- 24*3600   # simulation time (seconds)

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
    simulate <- function(param) {

      # identifiers for RA preambles
      preambles <- paste0("preamble_", 1:54)

      # IoT device logic
      meter <- trajectory() %>%
        trap("reading") %>%
        # sleep
        set_attribute("P", 0) %>%
        wait() %>%
        timeout(function() round(runif(1, 0, param[["backoff"]]), 3)) %>%
        # ra start
        simmer::select(preambles, policy="random") %>%
        seize_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
      trigger <- trajectory() %>%
        timeout(tx_period) %>%
        send("reading") %>%
        rollback(2, times=Inf)

      # simulation environment
      env <- simmer() %>%
        # 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
        env %>% add_resource(i, 1, 0, mon=FALSE)

      env %>%
        run(until=Tsim) %>%
        wrap()
    }

    # grid of scenarios
    cases <- expand.grid(meter_n = c(5e3, 1e4, 3e4), backoff = c(5, 10, 30, 60))

    # parallel simulation
    system.time({
      envs <- parallel::mclapply(split(cases, 1:nrow(cases)), simulate,
                                 mc.cores=nrow(cases), mc.set.seed=FALSE)
    })

Finally, the information is extracted, summarised and represented in a few lines of code:

    library(tidyverse)

    bp.vals <- function(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95)) {
      r <- quantile(x, probs=probs, na.rm=TRUE)
      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]")