Power Grid

extensions [vid]
breed [nodes node]
undirected-link-breed [edges edge]

breed [fires fire]

nodes-own [
  key
  available-power
  actual-power
  search-parent
  search-costs-from-start
]

edges-own [
  capacity
  load
  life-time
]

fires-own [
  fire-key1
  fire-key2
  fire-capacity
  fire-load
  fire-life-time
]

globals [mouse-first-node mouse-second-node search-open-node-list search-closed-node-list current-time is-grid-disrupted? is-grid-isolated? initial-number-of-edges _recording-save-file-name]

to setup-globals
  set-default-shape nodes "circle"
  set-default-shape fires "fire"
  set mouse-first-node nobody
  set mouse-second-node nobody
  set current-time 0
  set is-grid-disrupted? false
  set is-grid-isolated? false
  set initial-number-of-edges 0
end 

to startup
  setup
  let file "grid.dat"
  if file-exists? file [ read-grid-data-from-file file ]
end 

to setup
  clear-all
  setup-globals
  reset-ticks
  setup-fire
end 

to go
;  beep
;  user-message "The procedure 'go' is not implemented yet!"
;  stop
  if go-finished? [ stop ]
  ifelse ticks > 0
  [ go-failure-of-a-edge ]
  [
    set initial-number-of-edges count edges
    numerize-nodes
  ]
  optimize-flow
  go-calculate-life-time-of-edges
  if data-output != "no" [ write-data-to-files ]
  tick
end 

to-report go-finished?
  report is-grid-isolated? or count edges = 0 or count nodes = 0
end 

;; ============= Simulate network dynamics =================

to go-failure-of-a-edge
  if any? edges with [life-time > 0 ] [
    ask min-one-of edges with [life-time > 0 ] [life-time] [
      if output-level = 1 [ show "selected!" ]
      set current-time current-time + life-time
      go-create-fire-at-destroyed-edge
      die
    ]
  ]
end 

to go-calculate-life-time-of-edges
  ask edges [
    ifelse load > 0
      [ set life-time random-exponential capacity / load ]
      [ set life-time -1 ]
  ]
end 

to go-create-fire-at-destroyed-edge
  let x 0.5 * ( [xcor] of end1 + [xcor] of end2)
  let y 0.5 * ( [ycor] of end1 + [ycor] of end2)
  let key1 [key] of end1
  let key2 [key] of end2
  ask one-of fires [
    set xcor x
    set ycor y
    ifelse show-destroyed [ show-turtle ][ hide-turtle ]
    set fire-key1 key1
    set fire-key2 key2
    set fire-capacity [capacity] of myself
    set fire-load [load] of myself
  ]
end 

to setup-fire
  create-fires 1 [
    set size 2
    set color white
    hide-turtle
    set fire-key1 -1
    set fire-key2 -1
    set fire-capacity -1
    set fire-load -1
  ]
end 

;; ============= Draw nodes and edges ======================

to draw-structure
  ask nodes [ draw-node ]
  ask edges [ draw-edge ]
end 

to draw-node
  ifelse actual-power != 0 [ set size 1 + log abs actual-power 10 ] [ set size 1 ]
  let power available-power
  ifelse is-dangling? [ set color yellow ] [ set color green ]
  if power > 0 [ set color blue ]
  if power < 0 [ set color red ]
end 

to-report is-dangling?
  report count edges with [end1 = myself or end2 = myself] < 2
end 

to draw-edge
  ifelse load > 0
  [
    set color green
    set thickness (1 - exp (- 0.5 * load))
    if load > capacity
    [
      ifelse load > 2 * capacity
      [
        ifelse load > 4 * capacity
        [ set color red ]
        [ set color orange ]
      ]
      [ set color yellow ]
    ]
  ]
  [
     set color grey
     set thickness (1 - exp (- 0.5 * capacity))
  ]
end 

to numerize-nodes
  if any? nodes with [key < 0] [
    let counter 0
    ask nodes [
      set counter counter + 1
      set key counter
    ]
  ]
end 

to-report edge-mean
  report mean [count my-edges] of nodes
end 

to-report edge-variance
  report standard-deviation [count my-edges] of nodes
end 

;;============== Design the net-structure ===================

