This notebook simulates a taxi stand problem posed OR Stack Exchange. There is a single queue with Markovian arrivals and multiple taxis available. Taxis all have the same capacity and are only dispatched when full. Round trip travel times for taxis are constant. The measure of interest is average time between dispatches.

Parameters

These are the model parameters we will use.

nTaxis <- 5             # number of taxis
travelTime <- 20        # round trip travel time
loadSize <- 4           # number of passengers in a load
arrivalRate <- 0.7      # mean arrival rate of passengers

These are the simulation parameters.

seed <- 789              # random number seed
passengerLimit <- 10000  # number of passenger arrivals to simulate

Events

The simulation has three types of events, coded by the numbers 1, 2 and 3 in the simulation.

  1. Passenger arrival. This results in the queue being increased by one. If the queue length is loadSize and the taxi line contains at least one taxi, a dispatch event occurs immediately. All passenger arrivals are scheduled at the outset.

  2. Taxi dispatch. This results in loadSize passengers being removed from the queue and one taxi being removed from the line. A return event is scheduled for travelTime units in the future.

  3. Taxi return. This results in the taxi line count being increased by 1. If at least loadSize passengers are already in the queue, a dispatch event occurs immediately.

Event calendar

The event calendar is a data frame containing a time index, the type of event (1, 2, 3) scheduled at that time, the passenger queue size, and the size of the taxi line. We maintain a pointer to the row containing the next event to process.

We will need a function to insert an event into the calendar as well as a function to update the calendar when an event occurs.

# Inputs: t = time of the event, e = event type (1 - 3)
# Output; none (the calendar is expanded by one row)
schedule <- function(t, e) {
  # Find the indices of rows before/after the event.
  calendar <<- rbind(calendar[calendar$Time <= t, ],
                     c(t, e, 0, 0),
                     calendar[calendar$Time > t, ])
}

# Inputs: r = row index of an event in the calendar
# Output: none (the designated calendar row is updated, and a subsequent event may be scheduled)
process <- function(r) {
  if (calendar[r, "Event"] == 1) {
    # Passenger arrival. Copy the line length and bump the queue size.
    if (r == 1) {
      # Start of the simulation.
      calendar[r, "Queue"] <<- 1
    } else {
      calendar[r, "Line"] <<- calendar[r - 1, "Line"]
      calendar[r, "Queue"] <<- calendar[r - 1, "Queue"] + 1
      # If a taxi is available and can make load, schedule a dispatch.
      if ((calendar[r, "Queue"] == loadSize) & (calendar[r, "Line"] > 0)) {
        schedule(calendar[r, "Time"], 2)
      }
    }
  } else if (calendar[r, "Event"] == 2) {
    # Cab dispatch. Reduce the passenger queue and line contents and schedule a return.
    calendar[r, "Queue"] <<- calendar[r - 1, "Queue"] - loadSize
    calendar[r, "Line"] <<- calendar[r - 1, "Line"] - 1
    schedule(calendar[r, "Time"] + travelTime, 3)
  } else if (calendar[r, "Event"] == 3) {
    # Cab return. Bump the line by 1 and copy the queue size. If the queue is at least the load size, schedule an immediate dispatch.
    calendar[r, "Queue"] <<- calendar[r - 1, "Queue"]
    calendar[r, "Line"] <<- calendar[r - 1, "Line"] + 1
    if (calendar[r, "Queue"] >= loadSize) {
      schedule(calendar[r, "Time"], 2)
    }
  } else {
    # This should never happen.
    stop("Invalid event code: ", calendar[r, "Event"])
  }
}

Initialization

We initialize the event calendar by scheduling all passenger arrivals.

# Seed the random number generator.
set.seed(seed)
# Generate the interarrival times.
iat <- rexp(passengerLimit, arrivalRate)
# Set up the event calendar.
calendar <- data.frame(Time = cumsum(iat), Event = 1, Queue = 0, Line = nTaxis)
nextEvent <- 1  # row with next event to process
rm(iat)

Running the simulation.

We process the event pointed to by nextEvent (and increment the pointer) until it exceeds the number of rows in the data frame.

