FICO Xpress Optimization Examples Repository
 FICO Optimization Community FICO Xpress Optimization Home

Formulate and solve a Facility Location Problem

Description
Formulate and solve a Facility Location Problem

Further explanation of this example: Xpress R Reference Manual

Source Files
By clicking on a file name, a preview is opened at the bottom of this page.

facility_location.R

#####################################
# This file is part of the          #
# Xpress-R interface examples       #
#                                   #
#   (c) 2021 Fair Isaac Corporation #
#####################################
#' ---
#' title: "Facility Location Problem"
#' author: Gregor Hendel
#' date: Dec. 2020
#' ---
#'

#'
suppressMessages(library(xpress))

# for plotting the data
library(ggplot2)
library(magrittr)
library(reshape2)
library(dplyr)
library(RColorBrewer)

#'
#' Facility Location is an optimization problem which aims at connecting customers
#' to facilities at minimum costs. Costs in this optimization problem comprise an
#' initial cost for operating a facility, and additionally connection costs between
#' facilities and customers.
#'
#' The Facility Location Problem (FLP) is a classical application of mixed-integer
#' optimization and comes in many variants. We consider the capacitated facility
#' location problem, in which facilities can only serve customer demands up to a
#' given maximum capacity.
#'
#' This example is the baseline for many of the quick examples, during which we
#' will read the base model as an LP file and impose various additional
#' constraints.
#'
#' # Mathematical Formulation
#'
#' Let's review the mathematical formulation of Capacitated Facility Location. Let
#' $m,n \in \mathbb{N}$ and denote a set of facilities by $F = \{F_1,\dots,F_m\}$
#' and a set of customers by $C = \{C_1,\dots,C_n\}$. Furthermore, we have:
#'
#' - demands $d_j \geq 0$ for each customer $C_j$.
#' - capacities $u_i \geq 0$ for each facility $F_i$.
#' - opening costs $r_i \geq 0$ for each facility $F_i$,
#' - connection costs $s_{ij} \geq 0$ for connecting facility $F_i$ with customer $C_j$.
#'
#' With these parameters, we are ready to formulate the FLP:
#'
#' #' \begin{align} #' &\min & \sum\limits_{i = 1}^m r_i x_i + \sum\limits_{i=1,j=1}^n s_{ij}y_{ij}\\ #' & & \sum\limits_{i = 1}^n y_{ij} & = d_j & j = 1,\dots,n\\ #' & & \sum\limits_{j = 1}^n y_{ij} & \leq u_i x_i & i = 1,\dots,m\\ #' & & y_{ij} & \geq 0 & i=1,\dots m, \ j = 1,\dots,n\\ #' & & x_i & \in \{0,1\} & i=1,\dots m\\ #' \end{align} #'
#' We use binary variables $x_i$ to model the yes/no decision
#' whether to open facility $F_i$. Nonnegative continuous variables $y_{ij}$ are
#' used to model how much service customer $C_j$ receives from facility $F_i$. The
#' objective is formulated to incorporate the opening and service costs together.
#' The constraints ensure that the demand of each customer is satisfied, and that
#' no capacity is exceeded.
#'
#' Next, we turn this into an R model from standard data frames and solve it with
#' the FICO Xpress Optimizer.
#'
#' # Problem Data Preparation
#'
## ----Data Preparation---------------------------------------------------------

# six customers in a data frame
customer_df <- data.frame(
Demand=c(2,6,4,10,7,4)
)

# 5 facilities
facility_df <- data.frame(
Capacity=c(14,11,6,16,15),
Opening=c(3000, 4000,3500,8000,5000)
)

# some useful shorthand notation
m <- nrow(facility_df)
n <- nrow(customer_df)
mTimesN <- m * n

set.seed(811) # for reproducibility of the below random number generation
connection_df <- data.frame(
Facility=rep(1:m, each=n),
Customer=rep(1:n, m),
Cost=runif(mTimesN, min=1000, max = 2000)
)

#'
#' We visualize the generated data.
#'
## ----Visualizing the Generated Data, echo=FALSE-------------------------------
colors <- RColorBrewer::brewer.pal(10, "RdYlBu")[c(2,10)]
connection_df %>%
ggplot(aes(Customer, Facility, fill=Cost, label=format(Cost))) +
geom_tile() +
scale_fill_gradient(low = colors[2], high = colors[1]) +
geom_text() + ggtitle("Heat Map of Connection Cost")

