Business Intelligence (BI) software

Software packages covering most analytics task in a customized and modular way. Always contains data warehouse concepts as well as reporting and data analysis routines. Famous firms are IBM, Oracle, Microsoft, SAP, SAS, …

Definition

R is an open-source, interpreted (statistical) programming language and environment available at (www.r-project.org).

Frontend(s)

R is the name of a programming language (like Java) as well as programming environment (like Eclipse). Traditionally, R is just console and an editor for writing scripts. However, most users work with Rstudio as a more powerful evironment for R (www.rstudio.com)

• e-commerce sales data from www.olist.com
• 99441 orders of 32951 products
• data from 9-2016 until 10-2018

Data obtained from https://www.kaggle.com/olistbr/brazilian-ecommerce

Order data

Ordered items

Product data

Geoloc data

• Create an account for kaggle & login
• Search for the Brazilian E-Commerce project
• donwload zip-file including 9 csv-files \(\rightarrow\) unzip
• import csv-files in Rstudio via “import dataset” (upper-right panel) \(\rightarrow\) “From text”

Order data

head(orders.data, 4)
>                           order_id                      customer_id
> 1 e481f51cbdc54678b7cc49136f2d6af7 9ef432eb6251297304e76186b10a928d
> 2 53cdb2fc8bc7dce0b6741e2150273451 b0830fb4747a6c6d20dea0b8c802d7ef
> 3 47770eb9100c2d0c44946d9cf07ec65d 41ce2a54c0b03bf3443c3d931a367089
> 4 949d5b44dbf5de918fe9c16f97b45f8a f88197465ea7920adcdbec7375364d82
>   order_status order_purchase_timestamp   order_approved_at
> 1    delivered      2017-10-02 10:56:33 2017-10-02 11:07:15
> 2    delivered      2018-07-24 20:41:37 2018-07-26 03:24:27
> 3    delivered      2018-08-08 08:38:49 2018-08-08 08:55:23
> 4    delivered      2017-11-18 19:28:06 2017-11-18 19:45:59
>   order_delivered_carrier_date order_delivered_customer_date
> 1          2017-10-04 19:55:00           2017-10-10 21:25:13
> 2          2018-07-26 14:31:00           2018-08-07 15:27:45
> 3          2018-08-08 13:50:00           2018-08-17 18:06:29
> 4          2017-11-22 13:39:59           2017-12-02 00:28:42
>   order_estimated_delivery_date
> 1           2017-10-18 00:00:00
> 2           2018-08-13 00:00:00
> 3           2018-09-04 00:00:00
> 4           2017-12-15 00:00:00

Installation via menu or install.packages() functions. Be aware, this may take a while. This presentation can be found here

install.packages(c("tmap","sf", "tmaptools", "brazilmaps", "png", "plotly", "geosphere", "htmlwidgets","metR",
            "tidyverse", "lubridate", "widgetframe", "ROI","ROI.plugin.glpk", "ompr", "ompr.roi", "vis.network"))
library(tmap)
library(sf)
library(tmaptools)
library(brazilmaps)
library(png)
library(plotly)
library(geosphere)
library(htmlwidgets)
library(metR)
library(tidyverse)
library(lubridate)
library(widgetframe)
library(ROI)
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)
library(visNetwork)

Objects & workspace

Everything in R is an object. Objects are always in a environment and are categorized in numerous ways, e.g. by storage mode, or class. The most important is the global environment, also called workspace. We will distinguish data objects and functions. Functions are objects which create environemnts.

# shows current environement
environment()
# show all active environments
search()
# assign an object called 'da.1'
da.1 <- 5
# define a function
f.1 <- function() print("Hello world!")  
# function printing the environment
f.2 <- function(){
  print(environment()) 
} 
# objects in the workspace
objects()

environment()
> <environment: R_GlobalEnv>

f.1()
> [1] "Hello world!"

f.2() 
> <environment: 0x00000000378a3d68>

Definition