while (nextEvent <= nrow(calendar)) {
  process(nextEvent)
  nextEvent <- nextEvent + 1
}

Analyzing interdispatch times

We want to focus on the time intervals between dispatches.

# Calculate time intervals between consecutive dispatches.
tbd <- calendar[calendar$Event == 2, "Time"] |> diff()
# Print the sample mean time between dispatches.
cat("Mean TBD = ", mean(tbd), ".\n")
Mean TBD =  5.702537 .
# Print the ratio of load size to arrival rate.
cat("Load size / arrival rate = ", loadSize / arrivalRate, ".\n")
Load size / arrival rate =  5.714286 .
LS0tCnRpdGxlOiAiVGF4aSBTaW11bGF0aW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIG5vdGVib29rIHNpbXVsYXRlcyBhIFt0YXhpIHN0YW5kIHByb2JsZW1dKGh0dHBzOi8vb3Iuc3RhY2tleGNoYW5nZS5jb20vcXVlc3Rpb25zLzc3NzEvYXZlcmFnZS10aW1lLWJldHdlZW4tdHdvLWRpc3BhdGNoZXMtaW4tYS10YXhpLWZsZWV0LXByb2JhYmx5LWEtYmF0Y2gtcHJvY2Vzc2luZykgcG9zZWQgT1IgU3RhY2sgRXhjaGFuZ2UuIFRoZXJlIGlzIGEgc2luZ2xlIHF1ZXVlIHdpdGggTWFya292aWFuIGFycml2YWxzIGFuZCBtdWx0aXBsZSB0YXhpcyBhdmFpbGFibGUuIFRheGlzIGFsbCBoYXZlIHRoZSBzYW1lIGNhcGFjaXR5IGFuZCBhcmUgb25seSBkaXNwYXRjaGVkIHdoZW4gZnVsbC4gUm91bmQgdHJpcCB0cmF2ZWwgdGltZXMgZm9yIHRheGlzIGFyZSBjb25zdGFudC4gVGhlIG1lYXN1cmUgb2YgaW50ZXJlc3QgaXMgYXZlcmFnZSB0aW1lIGJldHdlZW4gZGlzcGF0Y2hlcy4KCiMgUGFyYW1ldGVycwoKVGhlc2UgYXJlIHRoZSBtb2RlbCBwYXJhbWV0ZXJzIHdlIHdpbGwgdXNlLgoKYGBge3J9Cm5UYXhpcyA8LSA1ICAgICAgICAgICAgICMgbnVtYmVyIG9mIHRheGlzCnRyYXZlbFRpbWUgPC0gMjAgICAgICAgICMgcm91bmQgdHJpcCB0cmF2ZWwgdGltZQpsb2FkU2l6ZSA8LSA0ICAgICAgICAgICAjIG51bWJlciBvZiBwYXNzZW5nZXJzIGluIGEgbG9hZAphcnJpdmFsUmF0ZSA8LSAwLjcgICAgICAjIG1lYW4gYXJyaXZhbCByYXRlIG9mIHBhc3NlbmdlcnMKYGBgCgpUaGVzZSBhcmUgdGhlIHNpbXVsYXRpb24gcGFyYW1ldGVycy4KCmBgYHtyfQpzZWVkIDwtIDc4OSAgICAgICAgICAgICAgIyByYW5kb20gbnVtYmVyIHNlZWQKcGFzc2VuZ2VyTGltaXQgPC0gMTAwMDAgICMgbnVtYmVyIG9mIHBhc3NlbmdlciBhcnJpdmFscyB0byBzaW11bGF0ZQpgYGAKCiMgRXZlbnRzCgpUaGUgc2ltdWxhdGlvbiBoYXMgdGhyZWUgdHlwZXMgb2YgZXZlbnRzLCBjb2RlZCBieSB0aGUgbnVtYmVycyAxLCAyIGFuZCAzIGluIHRoZSBzaW11bGF0aW9uLgoKMS4gUGFzc2VuZ2VyIGFycml2YWwuIFRoaXMgcmVzdWx0cyBpbiB0aGUgcXVldWUgYmVpbmcgaW5jcmVhc2VkIGJ5IG9uZS4gSWYgdGhlIHF1ZXVlIGxlbmd0aCBpcyBgbG9hZFNpemVgIGFuZCB0aGUgdGF4aSBsaW5lIGNvbnRhaW5zIGF0IGxlYXN0IG9uZSB0YXhpLCBhIGRpc3BhdGNoIGV2ZW50IG9jY3VycyBpbW1lZGlhdGVseS4gQWxsIHBhc3NlbmdlciBhcnJpdmFscyBhcmUgc2NoZWR1bGVkIGF0IHRoZSBvdXRzZXQuCgoyLiBUYXhpIGRpc3BhdGNoLiBUaGlzIHJlc3VsdHMgaW4gYGxvYWRTaXplYCBwYXNzZW5nZXJzIGJlaW5nIHJlbW92ZWQgZnJvbSB0aGUgcXVldWUgYW5kIG9uZSB0YXhpIGJlaW5nIHJlbW92ZWQgZnJvbSB0aGUgbGluZS4gQSByZXR1cm4gZXZlbnQgaXMgc2NoZWR1bGVkIGZvciBgdHJhdmVsVGltZWAgdW5pdHMgaW4gdGhlIGZ1dHVyZS4KCjMuIFRheGkgcmV0dXJuLiBUaGlzIHJlc3VsdHMgaW4gdGhlIHRheGkgbGluZSBjb3VudCBiZWluZyBpbmNyZWFzZWQgYnkgMS4gSWYgYXQgbGVhc3QgYGxvYWRTaXplYCBwYXNzZW5nZXJzIGFyZSBhbHJlYWR5IGluIHRoZSBxdWV1ZSwgYSBkaXNwYXRjaCBldmVudCBvY2N1cnMgaW1tZWRpYXRlbHkuCgojIEV2ZW50IGNhbGVuZGFyCgpUaGUgZXZlbnQgY2FsZW5kYXIgaXMgYSBkYXRhIGZyYW1lIGNvbnRhaW5pbmcgYSB0aW1lIGluZGV4LCB0aGUgdHlwZSBvZiBldmVudCAoMSwgMiwgMykgc2NoZWR1bGVkIGF0IHRoYXQgdGltZSwgdGhlIHBhc3NlbmdlciBxdWV1ZSBzaXplLCBhbmQgdGhlIHNpemUgb2YgdGhlIHRheGkgbGluZS4gV2UgbWFpbnRhaW4gYSBwb2ludGVyIHRvIHRoZSByb3cgY29udGFpbmluZyB0aGUgbmV4dCBldmVudCB0byBwcm9jZXNzLgoKV2Ugd2lsbCBuZWVkIGEgZnVuY3Rpb24gdG8gaW5zZXJ0IGFuIGV2ZW50IGludG8gdGhlIGNhbGVuZGFyIGFzIHdlbGwgYXMgYSBmdW5jdGlvbiB0byB1cGRhdGUgdGhlIGNhbGVuZGFyIHdoZW4gYW4gZXZlbnQgb2NjdXJzLgoKYGBge3J9CiMgSW5wdXRzOiB0ID0gdGltZSBvZiB0aGUgZXZlbnQsIGUgPSBldmVudCB0eXBlICgxIC0gMykKIyBPdXRwdXQ7IG5vbmUgKHRoZSBjYWxlbmRhciBpcyBleHBhbmRlZCBieSBvbmUgcm93KQpzY2hlZHVsZSA8LSBmdW5jdGlvbih0LCBlKSB7CiAgIyBGaW5kIHRoZSBpbmRpY2VzIG9mIHJvd3MgYmVmb3JlL2FmdGVyIHRoZSBldmVudC4KICBjYWxlbmRhciA8PC0gcmJpbmQoY2FsZW5kYXJbY2FsZW5kYXIkVGltZSA8PSB0LCBdLAogICAgICAgICAgICAgICAgICAgICBjKHQsIGUsIDAsIDApLAogICAgICAgICAgICAgICAgICAgICBjYWxlbmRhcltjYWxlbmRhciRUaW1lID4gdCwgXSkKfQoKIyBJbnB1dHM6IHIgPSByb3cgaW5kZXggb2YgYW4gZXZlbnQgaW4gdGhlIGNhbGVuZGFyCiMgT3V0cHV0OiBub25lICh0aGUgZGVzaWduYXRlZCBjYWxlbmRhciByb3cgaXMgdXBkYXRlZCwgYW5kIGEgc3Vic2VxdWVudCBldmVudCBtYXkgYmUgc2NoZWR1bGVkKQpwcm9jZXNzIDwtIGZ1bmN0aW9uKHIpIHsKICBpZiAoY2FsZW5kYXJbciwgIkV2ZW50Il0gPT0gMSkgewogICAgIyBQYXNzZW5nZXIgYXJyaXZhbC4gQ29weSB0aGUgbGluZSBsZW5ndGggYW5kIGJ1bXAgdGhlIHF1ZXVlIHNpemUuCiAgICBpZiAociA9PSAxKSB7CiAgICAgICMgU3RhcnQgb2YgdGhlIHNpbXVsYXRpb24uCiAgICAgIGNhbGVuZGFyW3IsICJRdWV1ZSJdIDw8LSAxCiAgICB9IGVsc2UgewogICAgICBjYWxlbmRhcltyLCAiTGluZSJdIDw8LSBjYWxlbmRhcltyIC0gMSwgIkxpbmUiXQogICAgICBjYWxlbmRhcltyLCAiUXVldWUiXSA8PC0gY2FsZW5kYXJbciAtIDEsICJRdWV1ZSJdICsgMQogICAgICAjIElmIGEgdGF4aSBpcyBhdmFpbGFibGUgYW5kIGNhbiBtYWtlIGxvYWQsIHNjaGVkdWxlIGEgZGlzcGF0Y2guCiAgICAgIGlmICgoY2FsZW5kYXJbciwgIlF1ZXVlIl0gPT0gbG9hZFNpemUpICYgKGNhbGVuZGFyW3IsICJMaW5lIl0gPiAwKSkgewogICAgICAgIHNjaGVkdWxlKGNhbGVuZGFyW3IsICJUaW1lIl0sIDIpCiAgICAgIH0KICAgIH0KICB9IGVsc2UgaWYgKGNhbGVuZGFyW3IsICJFdmVudCJdID09IDIpIHsKICAgICMgQ2FiIGRpc3BhdGNoLiBSZWR1Y2UgdGhlIHBhc3NlbmdlciBxdWV1ZSBhbmQgbGluZSBjb250ZW50cyBhbmQgc2NoZWR1bGUgYSByZXR1cm4uCiAgICBjYWxlbmRhcltyLCAiUXVldWUiXSA8PC0gY2FsZW5kYXJbciAtIDEsICJRdWV1ZSJdIC0gbG9hZFNpemUKICAgIGNhbGVuZGFyW3IsICJMaW5lIl0gPDwtIGNhbGVuZGFyW3IgLSAxLCAiTGluZSJdIC0gMQogICAgc2NoZWR1bGUoY2FsZW5kYXJbciwgIlRpbWUiXSArIHRyYXZlbFRpbWUsIDMpCiAgfSBlbHNlIGlmIChjYWxlbmRhcltyLCAiRXZlbnQiXSA9PSAzKSB7CiAgICAjIENhYiByZXR1cm4uIEJ1bXAgdGhlIGxpbmUgYnkgMSBhbmQgY29weSB0aGUgcXVldWUgc2l6ZS4gSWYgdGhlIHF1ZXVlIGlzIGF0IGxlYXN0IHRoZSBsb2FkIHNpemUsIHNjaGVkdWxlIGFuIGltbWVkaWF0ZSBkaXNwYXRjaC4KICAgIGNhbGVuZGFyW3IsICJRdWV1ZSJdIDw8LSBjYWxlbmRhcltyIC0gMSwgIlF1ZXVlIl0KICAgIGNhbGVuZGFyW3IsICJMaW5lIl0gPDwtIGNhbGVuZGFyW3IgLSAxLCAiTGluZSJdICsgMQogICAgaWYgKGNhbGVuZGFyW3IsICJRdWV1ZSJdID49IGxvYWRTaXplKSB7CiAgICAgIHNjaGVkdWxlKGNhbGVuZGFyW3IsICJUaW1lIl0sIDIpCiAgICB9CiAgfSBlbHNlIHsKICAgICMgVGhpcyBzaG91bGQgbmV2ZXIgaGFwcGVuLgogICAgc3RvcCgiSW52YWxpZCBldmVudCBjb2RlOiAiLCBjYWxlbmRhcltyLCAiRXZlbnQiXSkKICB9Cn0KYGBgCgoKIyBJbml0aWFsaXphdGlvbgoKV2UgaW5pdGlhbGl6ZSB0aGUgZXZlbnQgY2FsZW5kYXIgYnkgc2NoZWR1bGluZyBhbGwgcGFzc2VuZ2VyIGFycml2YWxzLgoKYGBge3J9CiMgU2VlZCB0aGUgcmFuZG9tIG51bWJlciBnZW5lcmF0b3IuCnNldC5zZWVkKHNlZWQpCiMgR2VuZXJhdGUgdGhlIGludGVyYXJyaXZhbCB0aW1lcy4KaWF0IDwtIHJleHAocGFzc2VuZ2VyTGltaXQsIGFycml2YWxSYXRlKQojIFNldCB1cCB0aGUgZXZlbnQgY2FsZW5kYXIuCmNhbGVuZGFyIDwtIGRhdGEuZnJhbWUoVGltZSA9IGN1bXN1bShpYXQpLCBFdmVudCA9IDEsIFF1ZXVlID0gMCwgTGluZSA9IG5UYXhpcykKbmV4dEV2ZW50IDwtIDEgICMgcm93IHdpdGggbmV4dCBldmVudCB0byBwcm9jZXNzCnJtKGlhdCkKYGBgCgojIFJ1bm5pbmcgdGhlIHNpbXVsYXRpb24uCgpXZSBwcm9jZXNzIHRoZSBldmVudCBwb2ludGVkIHRvIGJ5IGBuZXh0RXZlbnRgIChhbmQgaW5jcmVtZW50IHRoZSBwb2ludGVyKSB1bnRpbCBpdCBleGNlZWRzIHRoZSBudW1iZXIgb2Ygcm93cyBpbiB0aGUgZGF0YSBmcmFtZS4KCmBgYHtyfQp3aGlsZSAobmV4dEV2ZW50IDw9IG5yb3coY2FsZW5kYXIpKSB7CiAgcHJvY2VzcyhuZXh0RXZlbnQpCiAgbmV4dEV2ZW50IDwtIG5leHRFdmVudCArIDEKfQpgYGAKCiMgQW5hbHl6aW5nIGludGVyZGlzcGF0Y2ggdGltZXMKCldlIHdhbnQgdG8gZm9jdXMgb24gdGhlIHRpbWUgaW50ZXJ2YWxzIGJldHdlZW4gZGlzcGF0Y2hlcy4KCmBgYHtyfQojIENhbGN1bGF0ZSB0aW1lIGludGVydmFscyBiZXR3ZWVuIGNvbnNlY3V0aXZlIGRpc3BhdGNoZXMuCnRiZCA8LSBjYWxlbmRhcltjYWxlbmRhciRFdmVudCA9PSAyLCAiVGltZSJdIHw+IGRpZmYoKQojIFByaW50IHRoZSBzYW1wbGUgbWVhbiB0aW1lIGJldHdlZW4gZGlzcGF0Y2hlcy4KY2F0KCJNZWFuIFRCRCA9ICIsIG1lYW4odGJkKSwgIi5cbiIpCiMgUHJpbnQgdGhlIHJhdGlvIG9mIGxvYWQgc2l6ZSB0byBhcnJpdmFsIHJhdGUuCmNhdCgiTG9hZCBzaXplIC8gYXJyaXZhbCByYXRlID0gIiwgbG9hZFNpemUgLyBhcnJpdmFsUmF0ZSwgIi5cbiIpCmBgYAoK