#'
#' # The R Problem Data
#'
#' In order to load the model into Xpress at once, we fill a list object problemdata which
#' we pass into xprs_loadproblemdata.
#'
#' Each row of the facility_df data frame corresponds to one $x$ variable.
#' Each row of the connection_df data frame corresponds to one $y$ variable.
#'
## ----Problem Data, collapse=TRUE----------------------------------------------
# create an empty model as list object
problemdata <- list()

# Objective: opening costs for the x variables, connection costs for the y_ij variables
problemdata$objcoef <- c(facility_df$Opening, connection_df$Cost) # global entity types: binary for the x variables # # an alternative is to specify all column types # problemdata$columntypes <- c(rep('B', m), rep('C', mTimesN))
problemdata$coltype <- c(rep('B', m)) # the coltype approach requires us to specify the indices # of the binary variables explicitly. We have to use 0-based # indexing for Xpress! problemdata$entind <- 0:(m-1)

# demand constraints: 1 constraint for each customer
demand_matrix_transposed <- sapply(1:n, function(x) c(
rep(0,m),
connection_df$Customer == x)) # capacity constraint on each facility capacity_matrix_transposed <- sapply(1:m, function(x) c(ifelse( (1:m)==x, -facility_df$Capacity[x], 0
),
connection_df$Facility==x) ) # we transpose and glue both matrices together into the final A problemdata$A <- rbind(
t(demand_matrix_transposed),
t(capacity_matrix_transposed)
)

# right hand side declaration: customer demand for customer or 0
problemdata$rhs <- c( customer_df$Demand,
rep(0, m)
)

# lower and upper bounds on the variables
problemdata$lb=rep(0,m+mTimesN) problemdata$ub=c(rep(1,m), rep(1e+20, mTimesN))

# sense of the rows:
#'E'quations,
#'L'ess/'G'reater than (<=/>=),
#'R'ange rows, 'N'onbinding (neutral) rows
problemdata$rowtype <- c( rep("E", n), # n equalities to satisfy demands rep("L", m) # m inequalities not to exceed capacity ) # declare column names problemdata$colname <- c(
sprintf("x_%d", 1:m),
sprintf("y_%d_%d", connection_df$Facility, connection_df$Customer)
)

# declare row names
problemdata$rowname <- c( sprintf("Demand_%d", 1:n), sprintf("Capacity_%d", 1:m) ) #' #' Let's have a graphical look at the matrix to verify that we made no errors. For #' a better display, we have flipped columns onto the vertical axis. It appears #' that we translated the matrix correctly. Remember that variable$y_{35}$models #' the amount of service between facility$F_3$and customer$C_5$As expected #'$y_{35}$appears in 2 constraints with a coefficient of 1, namely the capacity #' constraint of$F_3$, and the demand constraint of$C_5$. #' ## ----echo=FALSE--------------------------------------------------------------- a_slim <- reshape2::melt(problemdata$A)
a_slim$rowname <- problemdata$rowname[a_slim$Var1] a_slim$colname <- problemdata$colname[a_slim$Var2]
a_slim %>% dplyr::filter(value != 0) %>%
ggplot(aes(rowname, colname, fill=value, label=value)) +
geom_raster() + geom_text() +
theme_bw() +
theme(axis.text.x = element_text(size = 5, angle = 30),
axis.text.y = element_text(size = 6))

#'
#' # Solving the Model with FICO Xpress
#'
#' It is now very easy to input and solve this model using FICO Xpress.
#'
## ----Solving The Problem------------------------------------------------------
problemdata$probname <- "FacilityLocation" # load the problem into a newly created problem p <- xprs_loadproblemdata(problemdata=problemdata) # enable output to console setoutput(p) # save the problem in an LP file writeprob(p, "flp.lp", "l") # call optimize summary(xprs_optimize(p)) #' ## ----Displaying the Solution-------------------------------------------------- # get solution and split it into facility and connection sol <- xprs_getsolution(p) facility_sol <- sol[1:m] connection_sol <- sol[m+(1:mTimesN)] # print the open facilities: print("Open facilities:") print(sprintf("F_%d", which(facility_sol > 0))) # print the connections: print("Connections:") connection_df$TotalCost <- connection_df\$Cost * connection_sol

# use a slight numeric tolerance to filter the actual connections
print(connection_df[connection_sol > 1e-9,])