Functions are code-wrappers. They consist of argument, body, and environment. The environment determines where function resides (i.e., from where it can be called). The body describes what happens to the arguments.

# show arguments of a function
formals(f.1) 
> NULL
# show body of a function
body(f.1) 
> print("Hello world!")
# show environment of a function
environment(f.1) 
> <environment: R_GlobalEnv>
# function with arguments
f.3 <- function(x,y) x+y
formals(f.3) # more to be seen
> $x
> 
> 
> $y

# most built-in functions have neither 
# attribute nor body --> external code
formals(sum);body(sum);environment(sum)
> NULL
> NULL
> NULL
# Dynamic lookup:
# In case an object is missing, R looks 
# in the functions' environment before
# reporting an error
x <- 5
y <- 7
f.4 <- function() x+y 
f.4()
> [1] 12

ggplotly(ggplot(products.data, aes(x=product_category_name)) + geom_bar() +  
  xlab("product category") + ylab("number of products") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7)))

Definition

Basic data objects (vectors,matrices and arrays) consist of scalars of the same storage type only (i.e. only numbers, characters, logicals) while advanced data objects (data frames and lists) can contain different types.

Name	Dim.	Built function
vector	1	c(),numeric()
matrix	2	matrix()
array	n	array()
data frame	2	data.frame()
list	1	list()

devopedia.org (2019)

Data objects are indexed by position, name, or logical identifier. Vectors are indexed by <vector>[<index>]. You can also negate the selection, i.e. extract all entries except the indexed entries by <vector>[-<index>]. Matrices are 2-dimensional and, thus, require 2 indices, i.e. <matrix>[<index row>, <index col.>]

x <- c(5, 6, 7)
# extract 1st entry 
x[1]              
> [1] 5
# extract 1st & 3rd entry
x[c(1,3)]         
> [1] 5 7
# logical index have same length as indexed vector 
x[c(T,F,T)]       
> [1] 5 7
# generate an indexing vector by logic expr.
y <- x <= 6       
# retain inices of those entries 
z <- which(y)     
# subvector of all entries of x which are > 6
x[!y]              
> [1] 7

x <- cbind(1:6, 2:7, 3:8) # column-wise binding
# first row, first column (just one element)
x[1,1]
> [1] 1
# first row (a vector)
x[1,]
> [1] 1 2 3
# logical index vector
y <- rep( x = c(T,F,T), times = 2) 
# mix row and column selections
x[y,c(1,3)]               
>      [,1] [,2]
> [1,]    1    3
> [2,]    3    5
> [3,]    4    6
> [4,]    6    8

All advanced data objects can have named elements. These namescan be used to extract the corresponding element by <object>$<element name>. Alternatively, elements can be accessed by <object>[[<element index>]] (single element) or <object>[<element index>] (multiple elements). Additionally, tibbles and data fames can be indexed like matrices

df1 <- data.frame(one = 1:2, two = c("nest","test") )
# access column one
df1$one           
> [1] 1 2
# a tibble (basically the same as data frame)
tb1 <- tibble(one = 1:2, two = c("nest","test")) 
tb1$two           # works too
> [1] "nest" "test"
# a list
l1 <- list(a = 1:3, b = "nest", c = TRUE)
#  for multiple components
l1[c(1,3)]
> $a
> [1] 1 2 3
> 
> $c
> [1] TRUE

# for multiple components by name
l1[c("a","c")]
> $a
> [1] 1 2 3
> 
> $c
> [1] TRUE
# names can be retrieved 
names(l1)
> [1] "a" "b" "c"
# ... and changed
names(l1) <- c("x","y","z")
l1$z
> [1] TRUE

# sum revenues (prices) per product over all orders 
df.prod.cat <- aggregate(price ~ product_id, FUN = sum, data = orders.items)
# add corresponding product categories
df.prod.cat <- merge(df.prod.cat, products.data, by = "product_id")
# create plot
ggplot(df.prod.cat, aes(x=product_category_name, y=log(price)) ) + geom_boxplot(outlier.size = 0.01) +   ylab("log. prices") +
  xlab("Product category") + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7))

