Discrete Event Simulation: Queues and Servers

;;------------------------------------------------------------------------------
;; Copyright 2017 Nick Bennett
;;
;; Licensed under the Apache License, Version 2.0 (the "License");
;; you may not use this file except in compliance with the License.
;; You may obtain a copy of the License at
;;
;;     http://www.apache.org/licenses/LICENSE-2.0
;;
;; Unless required by applicable law or agreed to in writing, software
;; distributed under the License is distributed on an "AS IS" BASIS,
;; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
;; See the License for the specific language governing permissions and
;; limitations under the License.
;;------------------------------------------------------------------------------

;;------------------------------------------------------------------------------
;; Version History
;; 2010-02-21 - v1.0:   M/M/n queueing simulation; up to 10 servers; standard
;;                      aggregate statistics & corresponding expected values.
;; 2017-07-26 - v1.5:   Updated for NetLogo 6. Licensed under Apache License
;;                      2.0.
;; 2017-07-29 - v1.6:   Added plot & more NetLogo 6-related updates, including
;;                      use of anonymous procedures and range reporter.
;; 2017-08-03 - v1.6.1: Fix info details.
;;------------------------------------------------------------------------------


;; The model has two types of agents: customers, who enter and wait in a queue
;; until a server is available; and servers, who serve each customer in
;; first-come, first-served order.

breed [customers customer]
breed [servers server]


;; Each customer records the time entering the system, and the time entering
;; service, so that average time-in-queue/time-in-system statistics can be
;; computed.

customers-own [
  time-entered-queue
  time-entered-service
]


;; Each server records the customer agent being served, and the scheduled
;; completion of that service. Since servers are homogenous, individual
;; utilization statistics aren't kept.

servers-own [
  customer-being-served
  next-completion-time
]


globals [
  ; Waiting line
  queue
  ; Arrival process
  arrival-count
  next-arrival-time
  ; Statistics for average load/usage of queue and servers
  stats-start-time
  total-customer-queue-time
  total-customer-service-time
  ; Statistics for average time-in-queue and time-in-service
  total-time-in-queue
  total-time-in-system
  total-queue-throughput
  total-system-throughput
  ; Physical layout parameters
  server-ycor
  queue-server-offset
  customer-xinterval
  customer-yinterval
  customers-per-row
  ; Saved slider values to allow detection of changes during a simulation run
  save-mean-arrival-rate
  save-mean-service-time
  ; Theoretical measures, computed analytically using classic queueing theory
  expected-utilization
  expected-queue-length
  expected-queue-time
  ; Anonymous procedures
end -run-task
  arrive-task
  complete-service-task
  reset-stats-task
]


;; Sets up the model when loaded.

to startup
  setup
end 


;; Initializes global variables and server agents.

to setup
  clear-all
  setup-globals
  setup-servers
  setup-tasks
  compute-theoretical-measures
  reset-ticks
  reset-stats
  schedule-arrival
end 


;; Resets statistics, initializes queue list, and sets agent shapes and other
;; display properties.

to setup-globals
  set queue []
  set next-arrival-time 0
  set arrival-count 0
  set-default-shape servers "server"
  set-default-shape customers "person"
  set server-ycor (min-pycor + 1)
  set queue-server-offset 1.5
  set customer-xinterval 0.5
  set customer-yinterval 1
  set customers-per-row (1 + (world-width - 1) / customer-xinterval)
  set save-mean-arrival-rate mean-arrival-rate
  set save-mean-service-time mean-service-time
end 


;; Creates server agents and arranges them horizontally, evenly spaced along
;; the bottom of the NetLogo world. This layout is purely cosmetic, and has no
;; functional purpose or impact.

to setup-servers
  let horizontal-interval (world-width / number-of-servers)
  create-servers number-of-servers [
    set color green
    setxy (min-pxcor - 0.5 + horizontal-interval * (0.5 + who)) server-ycor
    set size 2.75
    set label ""
    set customer-being-served nobody
    set next-completion-time 0
  ]