to add-nodes
  setup-nodes
  setup-spatially-clustered-network
  draw-structure
end 

to setup-nodes
  create-nodes number-of-nodes  [
    ; for visual reasons, we don't put any nodes *too* close to the edges
    setxy (random-xcor * 0.95) (random-ycor * 0.95)
    setup-node
  ]
end 

to setup-node
    set available-power 0
    set actual-power 0
    set key -1
end 

to setup-spatially-clustered-network
  let number-of-edges (average-node-degree * count nodes) / 2
  while [count edges < number-of-edges ]
  [
    ask one-of nodes [
      let choice ( min-one-of
        ( other nodes with [not edge-neighbor? myself] )
        [distance myself] )
      if choice != nobody [
        create-edge-with choice [
          set capacity 1
          set load 0
        ]
      ]
    ]
  ]
end 

to add-power
  let total-power total-power-flux
  while [total-power > 0] [
    let power 1
    if random-type = "float" [ set power power +   random-float 2 ]
    if random-type = "integer" [ set power power +   random 2 ]
    if power > total-power [ set power total-power ]
    set total-power total-power - power
    ask one-of nodes with [available-power <= 0] [
      set available-power available-power - power
    ]
    ask one-of nodes with [available-power >= 0] [
      set available-power available-power + power
    ]
  ]
  draw-structure
end 

to remove-power
  ask nodes [
    set available-power 0
    set actual-power 0
    draw-node
  ]
  ask edges [
    set load 0
    set capacity 1
    draw-edge
  ]
end 

to-report edge-level
  let value 100.
  if initial-number-of-edges > 0 [
    set value 100.0 * count edges / initial-number-of-edges
  ]
  report value
end 

to-report disrupted-level
  let value 0
  if is-grid-disrupted? [ set value 100 ]
  report value
end 

;;============== Editing the nodes and edges =========

to-report select-nearest-node
  report one-of nodes with-min [distancexy mouse-xcor mouse-ycor]
end 

to insert-node
  if mouse-down? and mouse-inside? [
    create-nodes 1 [
      setxy mouse-xcor mouse-ycor
      setup-node
      draw-node
    ]
    stop
  ]
end 

to delete-node
  if mouse-down? and mouse-inside? [
    ask select-nearest-node [ die ]
    draw-structure
    stop
  ]
end 

to move-node
  if mouse-down? and mouse-inside? [
    ask select-nearest-node [ setxy mouse-xcor mouse-ycor ]
    if not mouse-down? [ stop ]
  ]
end 

to set-selected
  set size size * 2
end 

to set-unselected
  set size size * 0.5
end 

to insert-edge
  if mouse-first-node = nobody [
     if mouse-down? and mouse-inside? [
       set mouse-first-node select-nearest-node
       ask mouse-first-node [ set-selected ]
     ]
  ]
  if mouse-second-node = nobody or mouse-second-node = mouse-first-node [
    if mouse-down? and mouse-inside? [ set mouse-second-node select-nearest-node ]
  ]
  if mouse-first-node != nobody and mouse-second-node != nobody
  and mouse-second-node != mouse-first-node [
    ask mouse-first-node [
      set-unselected
      create-edge-with mouse-second-node [
          set capacity 1
          set load 0
        ]
    ]
    draw-structure
    set mouse-first-node nobody
    set mouse-second-node nobody
    stop
  ]
end 

to-report edges-between [node1 node2]
  report edges with [end1 = node1 and end2 = node2]
end 

to delete-edge
  if mouse-first-node = nobody [
     if mouse-down? and mouse-inside? [
       set mouse-first-node select-nearest-node
       ask mouse-first-node [ set-selected ]
     ]
  ]
  if mouse-second-node = nobody or mouse-second-node = mouse-first-node [
    if mouse-down? and mouse-inside? [ set mouse-second-node select-nearest-node ] ]
  let found-edges edges-between mouse-first-node mouse-second-node
  if any? found-edges [
    ask one-of found-edges [ die ]
    ask mouse-first-node [ set-unselected ]
    set mouse-first-node nobody
    set mouse-second-node nobody
    draw-structure
    stop
  ]
end 