Definition

Computer simulation means to model a real-world system by a computer model. Thereby, the modelled entities and their relations among each other as well as external factors are usually described by mathematical expressions.

Naturally, R is often useful for stochastic simulation as it allows one to easily model stochastic processes. One core element are distribution functions. For a given distribution R offers usually the density (d), probability (p), quantile (q) and random number generator (r). The syntax for the normal distribution is e.g. d/p/q/rnorm(x). See ?Distributions

# sample from a vector
sample(1:10, 5)
> [1]  6 10  2  7  5
# .. also with relacements
sample(1:5, 5, replace = TRUE)
> [1] 4 1 2 5 2
# 95%-qunatile of the st. normal dist.
qnorm(0.95)
> [1] 1.644854
# 5 Bernoulli trial with prob. 0.5
rbinom(n=5, size=1 , prob=0.5)
> [1] 0 0 0 1 0
# ... same with replicate
replicate(5, rbinom(n=1, size=1 , prob=0.5))
> [1] 1 0 0 1 0

# most often ordered product
x.max <- rownames(df.freq.prod)[1]
# identify orders with x.max
orders.items.x.max <- orders.items[which(orders.items$product_id == x.max), ]
# add order information
orders.x.max <- merge(orders.items.x.max, orders.data, by="order_id")
# reorder according to date
orders.x.max <- orders.x.max[order(orders.x.max$order_purchase_timestamp), ]
# convert dates 
orders.x.max$order_purchase_timestamp <- as_datetime(orders.x.max$order_purchase_timestamp)
# add date components and order quantity column
orders.x.max <- data.frame(orders.x.max, 
          year.order = year(orders.x.max$order_purchase_timestamp),
          day.order = yday(orders.x.max$order_purchase_timestamp),
          ord.quant = rep(1,nrow(orders.x.max)))
# convert days in factor 
orders.x.max$day.order <- as.factor(orders.x.max$day.order)
# add missing factor levels
levels(orders.x.max$day.order) <- c(
  levels(orders.x.max$day.order),
  as.character(1:365)[!as.character(1:365) %in% levels(orders.x.max$day.order)])
# sort factor levels
orders.x.max$day.order <- factor(orders.x.max$day.order,levels(orders.x.max$day.order)[order(as.numeric(levels(orders.x.max$day.order)))])

# calculate daily demand
day.dem.x.max <- aggregate(ord.quant ~ day.order + year.order, orders.x.max, FUN=sum, drop = F)
# cut off days before product entered portfolio
id.nan <- which( !is.na(day.dem.x.max$ord.quant) )
day.dem.x.max <- day.dem.x.max[min(id.nan):max(id.nan), ]
# each day without entry has  demand of 0
day.dem.x.max$ord.quant[is.na(day.dem.x.max$ord.quant)] <- 0

# recompose date
tmp <- make_date(year = day.dem.x.max$year.order)
yday(tmp) <- as.numeric(day.dem.x.max$day.order)

# add to data frame
day.dem.x.max <- data.frame(day.dem.x.max,
                            date = tmp)

# make time series plot
ggplot(day.dem.x.max, aes(x=date, y=ord.quant)) + 
  geom_point() + geom_bar(stat='identity') +
  ylab("order quantity")