end 


;; Sets up anonymous procedures for event queue entries.

to setup-tasks
  set end-run-task [[?ignore] -> end-run]
  set arrive-task [[?ignore] -> arrive]
  set complete-service-task [[?server] -> complete-service ?server]
  set reset-stats-task [[?ignore] -> reset-stats]
end 


;; Updates statistics (which also advances the clock) and invokes the
;; relevant anonymous procedure for the next event scheduled.

to go
  ifelse (ticks < max-run-time) [
    let next-event []
    let event-queue (list (list max-run-time end-run-task nobody))
    let next-server-to-complete next-server-complete
    set event-queue (
      fput (list next-arrival-time arrive-task nobody) event-queue)
    if (is-turtle? next-server-to-complete) [
      set event-queue (fput
        (list
          ([next-completion-time] of next-server-to-complete)
          complete-service-task
          next-server-to-complete)
        event-queue)
    ]
    if (stats-reset-time > ticks) [
      set event-queue (
        fput (list stats-reset-time reset-stats-task nobody) event-queue)
    ]
    set event-queue (sort-by [[?1 ?2] -> first ?1 < first ?2] event-queue)
    set next-event (first event-queue)
    update-usage-stats (first next-event)
    set next-event (but-first next-event)
    (run (first next-event) (last next-event))
  ]
  [
    stop
  ]
end 


;; Ends the execution of the simulation. In fact, this procedure does nothing,
;; but is still necessary. When the associated event is the first in the event
;; queue, the clock will be updated to the simulation end time prior to this
;; procedure being invoked; this causes the go procedure to stop on the next
;; iteration.

to end-run
  ; Do nothing
end 


;; Creates a new customer agent, adds it to the queue, and attempts to start
;; service.

to arrive
  let color-index (arrival-count mod 70)
  let main-color (floor (color-index / 5))
  let shade-offset (color-index mod 5)
  create-customers 1 [
    set color (3 + shade-offset + main-color * 10)
    set time-entered-queue ticks
    move-forward length queue
    set queue (lput self queue)
    set time-entered-queue ticks
  ]
  set arrival-count (arrival-count + 1)
  schedule-arrival
  begin-service
end 


;; Samples from the exponential distribution to schedule the time of the next
;; customer arrival in the system.

to schedule-arrival
  set next-arrival-time (ticks + random-exponential (1 / mean-arrival-rate))
end 


;; If there are customers in the queue, and at least one server is idle, starts
;; service on the first customer in the queue, using a randomly selected
;; idle server, and generating a complete-service event with a time sampled
;; from the exponential distribution. Updates the queue display, moving each
;; customer forward.

to begin-service
  let available-servers (servers with [not is-agent? customer-being-served])
  if (not empty? queue and any? available-servers) [
    let next-customer (first queue)
    let next-server one-of available-servers
    set queue (but-first queue)
    ask next-customer [
      set time-entered-service ticks
      set total-time-in-queue
        (total-time-in-queue + time-entered-service - time-entered-queue)
      set total-queue-throughput (total-queue-throughput + 1)
      move-to next-server
    ]
    ask next-server [
      set customer-being-served next-customer
      set next-completion-time (ticks + random-exponential mean-service-time)
      set label precision next-completion-time 3
      set color red
    ]
    (foreach queue (range (length queue)) [
      [?1 ?2] ->
      ask ?1 [
        move-forward ?2
      ]
    ])
  ]
end 


;; Updates time-in-system statistics, removes current customer agent, returns
;; the server to the idle state, and attempts to start service on another
;; customer.

to complete-service [?server]
  ask ?server [
    set total-time-in-system (total-time-in-system + ticks
      - [time-entered-queue] of customer-being-served)
    set total-system-throughput (total-system-throughput + 1)
    ask customer-being-served [
      die
    ]
    set customer-being-served nobody
    set next-completion-time 0
    set color green
    set label ""
  ]
  begin-service