to increase-power
  if mouse-down? and mouse-inside? [
    ask select-nearest-node [
      set available-power available-power + 1
      draw-node
    ]
    stop
  ]
end 

to decrease-power
  if mouse-down? and mouse-inside? [
    ask select-nearest-node [
      set available-power available-power - 1
      draw-node
    ]
    stop
  ]
end 

to radial-layout
  if mouse-down? and mouse-inside? [
    layout-radial nodes edges select-nearest-node
    stop
  ]
end 

to spring-layout
  let spring-force 1
  let spring-length world-width / (sqrt count nodes)
  let repulsion-force 1
  repeat 30 [ layout-spring nodes edges spring-force spring-length repulsion-force ]
end 

to circle-layout
  let radius 0.4 * min (list world-width world-height)
  let node-set max-n-of 3 nodes [count edge-neighbors ]
  repeat 10 [ layout-tutte node-set edges radius ]
end 

;;============== Save and load net-structure ===================

to write-data-to-files
  if data-output = "each" or data-output = "all" [
    let file-name (word "grid-sim-" ticks)
    let network-file (word file-name ".dat")
    if is-string? network-file [
      if file-exists? network-file [ file-delete network-file ]
      write-grid-data-to-file network-file
    ]
  ]
  if data-output = "for R" or data-output = "all" [ write-data-to-R-files ]
end 

to write-data-to-R-files
  numerize-nodes
  let file-name (word "grid-nodes-R-" number-of-run ".dat")
  if ticks = 0 and file-exists? file-name [ file-delete file-name ]
  file-open file-name
  if ticks = 0 [ file-print "key current-time xcor ycor available-power actual-power" ]
  ask nodes [
    file-write key
    file-write current-time
    file-write xcor
    file-write ycor
    file-write available-power
    file-write actual-power
    file-print " "
  ]
  file-close
  set file-name (word "grid-edges-R-" number-of-run ".dat")
  if ticks = 0 and file-exists? file-name [ file-delete file-name ]
  file-open file-name
  if ticks = 0 [ file-print "key1 key2 current-time capacity load life-time deleted?" ]
  ask fires [
    if fire-key1 >= 0 [
      file-write fire-key1
      file-write fire-key2
      file-write current-time
      file-write fire-capacity
      file-write fire-load
      file-write fire-life-time
      file-print " 1 "
    ]
  ]
  ask edges [
    file-write [key] of end1
    file-write [key] of end2
    file-write current-time
    file-write capacity
    file-write load
    file-write life-time
    file-print " 0 "
  ]
  file-close
end 

to-report check-file-name [file-name file-tag]
  if is-string? file-name [
    let found substring file-name (length file-name - length file-tag) length file-name
    if found != file-tag [ set file-name (word file-name file-tag)     ]
  ]
  report file-name
end 

to save-grid
  let network-file check-file-name user-new-file ".dat"
  if is-string? network-file [
    if file-exists? network-file [ file-delete network-file ]
    write-grid-data-to-file network-file
  ]
end 

to write-grid-data-to-file [network-file]
  numerize-nodes
  file-open network-file
  file-print count nodes
  file-print "* node data key label x y available-power actual-power"
  foreach sort-on [key] nodes  [ [?1] ->
    ask ?1 [
      if empty? label [ set label (word key)]
      file-write key
      file-write label
      file-write xcor
      file-write ycor
      file-write available-power
      file-write actual-power
      file-print " "
    ]
  ]
  file-print "* edge data key1 key2 capacity load"
  ask edges [
    file-write [key] of end1
    file-write [key] of end2
    file-write capacity
    file-write load
    file-print " "
  ]
  file-close
end 

to load-grid
  setup
  let network-file user-file
  if is-string? network-file and file-exists? network-file [
    read-grid-data-from-file network-file
  ]
end 

to read-grid-data-from-file [network-file]
  file-open network-file
  let counter file-read
  let dummy file-read-line
  while [counter > 0] [
    create-nodes 1 [
      set color green
      set key file-read
      set label file-read
      setxy file-read file-read
      set available-power file-read
      set actual-power file-read
    ]
    set counter counter - 1
  ]
  set dummy file-read-line
  while [not file-at-end? ] [
    let token file-read
    let next-token file-read
    let first-node one-of nodes with [key = token]
    let second-node one-of nodes with [key = next-token]
    ask first-node [
      create-edge-with second-node [
        set capacity file-read
        set load file-read
      ]
    ]
  ]
  file-close
  draw-structure