sq.inv <- function(d, s, q, l.ini, L = 1  ){
        # number of periods
        n <- length(d)
       # initialize result matrix
        tmp.mat <- matrix(0, nrow = n + 1 + L, ncol = 5 )
        colnames(tmp.mat) <- c("demand","inventory","disp.inv", "order", "delivery")
        # initialize matrix
        tmp.mat[,"demand"] <- c(0, d, rep(0,L))
        tmp.mat[1,"inventory"] <- l.ini
        tmp.mat[1,"disp.inv"] <- l.ini
        # inventory calculation
        for(i in 2:(n+1)){
            # check whether reorder level is reached
            if(tmp.mat[i-1, "disp.inv"] <= s){
                tmp.mat[i, "order"] <- q
                tmp.mat[i+L, "delivery"] <- q
            }
        # inventory balances  
        tmp.mat[i, "disp.inv"] <- tmp.mat[i-1, "disp.inv"] - tmp.mat[i, "demand"] + tmp.mat[i, "order"]
        tmp.mat[i, "inventory"] <- tmp.mat[i-1, "inventory"] - tmp.mat[i, "demand"] + tmp.mat[i, "delivery"]
        }
        # return results
        return(as.data.frame(tmp.mat[2:(n+1),]))
    }

eval.inv <- function(result, c.order, c.inv, c.stockout = 0){
  # determine periods with deliveries 
  tmp.id <- result[,"delivery"] > 0
  #  number of orders
  nb.orders <- sum(tmp.id)
  # total stock
  tot.stock <- result %>% filter(inventory > 0) %>% pull(inventory) %>% sum
  # total backlog
  tot.backlog <- abs(result %>% filter(inventory < 0) %>% pull(inventory) %>% sum)
  # total cost per period
  total.cost <- (nb.orders * c.order +  tot.stock * c.inv + tot.backlog * c.stockout)/nrow(result)
  # retain stock levels befor delivery arrives
  tmp.min.stock <- result[which(tmp.id) - 1,"inventory"]
  # alpha-service level (per period, not order cycle)
  alpha.sl <- 100 * ( 1 - sum( tmp.min.stock < 0 ) / nrow(result) )
  # Beta-service level
  beta.sl <- (100 * (1 - sum( -tmp.min.stock[ tmp.min.stock < 0]) / sum(result[,"demand"])))
  # return results
  return(list(  total.per.cost = total.cost, alpha.sl = alpha.sl, beta.sl = beta.sl,
    nb.orders = nb.orders, total.stock = tot.stock, total.backlog = tot.backlog
  ))}

For the order quantity, Harris-Andlers formula can be used: \[ q = \sqrt{\frac{2 \cdot \mu_d \cdot c_o}{c_i}}\] For \(s\) holds that it controlls the stock-out risk. To assess the risk, the demand in the lead time is the sum of \(L\) negative binomially distributed variables \(d_t \sim NB(r,p)\), for which holds \[D^L = \sum_{t=1}^L d_t \sim NB(L \cdot r, p)\] For controlling the \(\alpha\) service level, the \(\alpha\) quantile of \(D^L\) can be used. For the \(\beta\) service level, the corresponding loss function \(V_{NB}(x, r, p)\) can be used as follows

\[ s = V^{-1}_{NB}\left((1-\beta)\cdot \mu_{d} \cdot L ,r,p\right)\]

Most multi-purpose optimization routines in R are dedicated to continuous optimization as this is most often encoutered in statistics. Additionally, also constrained optimization problems appear which are notoriously difficult to solve in their most general form. A comprehensive overview can be found here.

For unconstraint (or at most box-constraint) general-purpose optimization, R offers the function optim(). The syntax of both functions is identical: optim(par = <initial parameter>, fn = <obj. function>, method = <opt. routine>). Additionally, upper and lower values for the parameters can be set by option lower = <vector of lower bounds> and upper = <vector of upper bounds>. In this case, method = "L-BFGS-B" must be selected.

Example: Continuous location planning

\[\min\limits_{(u,v) \in \mathbb R^2} \rightarrow \sum_{i \in \mathcal{C} } a_i \cdot (|x_i - u| + |y_i - v|)\]

loc.mh <- function(loc, a, x, y) sum(a * (abs(x - loc[1]) + abs(y - loc[2]) ) ) 
n <- 100
df.mh.loc <- data.frame(
  weight = sample(1:100, size = n),    # sample weights for each point/customer
  x = rnorm(n),                   # sample x coordinates
  y = rnorm(n))                   # sample y coordinates
