Collaborative demand forecasting

;; forecasting is done in R environment
;; Before using the model it is important to install NetLogo 5.3.1, R code above 3.2, extension R (if it is not present) and the newest rJava package in R
;; Configuration  on windows environment variables:
;; create new variable - name: R_HOME, path: C:\Program Files\R\R-3.3.3
;; create new variable - name: JRI_HOME, path: C:\Users\user\Documents\R\win-library\3.3\rJava\jri
;; Add to existing variable  - name: Path, path: C:\Program Files\R\R-3.3.3\bin\x64
;; More information can be found in: https://ccl.northwestern.edu/netlogo/docs/r.html
extensions [r]

globals [
  first_demand ;; create a first demand variable to generate initial inputs and prepare demand history for forecasting
  mean_demand ;; distributions are used from multiple product category, Therefore there is extreme variation in the demand. The mean_demand is used as a threshold to limit this variation.
        ]

breed [distributors distributor]
breed [households household]

distributors-own [
  order_to_supplier
  order_from_household
  total_order_from_household
  available_stock
  demand_history ;; a list of all previous demand history
  group_demand_history ;; list of lists of all distributor demand history
  forecasting ;; variable to save forecasting result with all statistics
  forecast ;; variable for final forecasted value
  forecasting_history ;;save forecasting values to history
  plot_forecast ;; save previous tick forecasting result for plotting with delay
  backlog ;; level of unfulfilled demand
  total_backlog ;; level of unfulfilled demand for all distributors
  unfulfilled_demand ;; the backlog level is added to unfulfilled demand and is necessary to deliver once available inventory is present
  safety_stock ;; safety stock level based on previous 7 tick demand average and safety_level set by slider
]

households-own [
  order_to_distributor
  distribution_type ;; randomly assigns distribution type to select a distribution function
  product_quantity ;; the distributions are generated from ratio of sales by day. Therefore, the generated ratio is multiplied by different product quantity, to simulate the variation of consumer shopping habits.
]

to setup
  __clear-all-and-reset-ticks
  ask patches [set pcolor white]
  ;; assign the initial demand of all household
  random-seed 100

  create-distributors Distributor_N [
    set size 5
    set color blue
    set label (word "Distr." [who] of self "     ")
    set forecast 0
    set forecasting_history []
  ]

  create-households Households_N [
   set size 1
   set color 3
   set order_to_distributor 0
   set distribution_type []
  ]

  set-default-shape distributors "house two story"
  set-default-shape households "person"

  ask turtles [find_patch]
end 

to go
  reset_variables

  if ticks = 0 [ set_first__variables ]
  ifelse consumer_loyalty = FALSE
  [
    create-link
    interaction_supplier
  ]
  [
    if ticks = 0 [ create-link ]
    interaction_supplier
  ]

  receive_order_from_supplier
  select_demand_distribution
  interaction_household
  check_mean_demand
  update_demand_history
  set_group_history
  update_inventory

  if consumer_loyalty = FALSE [ ask links [die] ]

  tick
end 

;; randomly assign all distributors and household to individual patch

to find_patch
  setxy random-xcor random-ycor
  while [any? other turtles-here]
      [find_patch]
end 

;; save individual distributor sales history to a list of lists, where the first element in the index of the distributor
;; and the second is a list of demand history of that distributor

to set_group_history
  ask turtles with [ breed = distributors ]
  [
  set group_demand_history (list (list who) (list demand_history))
  ]
end 

to set_first__variables
  select_demand_distribution
  ask turtles with [breed = households]
    [
      ;; the distribution functions are obtained from ratio of real-data. The purpose of this is to simulate seasonality of the real-world environment.
      ;; the ratio is then multiplied by the product quantity which household buy. This way the customer is also distinguished by number of people in the household.
      set product_quantity random (Max_product_quantity - 1) + 1
      set order_to_distributor round product_quantity * order_to_distributor
      if order_to_distributor < 0 [ set order_to_distributor (- order_to_distributor)]
      let demand []
      set demand order_to_distributor
      set demand round demand
      create-link
      ask link-neighbors with [ breed = distributors ] [
        set order_from_household demand
        set total_order_from_household (total_order_from_household + order_from_household)
      ]
    ]

  ask turtles with [ breed = distributors ] [
    set first_demand total_order_from_household
    set demand_history (list (first_demand) (first_demand) (first_demand) (first_demand) (first_demand))
    set group_demand_history (list (list who) (list demand_history))
    set available_stock first_demand * 3
     ]
end 

;; function to get first n elements of a list

to-report take_first [n xs]
  report sublist xs 0 min list n (length xs)