end 


;; Reports the busy server with the earliest scheduled completion.

to-report next-server-complete
  report (min-one-of
    (servers with [is-agent? customer-being-served]) [next-completion-time])
end 


;; Sets all aggregate statistics back to 0 - except for the simulation start
;; time (used for computing average queue length and average server
;; utilization), which is set to the current time (which is generally not 0,
;; for a reset-stats event).

to reset-stats
  set total-customer-queue-time 0
  set total-customer-service-time 0
  set total-time-in-queue 0
  set total-time-in-system 0
  set total-queue-throughput 0
  set total-system-throughput 0
  set stats-start-time ticks
end 


;; Updates the usage/utilization statistics and advances the clock to the
;; specified event time.

to update-usage-stats [event-time]
  let delta-time (event-time - ticks)
  let busy-servers (servers with [is-agent? customer-being-served])
  let in-queue (length queue)
  let in-process (count busy-servers)
  let in-system (in-queue + in-process)
  set total-customer-queue-time
    (total-customer-queue-time + delta-time * in-queue)
  set total-customer-service-time
    (total-customer-service-time + delta-time * in-process)
  tick-advance (event-time - ticks)
  update-plots
end 


;; Moves to the specified queue position, based on the global spacing
;; parameters. This queue display is purely cosmetic, and has no functional
;;  purpose or impact.

to move-forward [queue-position]
  let new-xcor
    (max-pxcor - customer-xinterval * (queue-position mod customers-per-row))
  let new-ycor (server-ycor + queue-server-offset
    + customer-yinterval * floor (queue-position / customers-per-row))
  ifelse (new-ycor > max-pycor) [
    hide-turtle
  ]
  [
    setxy new-xcor new-ycor
    if (hidden?) [
      show-turtle
    ]
  ]
end 


;; Checks to see if the values of the mean-arrival-rate and mean-service-time
;; sliders have changed since the last time that the theoretical system measures
;; were calculated; if so, the theoretical measures are recalculated.
;; (Note: This reporter is invoked by the expected queue length monitor. Since
;; monitors are updated periodically, even when a forever button isn't pressed,
;; this allows the theoretical statistics to be updated when the sliders
;; change, even if the model isn't running.

to-report sliders-changed?
  let changed? false
  if ((save-mean-arrival-rate != mean-arrival-rate)
      or (save-mean-service-time != mean-service-time)) [
    set changed? true
    set save-mean-arrival-rate mean-arrival-rate
    set save-mean-service-time mean-service-time
    compute-theoretical-measures
  ]
  report changed?
end 


;; Computes the expected utilization, queue length, and time in queue for M/M/n
;; queueing system.

to compute-theoretical-measures
  let balance-factor (mean-arrival-rate * mean-service-time)
  let n (count servers)
  ifelse ((balance-factor / n) < 1) [
    let k 0
    let k-sum 1
    let power-product 1
    let factorial-product 1
    let busy-probability 0
    foreach (range 1 n) [
      [?] ->
      set power-product (power-product * balance-factor)
      set factorial-product (factorial-product * ?)
      set k-sum (k-sum + power-product / factorial-product)
    ]
    set power-product (power-product * balance-factor)
    set factorial-product (factorial-product * n)
    set k (k-sum / (k-sum + power-product / factorial-product))
    set busy-probability ((1 - k) / (1 - balance-factor * k / n))
    set expected-utilization (balance-factor / n)
    set expected-queue-length
      (busy-probability * expected-utilization / (1 - expected-utilization))
    set expected-queue-time
      (busy-probability * mean-service-time / (n * (1 - expected-utilization)))
  ]
  [
    set expected-utilization 1
    set expected-queue-length "N/A"
    set expected-queue-time "N/A"
  ]
end