res <- optim(par = c(0,0), fn = loc.mh, method = "BFGS", a = df.mh.loc$weight, x = df.mh.loc$x , y= df.mh.loc$y) 
res$par   
> [1] -0.1303003 -0.1399569

# generate grid of potential locations
emat <- expand.grid(x=seq(-3,3, length.out = 100), y=seq(-3,3, length.out = 100))
emat <- emat %>% add_column(objective = apply(emat, 1, function(x) 
  loc.mh(x, a = df.mh.loc$weight, x = df.mh.loc$x, y = df.mh.loc$y)))
# generate plot
ggplot(NULL, aes(x, y))  + geom_point(data = df.mh.loc , aes(size=weight)) + 
  geom_point(aes(x=res$par[1], y=res$par[2]), colour="blue", shape = 17, size= 3) + 
  geom_contour(data = emat, aes(z = objective))  +  geom_text_contour(data = emat, aes(z = objective), stroke=.15)

# wrapper function for simulating & evaluating (s,q) setting
sim.inv <- function(x, d.vec = day.dem.x.max[,"ord.quant"]){
  tmp <- sq.inv(d = d.vec, s = x[1],  q = x[2], l.ini = 10, L = 3)
eval.inv(tmp, c.inv=0.01, c.order=100, c.stockout = 0.5)$total.per.cost
} 

For finding the cost-optimal setting of \(s\) and \(q\), the following constraints are modelled: \[ \begin{align*} -q + s & \geq 0\\ q & > 0\\ s & \geq 0 \end{align*} \]

c.mat <- cbind( c(-1,1,0), c(1,0,1))
c.vec <- c(0,-100,0)
res <- constrOptim(theta = c(12,151), f = sim.inv, ui = c.mat, ci=c.vec, grad=NULL) 
head(res,2)
> $par
> [1]   8.000324 131.000082
> 
> $value
> [1] 1.727231

# generate grid of parameter constellations
emat <- expand.grid(s=seq(6,10, length.out = 50), q=seq(110,150, length.out = 50))
emat <- emat %>% add_column(objective = apply(emat, 1, function(x) sim.inv(x)))
# generate plot
ggplot(emat, aes(x=s, y=q))  + geom_contour(aes(z=objective))  +  geom_text_contour(aes(z=objective), stroke=.15) +
geom_point(aes(x=res$par[1], y=res$par[2]), colour="red", shape = 17, size= 3)

Mixed-integer linear optimization problems (MILP) are characterized by linear objective functions and constraints w.r.t. the decision variables. However, some or all decision variables are integer and/or binary variables. In general, the canonical form of a MILP can be written as

\[ \min\limits_{ \mathbf{x}, \mathbf{y}, \mathbf{z}} \rightarrow \mathbf{c}_x^T \cdot \mathbf{x} + \mathbf{c}_y^T \cdot \mathbf{y} + \mathbf{c}_z^T \cdot \mathbf{z}\] \[ \begin{align} \mathbf{A}_x \cdot \mathbf{x} + \mathbf{A}_y \cdot \mathbf{y} + \mathbf{A}_z \cdot \mathbf{z} & \leq \mathbf{b} \\ \mathbf{x} & \in \mathbb{R} \\ \mathbf{y} & \in \mathbb{Z} \\ \mathbf{z} & \in \{0,1\} \\ \end{align} \]

Hence, a MILP basically consists of five parts:

1. coefficient vector \(c = (c_x, c_y, c_z)\)
2. constraint matrix \(A = (A_x,A_y,A_z)\)
3. right-hand side vector \(\mathbf{b}\)
4. direction of the constraints
5. the domain declarations of the decision variables.

Although any optimization problem consists of objective function, variables, and constraints, there are numerous ways to formalize these components. The package ROI provides a unified framework for seting up and solving optimization problems. An extensive description of the package can be found here.

library(ROI)            # load package
ROI_available_solvers() # check available solvers
ROI_installed_solvers() # check installed solvers

Theussl et al. (2017)