end 

;;============== Export net-structure in various formats ===

to export
  numerize-nodes
  if format = "NET"
  [
    let network-file check-file-name user-new-file".net"
    if is-string? network-file [
      if file-exists? network-file [ file-delete network-file ]
      write-NET-data-to-file network-file
    ]
  ]
  if format = "VNA"
  [
    let network-file check-file-name user-new-file ".vna"
    if is-string? network-file [
      if file-exists? network-file [ file-delete network-file ]
      write-VNA-data-to-file network-file
    ]
  ]
  if format = "R"
  [
    let node-file check-file-name user-new-file ".nodes.imp"
    let edge-file (word (remove ".nodes.imp" node-file ) ".edges.imp")
    if is-string? node-file and is-string? edge-file [
      if file-exists? node-file [ file-delete node-file ]
      if file-exists? edge-file [ file-delete edge-file ]
      write-R-node-data-to-file node-file
      write-R-edge-data-to-file edge-file
    ]
  ]
end 

to write-NET-data-to-file [network-file]
  file-open network-file
  file-type "*Vertices " file-print count nodes
  foreach sort-on [key] nodes  [ [?1] ->
    ask ?1 [
      file-write key
      file-write label
      file-write xcor
      file-write ycor
      file-print " "
    ]
  ]
  file-print "*Arcs"
  ask edges [
    file-write [key] of end1 file-type " "
    file-write [key] of end2 file-type " "
    file-write 1 + load
    file-print " "
  ]
  file-close
end 

to write-VNA-data-to-file [network-file]
  file-open network-file
  file-print "*Node data"
  file-print "id available-power actual-power"
  foreach sort-on [key] nodes [ [?1] ->
    ask ?1 [
      if empty? label [ set label (word key)]
      file-write key
      file-write precision available-power 2
      file-write precision actual-power 2
      file-print " "
    ]
  ]
  let size-factor 10
  file-print "*Node properties"
  file-print "id x y color shape size shortlabel"
  let vshape 1
  foreach sort-on [key] nodes [ [?1] ->
    ask ?1 [
      file-write key
      file-write precision (size-factor * (xcor - min-pxcor)) 0
      file-write precision (size-factor * (ycor - min-pycor)) 0
      file-write integer-color
      file-write vshape
      file-write precision (size-factor * size) 0
      file-write label
      file-print " "
    ]
  ]
  file-print "*Tie data"
  file-print "from to strength load capapcity"
  ask edges [
    file-write [key] of end1
    file-write [key] of end2
    file-write 1
    file-write load
    file-write capacity
    file-print " "
    file-write [key] of end2
    file-write [key] of end1
    file-write 1
    file-write load
    file-write capacity
    file-print " "
  ]
  file-print "*Tie properties"
  file-print "from to color size headcolor headsize active"
  let headsize 0
  let active -1
  ask edges [
    let lcolor integer-color
    file-write [key] of end1
    file-write [key] of end2
    file-write lcolor
    file-write precision (size-factor * thickness) 0
    file-write lcolor
    file-write headsize
    file-write active
    file-print " "
    file-write [key] of end2
    file-write [key] of end1
    file-write lcolor
    file-write precision (size-factor * thickness) 0
    file-write lcolor
    file-write headsize
    file-write active
    file-print " "
  ]
  file-close
end 

to write-R-node-data-to-file [network-file]
  file-open network-file
  file-print "key x y available actual"
  ask nodes [
    file-write key
    file-write xcor
    file-write ycor
    file-write available-power
    file-write actual-power
    file-print " "
  ]
  file-close
end 

to write-R-edge-data-to-file [network-file]
  file-open network-file
  file-print "key1 key2 capacity load "
  ask edges [
    file-write [key] of end1
    file-write [key] of end2
    file-write capacity
    file-write load
    file-print " "
  ]
  file-close
end 

to-report integer-color
  let value 0
  let color-list extract-rgb color
  let red-value item 0 color-list
  let green-value item 1 color-list
  let blue-value item 2 color-list
  set value red-value + 256 * green-value + 256 * 256 * blue-value
  report value