end 

;; the order from supplier is added to the available stock. The value is obtained from previous tick, therefore there is a delay of delivery of 1 day

to receive_order_from_supplier
  ask turtles with [ breed = distributors ]
  [
    set safety_stock demand_history
    set safety_stock reverse safety_stock
    set safety_stock filter [? != 0] take_first 7 safety_stock
    set safety_stock (mean safety_stock) * safety_level
    ifelse available_stock < safety_stock
    [ set available_stock (available_stock + order_to_supplier + safety_stock) ]
    [ set available_stock (available_stock + order_to_supplier) ]
  ]
end 

;; the demand is forecasted based on individual history or group history of companies. The distributor is ordering supplies 1 tick before the household orders.

to interaction_supplier
   ifelse Share_information = false
     ;; individual demand forecasting
     [
       let distributors_index n-values count  turtles with [breed = distributors] [?]
       ask turtles with [breed = distributors] [
         ;; need to iterate trough individual turtle, otherwise demand_history is added for all distributors not only 1 (list of lists problem)
         foreach distributors_index [

         ;; machine learning/neural network algorithm for individual forecasting
         ;; clear the R workspace
         r:clear
         ;; command to clear console when printing variables use r:eval "print(variable_name)" or show r:get "variable_name", but some formats are not supported by NetLogo
         r:eval "cat(rep('\n',50))"
         ;; load extreme machine learning package
         r:eval "library('elmNN')"
         r:eval " splitWithOverlap <- function(vec, seg.length) { overlap = seg.length - 1 ; starts = seq(1, length(vec), by=seg.length-overlap); ends   = starts + seg.length - 1 ;ends[ends > length(vec)] = length(vec) ;y <- lapply(1:length(starts), function(i) vec[starts[i]:ends[i]]) ;y <- y[unlist(lapply(y, length) == seg.length)] ;y <- as.data.frame(t(as.data.frame(y))) ;rownames(y) <- c() ;colnames(y)[length(colnames(y))] <- 'y' ;return(y) } "

         ;; put demand_history variable in to R environment and create a sequence column to identify the time series
         r:put "r_demand_history" [demand_history] of turtle ?
         r:eval "train <- splitWithOverlap(r_demand_history,4)"
         ;; in the beginning the neural network when facing with spikes over predicts. Manually we set the over predictions to the mean. The threshold is decided by input
         r:eval "mean_demand <- mean(r_demand_history)"
         r:put "mean_threshold" mean_threshold

         ;; train the extreme machine learning algorithm
         r:eval "model <- elmtrain(y ~ . , data=train, nhid=1, actfun='sig')"
         ;; create an empty variable for the next tick prediction
         r:eval "predict_df <- tail(r_demand_history,3)"
         r:eval "predict_df <- as.data.frame(t(as.data.frame(predict_df)))"
         r:eval "predict_df$y <- 0"
         r:eval "colnames(predict_df) <- c('v1','v2','v3','y')"

         ;; predict the next tick's demand
         r:eval "prediction <- predict.elmNN(model, newdata = predict_df)"
         r:eval "prediction <- as.numeric(prediction)"
         ;; check if the predicted demand is not above the mean with a certain threshold
         if mean_limit_status = true [ r:eval "if (prediction > (mean_demand * mean_threshold)) { prediction <- mean_demand * mean_threshold}" ]
         ;; save the demand to NetLogo variable
         set forecasting r:get "prediction"
         ;; save the forecasting results and prepare order from supplier
         set forecast round forecasting
         set order_to_supplier forecast
         if order_to_supplier < 0 [set order_to_supplier 0]
         ifelse order_to_supplier + unfulfilled_demand < available_stock
         [ set order_to_supplier 0 ]
         [ set order_to_supplier (order_to_supplier + unfulfilled_demand - available_stock) ]
         ]
       ]
     ]
      ;; collaborative demand forecasting (forecast with machine learning from all distributor information, not only individual)
      [
        ;; machine learning/neural network algorithm for individual forecasting
        ;; clear the R workspace
        r:clear
        ;; command to clear console when printing variables use r:eval "print(variable_name)" or show r:get "variable_name", but some formats are not supported by NetLogo
        r:eval "cat(rep('\n',50))"
        ;; load extreme machine learning package
        r:eval "library('elmNN')"
        r:eval "splitWithOverlap <- function(vec, seg.length) { overlap = seg.length - 1 ; starts = seq(1, length(vec), by=seg.length-overlap); ends   = starts + seg.length - 1 ;ends[ends > length(vec)] = length(vec) ;y <- lapply(1:length(starts), function(i) vec[starts[i]:ends[i]]) ;y <- y[unlist(lapply(y, length) == seg.length)] ;y <- as.data.frame(t(as.data.frame(y))) ;rownames(y) <- c() ;colnames(y)[length(colnames(y))] <- 'y' ;return(y) } "
        r:eval "splitWithOverlap_combine <- function(x, seg.length = 4) { ;for (i in 1:ncol(x)) { ;active_df <- x[i] ;active_df <- t(as.vector(active_df)) ;ts <- splitWithOverlap(active_df,seg.length) ;colnames(ts) <- paste(colnames(x)[i], colnames(ts), sep = '_') ;if (exists('ts_final') == FALSE) { ;ts_final <- ts ;} else { ;ts_final <- cbind(ts_final,ts) ;} ;} ;return(ts_final) ;}"
        r:eval "create_predict_df <- function(x) { ;for (i in 1:ncol(x)) { ;active_df <- x[i] ;active_df <- tail(active_df,3) ;active_df <- as.data.frame(t(as.vector(active_df))) ;colnames(active_df) <- c('V1','V2','V3') ;active_df$y <- 0 ;colnames(active_df) <- paste(rownames(active_df), colnames(active_df), sep = '_') ;if (exists('df_final') == FALSE) { ;df_final <- active_df ;} else { ;df_final <- cbind(df_final,active_df) ;} ;} ;return(df_final);}"

        (r:putagent "r_demand_history" turtles with [breed = distributors] "group_demand_history")

        r:eval "r_demand_history <- as.data.frame(matrix(data = as.vector(unlist(r_demand_history)), nrow = length(unlist(r_demand_history$group_demand_history[[1]])), ncol = length(r_demand_history[[1]])))"
        r:eval "colnames(r_demand_history) <- r_demand_history[1, ]"
        r:eval "r_demand_history <- r_demand_history[-1, ] "
        r:eval "row.names(r_demand_history) <- NULL"
        r:eval "colnames(r_demand_history) <- paste('Distributor', colnames(r_demand_history), sep = '_')"
        r:eval "column_names <- colnames(r_demand_history)"
        r:eval "train <- splitWithOverlap_combine(r_demand_history,4)"

        ;; in the beginning the neural network when facing with spikes over predicts. Manually we set the over predictions to the mean. The threshold is decided by input
        r:eval "mean_demand <- mean(colMeans(r_demand_history[1:(ncol(r_demand_history)-1)]))"
        r:put "mean_threshold" mean_threshold

        ;; train the extreme machine learning algorithm
        ;;create a variable of the distributor index
        let distributors_index n-values count  turtles with [breed = distributors] [?]

        foreach distributors_index [

        r:put "active_distributors_index" ?
        r:eval "active_distributors_index <- as.numeric(active_distributors_index)"
        r:eval "active_distributor <- paste0('Distributor_',active_distributors_index, '_y')"
        r:eval "formula <- paste(active_distributor, ' ~ ' ,'.')"
        r:eval "model <- elmtrain(as.formula(formula) , data=train, nhid=1, actfun='sig')"

        ;; create an empty variable for the next tick prediction
        r:eval "col_names <- colnames(train)"
        r:eval "predict_df <- create_predict_df(r_demand_history)"

        ;; predict the next tick's demand
        r:eval "prediction <- predict.elmNN(model, newdata = predict_df)"
        r:eval "prediction <- as.numeric(prediction)"
        ;; check if the predicted demand is not above the mean with a certain threshold
        if mean_limit_status = true [ r:eval "if (prediction > (mean_demand * mean_threshold)) { prediction <- mean_demand * mean_threshold}" ]

        ;; save the demand to NetLogo variable
        ask turtle ? [
        set forecasting r:get "prediction"
        ;; save the forecasting results and prepare order from supplier
        set forecast round forecasting
        set order_to_supplier forecast
        if order_to_supplier < 0 [set order_to_supplier 0]
        ifelse order_to_supplier + unfulfilled_demand < available_stock
          [ set order_to_supplier 0 ]
          [ set order_to_supplier (order_to_supplier + unfulfilled_demand - available_stock) ]
        ]
        ]
      ]