To formulate an optimization problem, the function OP(objective, constraints, types, bounds). All components are generated by creator functions. The types of variables are defined by a character vector of length \(n\) consisting of the characters ‘C’, ‘I’ or ‘B’ indicating continuous, integer or binary variables, respectively. Corresponsindgly, upper and lower bounds of the variables can be given as a named list consisting of two vectors of length \(n\) and named “lower” and “upper”.

For the objective function these creator functions are defined as follows:

The second component of optimization problems are the constraints. In general, the constraints consist of the left-hand side which depends on the decision variables \(x\), the right-hand side vector (\(\mathbf{b}\), rhs), and the direction vector (dir), e.g. \(lhs(\mathbf{x})\leq \mathbf{b}\). The left-hand side is created in a similar way as the objective function using the following creators

\[min \rightarrow \sum_{i,j =1}^n x_{ij} \cdot c_{ij}\] \[\begin{align}\small \sum_{i=1}^n x_{ij} &= 1 & \forall j = 1,...,n\\ \sum_{j=1}^n x_{ij} &= 1 & \forall i = 1,...,n\\ x_{ij} & \in \left\{0,1\right\} \end{align}\]

For a linear assignment problem with 3 rooms follows \(\small \mathbf{x} = (x_{11},x_{12},...,x_{1n},x_{21},...,x_{2n},...,x_{nn})\) and \[ \small \mathbf A = \left[ \begin{array}{*{10}{r}} 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\ \end{array} \right] \]

# create problem dimension
n <- 10
# create cost vector
c.vec <- rpois(n^2, 100)    # sample costs from Poisson distribution with mean 100
# create empty constraint matrix
L.mat <- matrix(0, ncol = n^2, nrow = 2*n)
# first n rows
L.mat[1:n,] <- do.call(cbind, lapply(1:n, function(x) diag(n)))
# last n rows
L.mat[(n+1):(2*n),] <- t(sapply(1:n, function(x) c(rep(0,(x-1)*n), rep(1,n),rep(0,(n-x)*n) )))

# create optimization problem
copt <- OP(
  objective = L_objective(L = c.vec),
  constraints = L_constraint(L = L.mat, dir = rep("==", 2*n), rhs = rep(1, 2*n)),
  types = rep("B", n^2)
)

# solve the problem
copt_sol <- ROI_solve(copt)        

# display solution in matrix form
sol.mat <- matrix(copt_sol$solution, ncol = n)
# nodes
nodes <- data.frame(id = 1:20,  label = c(paste("room", 1:10), paste("element", 1:10)),
  group = rep(c("rooms", "elements"), each=10),  title = c(paste0("<p><b>room ", 1:10,"</b><br></p>"),
            paste0("<p><b>element ", 1:10,"</b><br></p>")),  color = rep(c("darkred","darkblue"), each=10))
# edges
edges <- data.frame(from = 1:10, to = 10+apply(sol.mat, 1, function(x) which(x>0)), label = matrix(c.vec, ncol = n, byrow=T)[which(t(sol.mat)>0)] )
# visNetwork(nodes, edges, height = "200px", width = "99.999%") %>%
#  visLayout(hierarchical = TRUE) 

To formulate large (mixed-)integer optimization problems. algebraic modelling languages are usually used. Commercial solvers often provide a tailored algebraic programming language like ILOG OPL. The package ompr provides a similar access for modelling MILPs in a algebraic way. The package uses the ROI framework for accessing MILP solvers.

Models in ompr are set up by adding pats of the model via piping with the following functions:

Note that all expressions must be vectorized, i.e., should be able to process vectors and return results as a vector.

Given a set of nodes \(\mathcal{N}\) with \(n = |\mathcal{N}|\) and pairwise distances \(d_{ij}\) between any pair of nodes (\(i,j \in \mathcal{N}\)), the TSP is to find the shortest tour which visits all nodes and returns to the origin. To avoid subtours the Miller-Tucker-Zemlin formulation is used:

