Discrete event with event queue

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

;-------------------------------------------------------------------------------
; v2.0 Copyright 2018, Russ Abbott, California State University, Los Angeles
;;                     (Russ.Abbott@gmail.com)
;;
;; 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.
;;
;; Contact Russ Abbott: Russ.Abbott@gmail.com
;-------------------------------------------------------------------------------

;-------------------------------------------------------------------------------
; Version History
; 2018-03-18 - v2.0  Uses and displays an actual event-queue.
;                    Updated to NetLogo 6.02.
;-------------------------------------------------------------------------------


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

; 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.
breed [customers customer]
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.
breed [servers server]
servers-own [
  customer-being-served
  next-completion-time
]


globals [
  ; Event-queue
  event-queue
  ; Waiting line
  customer-queue
  ; Arrival process
  arrival-count
  ; Statistics for average load/usage of customer-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
  customer-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
  ; max queue length
  max-queue
]

to startup
  setup
end 

; Initializes global variables and server agents.

to setup
  clear-all
  setup-globals
  setup-servers
  compute-theoretical-measures
  reset-ticks
  begin-stats-collection
  set-plot-pen-color magenta
  insert-event-into-event-queue (list begin-stats-collection-time (word "begin-stats-collection"))
  schedule-next-arrival
  print-event-queue
end 


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

to setup-globals
  set event-queue []
  set customer-queue []
  set arrival-count 0
  set-default-shape servers "server"
  set-default-shape customers "person"
  set server-ycor (min-pycor + 2)
  set customer-queue-server-offset 1.75
  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
  set max-queue 0
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 4
    set label ""
    set customer-being-served nobody
    set next-completion-time 0
  ]
end 


; The main loop.

to go
  update-theoretical-values
  ;; Event queue will eventually become empty because no events
  ;; are added with a scheduled time past max-run-time.
  if empty? event-queue [stop]
  let next-event first event-queue
  set event-queue but-first event-queue
  let event-time first next-event
  update-usage-stats event-time
  run second next-event
  print-event-queue
end 


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

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


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

to begin-service
  let available-servers (servers with [not is-agent? customer-being-served])
  if (not empty? customer-queue and any? available-servers) [
    let next-customer (first customer-queue)
    let next-server one-of available-servers
    set customer-queue (but-first customer-queue)
    ;; Move each customer in the queue to the position of the customer in front of it.
    (foreach customer-queue (range length customer-queue)
          [ [cust pos] -> ask cust [  move-to-queue-position pos ]  ] )
    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 precision (ticks + random-exponential mean-service-time) 2
      set label (word precision next-completion-time 2 " ")
      insert-event-into-event-queue (list next-completion-time (word "complete-service " who))
      set color red
    ]
  ]
end 


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

to begin-stats-collection
  set-plot-pen-color black
  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 

to-report begin-stats-collection-time
  report round (max-run-time * begin-stats-% / 100)
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-id]
  ask (server server-id) [
    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 

to-report event-to-string [event]
  report (word first event ": " second event)
end 


;; Don't insert events at or past max-run-time into the event queue.

to insert-event-into-event-queue [an-event]
  if first an-event <= max-run-time [
    set event-queue (insert-item-into-queue an-event event-queue)
  ]
end 

to-report insert-item-into-queue [an-event the-queue]
  if empty? the-queue [report (list an-event )]
  let first-event first the-queue
  report
    ifelse-value (first an-event <= first first-event)
           [fput an-event the-queue]
           [fput first-event insert-item-into-queue an-event but-first the-queue]
end 


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

to move-to-queue-position [queue-position]
  let new-xcor
    (max-pxcor - customer-xinterval * (queue-position mod customers-per-row))
  let new-ycor (server-ycor + customer-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 


; Reports the busy server with the earliest scheduled completion.

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

to print-event-queue
  clear-output
  if not empty? event-queue [
    let n (min (list 10 length event-queue))
    let events (take n event-queue)
    foreach but-last events [ s -> output-print event-to-string s ]
    output-type event-to-string last events
  ]
end 


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

to schedule-next-arrival
  let increment random-exponential (1 / mean-arrival-rate)
  insert-event-into-event-queue (list precision (ticks + increment) 2 "arrive")
end 

to-report second [a-list]
  report first but-first a-list
end 

;; the first n elements of a list

to-report take [n a-list]
  report ifelse-value (n <= 0 or empty? a-list)
             [ [] ]
             [ fput first a-list take (n - 1) but-first a-list ]
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 customer-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
  set max-queue max list max-queue length customer-queue
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Theoretical values calculations ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Compute the expected utilization, customer-queue length, and time in customer-queue for an M/M/n
; customer-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 (n-values (n - 1) [ ?1 -> ?1 + 1 ]) [ ?1 ->
      set power-product (power-product * balance-factor)
      set factorial-product (factorial-product * ?1)
      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 

; Checks to see if the values of the mean-arrival-rate and mean-service-time
; sliders have changed since the last time the theoretical system measures
; were calculated; if so, the theoretical measures are recalculated.

to update-theoretical-values
  if ((save-mean-arrival-rate != mean-arrival-rate)
      or (save-mean-service-time != mean-service-time)) [
    set save-mean-arrival-rate mean-arrival-rate
    set save-mean-service-time mean-service-time
    compute-theoretical-measures
  ]
end