end 

;;============== provide characteristic values for net-structure ==

to-report power-supply-level
  let value 0
  let total power-supply
  if total > 0 [
    set value 100.0 * power-input / total
  ]
  report value
end 

to-report power-demand-level
  let value 0
  let total power-demand
  if total > 0 [
    set value 100.0 * power-output / total
  ]
  report value
end 

to-report power-supply
  let power 0
  ask nodes with [available-power > 0] [
    set power power + available-power
  ]
  report power
end 

to-report power-demand
  let power 0
  ask nodes with [available-power < 0] [
    set power power + available-power
  ]
  report power * -1
end 

to-report power-input
  let power 0
  ask nodes with [available-power > 0] [
    set power power + actual-power
  ]
  report power
end 

to-report power-output
  let power 0
  ask nodes with [available-power < 0] [
    set power power + actual-power
  ]
  report power * -1
end 

to-report power-variance
  let value 0
  ask nodes [
    let delta available-power - actual-power
    set value value + delta * delta
  ]
  report sqrt value
end 

to-report flux-variance
  let value 0
  ask edges [
    let delta capacity - load
    set value value + delta * delta
  ]
  report sqrt value
end 

;; ================= Optimize the flux in net-structure ==========

to reset-structure
  ask nodes [
    set actual-power 0
    draw-node
  ]
  ask edges [
    set load 0
    draw-edge
  ]
end 

to optimize-flow
  reset-structure
  let is-isolated? true
  let is-disrupted? false
  let targets nodes with [available-power < 0]
  ask targets [
    let sources nodes with [can-provide?]
    ask min-n-of (count sources) sources [distance myself] [
      if [is-needing?] of myself [
        set is-disrupted? not update-net-structure self myself
        set is-isolated? is-isolated? and is-disrupted?
        set is-grid-disrupted? is-disrupted? or is-grid-disrupted?
      ]
    ]
  ]
  set is-grid-isolated? is-isolated?
  draw-structure
end 

to-report is-needing?
  report available-power < 0 and actual-power > available-power
end 

to-report can-provide?
  report available-power > 0 and actual-power < available-power
end 

to change-power [this-load]
  if available-power > 0 [ set actual-power actual-power + this-load ]
  if available-power < 0 [ set actual-power actual-power - this-load ]
end 

to-report calculate-net-flow [start-node target-node]
  let this-flow  0
  ask start-node [
    set this-flow available-power - actual-power
  ]
  ask target-node [
    let that-flow actual-power - available-power
    if that-flow < this-flow [ set this-flow that-flow ]
  ]
  report this-flow
end 

to change-flow-structure [start-node target-node edge-list this-load]
  if not empty? edge-list [
    ask start-node [
      change-power this-load
      draw-node
    ]
    ask target-node [
      change-power this-load
      draw-node
    ]
    foreach edge-list [ [?1] ->
      ask ?1 [
        set load load + this-load
        draw-edge
      ]
    ]
  ]
end 

to-report update-net-structure [start-node target-node]
  let path-found? true
  let edge-list search-go start-node target-node
  let found-path? not empty? edge-list
  ifelse found-path?
  [
    let this-load calculate-net-flow start-node target-node
    if output-level = 1 [ show (word "Power flow: " this-load " with " (length edge-list)
      " edges between " start-node " and " target-node) ]
    if this-load > 0 [ change-flow-structure start-node target-node edge-list this-load ]
  ]
  [
    if output-level = 1 [ show (word "No path found between " start-node " and " target-node) ]
    set path-found? false
  ]
  report path-found?
end 

to adjust-capacities
  ask edges [
    if load > capacity [ set capacity load ]
    draw-edge
  ]
end 

to remove-unused
  ask edges with [load = 0] [ die ]
  ask nodes with [ count my-edges = 0] [die]
  draw-structure
end 

to optimize-structure
  repeat optimize-steps [
    optimize-flow
    adjust-capacities
  ]
  optimize-flow
end 

;; ================= Search shortest path in net-structure =======

to-report search-go [start-node target-node]
  let node-list search-path start-node target-node
  let edge-list search-transfer-node-list-to-edge-list node-list
  report edge-list