\[\min \rightarrow \sum_{i=1}^n \sum_{j=1}^n d_{ij} \cdot x_{ij}\]

\[ \begin{align} \sum_{i \in \mathcal{N}: i\neq j} x_{ij} & = 1 & \forall j \in \mathcal{N} \\ \sum_{j \in \mathcal{N}: j\neq i} x_{ij} & = 1 & \forall i \in \mathcal{N} \\ u_j & \geq u_i + 1 - n \cdot (1-x_{ij}) & \forall i,j \in \mathcal{N}, i \neq j, i,j>1 \\ u_1 & = 1 & \\ u_i & \in \left[2,n\right] & \forall i \in \mathcal{N}, i>1\\ x_{ii} & = 0 & \forall i \in \mathcal{N} \\ x_{ij} & \in \left\{0,1\right\} & \forall i,j \in \mathcal{N} \end{align} \]

# set number of nodes  & sample node coordinates
n <- 10; cities <- data.frame(id = 1:n, x = runif(n), y = runif(n))
# calculate distances
distance <- as.matrix(dist(cities[,c("x","y")]))       
# create vectorized function for accessing distance values
dist_fun <- function(i, j)  vapply(seq_along(i), function(k) distance[i[k], j[k]], numeric(1L))
# remove variables with same name as your decision var.
rm(x); rm(u)
# create model
model <- MIPModel() %>%
  # create  variables x and u
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "integer", lb = 0, ub = 1) %>%
  add_variable(u[i], i = 1:n, lb = 1, ub = n) %>% 
  # minimize travel distance
  set_objective(sum_expr(dist_fun(i, j) * x[i, j], i = 1:n, j = 1:n), "min") %>%
  # leave each node
  add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = 1:n) %>%
  # visit each node
  add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 1:n) %>%
  # ensure no subtours (arc constraints)
  add_constraint(u[1] == 1) %>% 
  add_constraint(u[j] >= u[i] + 1 - n * (1 - x[i, j]), i = 1:n, j = 2:n) %>%
  # exclude self-cycles
  set_bounds(x[i, i], ub = 0, i = 1:n)

For delivering the customers in Brazil, 5 distribution centers (DC) shall be installed such that the total weighted distance between the customers and nearest DC is minimized. The customers’ weights are the number of orders in the dat set. The optimization problem can be set up as follows \[ \min \rightarrow \sum_{i \in \mathcal C} w_i \cdot \sum_{j \in \mathcal D} x_{ij} \cdot d_{ij}\] \[ \begin{align*} \sum_{j \in \mathcal D} x_{ij} & = 1 & \forall i \in \mathcal C\\ x_{ij} & \leq y_j & \forall i \in \mathcal C, j \in \mathcal D\\ \sum_{j \in \mathcal D} y_j & \leq K & \\ x_{ij}, y_{j} & \in \{0,1\} & \forall i \in \mathcal C, j \in \mathcal D \end{align*} \]

Customer locations are identified by the first 5 digits (out of 8) of the postal zip number.

# change name of zip code name for subsequent merging
colnames(customer.data)[3] <- "zip_code"
colnames(geo.data)[1:3] <- c("zip_code","lat","lng")
# calculate median coordinates of zip areas
geo.data.sub <- data.frame(
  aggregate(lat ~ zip_code, FUN=median, 
            data = geo.data), 
  lng = aggregate(lng ~ zip_code, FUN=median, 
                  data = geo.data)[,2])
# merge geolocation data to customer data 
cust.data.add <- merge(customer.data, geo.data.sub, 
                       by="zip_code")
# remove outlier
cust.data.add <- filter(cust.data.add, lng < -30)
# generate plot
data(World)
# obtain brazil map
bra_map <- get_brmap(geo = "State", class = "sf")
# convert customer coordinates
cust.points <- st_as_sf(cust.data.add, 
                        coords = c("lng", "lat"))
# plot brazil map
 tm_shape(bra_map) + tm_borders() +
   tm_shape(cust.points) + tm_dots()