end 

;; pareto distribution (pareto from R package: tdistrplus)

to-report random-pareto [alpha mm]
  report mm / ( random-float 1 ^ (1 / alpha) )
end 

;; lognormal distribution (lnorm from R package: tdistrplus)

to-report  log-normal [mu sigma]
  let beta ln (1 + ((sigma ^ 2) / (mu ^ 2)))
  let x exp (random-normal (ln (mu) - (beta / 2)) sqrt beta)
  report x
end 

;; log-logistics distribution (llogis from R package: tdistrplus)

to-report random-loglogistics [alpha beta]
    let r random-float 1
    let q r / (1 - r)
    let x alpha * exp (ln (q) / beta)
    report x
end 

;; burr distribution (burr from R package: tdistrplus)

to-report random-burr [c k]
    let r random-float 1
    let q exp ( - ln(r) / k) - 1
    ;;if q < 0 [ set q (- q) ]
    let x exp ( ln(q) / c )
    report x
end 

to select_demand_distribution
     ask turtles with [breed = households]
        [
     ;; the household decides how much to order based on distributions fitted on real-data.
     ;; if the market type is perfect competition it means that there are distributors who are not involved in sharing information
     ;; therefore some household randomly buy products from other distributors
     ;; if the market is oligopoly it means that all distributors from the market are involved in sharing information
     if market_type = "Oligopoly" [set distribution_type random (7 - 1) + 1]
     if market_type = "Perfect competition" [set distribution_type random (28 - 1) + 1]
     if distribution_type = 1  [ set order_to_distributor random-pareto 8.021208 10.829160  ]
     if distribution_type = 2  [ set order_to_distributor random-loglogistics 0.5279758 1.0251284 ]
     if distribution_type = 3  [ set order_to_distributor random-burr 0.8406488 2.1144042]
     if distribution_type = 4  [ set order_to_distributor random-burr 0.888256 1.807287]
     if distribution_type = 5  [ set order_to_distributor random-pareto 8.021208 10.829160 ]
     if distribution_type = 6  [ set order_to_distributor random-loglogistics 0.5594916 0.9898947 ]
     if distribution_type = 7  [ set order_to_distributor log-normal 0.5594916 0.9898947 ]
     if distribution_type > 7  [ set order_to_distributor 0 ]
   ]