end 

to-report search-path [start-node target-node]
  if output-level = 2 [ show ( word "Search path between " start-node " and " target-node ) ]
  let new-path ( list )
  search-init start-node
  if search-do target-node [
    set new-path search-path-back target-node
  ]
  report new-path
end 

to search-init [ start-node ]
  if output-level = 2 [ show ( word "Init search from " start-node ) ]
  ask nodes [
    set search-parent nobody
    set search-costs-from-start 0
  ]
  set search-open-node-list fput start-node ( list )
  set search-closed-node-list ( list )
end 

to-report search-rank
  report search-costs-from-start
end 

to-report search-do [target-node]
  if output-level = 2 [ show ( word "Do search to " target-node ) ]
  let current-node nobody
  while [target-node != current-node] [
    if empty? search-open-node-list [
      if output-level = 2 [ show ( word "No path to " target-node ) ]
      report false
    ]
    ; remove lowest rank item from open list of patches and add it to the closed list
    set search-open-node-list sort-by [ [?1 ?2] -> [ search-rank ] of ?1 < [ search-rank ] of ?2 ] search-open-node-list
    set current-node first search-open-node-list
    set search-open-node-list but-first search-open-node-list
    set search-closed-node-list fput current-node search-closed-node-list
    if output-level = 2 [ show ( word "Current node " current-node ) ]
    ; check adjacent nodes
    if target-node != current-node [
      ask current-node [ search-handle-neighbors self target-node]
    ]
  ]
  if output-level = 2 [ show ( word "Found target " current-node ) ]
  report true
end 

to search-handle-neighbors [parent-node target-node]
  ask my-edges [
    let costs [ search-costs-from-start ] of parent-node + 1
    if load > capacity [ set costs costs + over-capacity-costs ]
    ask other-end [
      if member? self search-open-node-list and costs < search-costs-from-start [
        set search-open-node-list remove self search-open-node-list
        if output-level = 2 [ show ( word "Neighbor node " self
            " removed from open " search-open-node-list ) ]
      ]
      if member? self search-closed-node-list and costs < search-costs-from-start [
        set search-closed-node-list remove self search-closed-node-list
        if output-level = 2 [ show ( word "Neighbor node " self
            " removed from closed " search-closed-node-list ) ]
      ]
      if ( not member? self search-open-node-list )
      and ( not member? self search-closed-node-list ) [
        if output-level = 2 [ show ( word "Neighbor node " self
            " with costs=" costs " to parent " parent-node ) ]
        set search-parent parent-node
        set search-costs-from-start costs
        set search-open-node-list fput self search-open-node-list
      ]
    ]
  ]
end 

to-report search-path-back [target-node]
  let found-path fput target-node ( list )
  let current-node target-node
  if output-level = 2 [ show ( word "Revert search " current-node ) ]
  while [ [ search-parent ] of current-node != nobody ] [
    set current-node [ search-parent ] of current-node
    set found-path fput current-node found-path
    if output-level = 2 [ show ( word "Revert search " current-node ) ]
  ]
  report found-path
end 

to-report search-edge-for-nodes [that-node this-node]
  report one-of edges with [ (end1 = that-node and end2 = this-node)
    or (end2 = that-node and end1 = this-node) ]
end 

to-report search-transfer-node-list-to-edge-list [node-list]
  let edge-list (list)
  let last-node nobody
  foreach node-list [ [?1] ->
    let current-node ?1
    if last-node != nobody [
      let found-edge search-edge-for-nodes current-node last-node
      if output-level = 2 [ show (word "Found " found-edge " of " current-node " " last-node) ]
      set edge-list fput found-edge edge-list
    ]
    set last-node current-node
  ]
  report edge-list
end 

;; ============== plotting ===========================

to plot-histogram-of [that-distribution]
  ifelse length that-distribution > 0 [
    let x-max ( ceiling max that-distribution )
    if x-max <= 0 [ set x-max 1.0 ]
    let y-max length that-distribution
    set-plot-x-range 0 x-max
    set-plot-y-range 0 y-max
    set-histogram-num-bars 20
    histogram that-distribution
  ] [
    clear-plot
  ]
end