end 

to create-link
   ask turtles with [breed = households] [create-link-with one-of other turtles with [ breed = distributors ] ]
end 

;;households randomly selects a distributor and sends him an order based on demand distribution type

to interaction_household
   ask turtles with [breed = households]
   [
     ;; the distribution functions are obtained from ratio of real-data. The purpose of this is to simulate seasonality of the real-world environment.
     ;; the ratio is then multiplied by the product quantity which household buy. This way the customer is also distinguished by number of people in the household.
     set product_quantity random (Max_product_quantity - 1) + 1
     if order_to_distributor < 0 [ set order_to_distributor (- order_to_distributor)]
     set order_to_distributor round product_quantity * order_to_distributor
     let demand []
     set demand order_to_distributor
     set demand round demand

     ask link-neighbors with [ breed = distributors ] [

       set order_from_household demand
       set total_order_from_household (total_order_from_household + order_from_household)
     ]

   ]
end 

;; determines the mean demand

to check_mean_demand
  if mean_limit_status = True [
     ask turtles with [breed = distributors]
   [
     let distributors_index n-values count  turtles with [breed = distributors] [?]
     foreach distributors_index [
     r:put "mean_demand" [demand_history] of turtle ?
     r:eval "mean_demand <- mean(mean_demand)"
     r:eval "mean_demand <- as.numeric(mean_demand)"
     set mean_demand r:get "mean_demand"
     if mean_limit_status = true [ if total_order_from_household > mean_demand * mean_threshold [set total_order_from_household mean_demand * mean_threshold] ]
     ]
   ]
   ]
end 

;; send orders to household by subtracting total order from household from available stock

to update_inventory
   ask turtles with [breed = distributors]
     [
       ifelse unfulfilled_demand != 0
       [
         set available_stock (available_stock - total_order_from_household - unfulfilled_demand)
         set unfulfilled_demand 0
       ]
       [
         set available_stock (available_stock - total_order_from_household)
       ]

       ;; evaluate  the unfulfilled orders if the stock level was negative
       if available_stock < 0
      [
        set backlog (- available_stock)
        set total_backlog total_backlog + backlog
        set unfulfilled_demand total_backlog
      ]
     ]
end 

;; add a total demand variable to the list of demand history

to update_demand_history
   ask turtles with [breed = distributors]
     [
       set demand_history lput total_order_from_household demand_history
       set forecasting_history lput forecast forecasting_history
       if ticks != 0 [
       set plot_forecast item (length(forecasting_history) - 2) forecasting_history
       ]
     ]
end 

;; reset save variables

to reset_variables
   ask turtles with [ breed = distributors ]
   [
       set total_order_from_household 0
       set total_backlog 0

   ]
end