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

Solve a modified bin packing problem with the addition of a simple greedy heuristic

Description

In this example, 16 boxes are required to be loaded into 3 railway wagons each with capacity of 100 quintals. This example is similar to the bin packing problem, but instead of minimizing the number of bins used, here we would like to minimize the heaviest wagon load. This example also shows how to add a user solution found by a simple greedy heuristic.

Further explanation of this example: This is a conversion of the Mosel example 'Wagon Load Balancing'

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

#####################################
# This file is part of the          #
# Xpress-R interface examples       #
#                                   #
#   (c) 2021 Fair Isaac Corporation #
#####################################
#' ---
#' author: Y. Gu
#' date: Jun.2021
#' ---
#'
#'
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(results = "hold")
knitr::opts_chunk$set(warning = FALSE,message = FALSE) #' #' #' ## Brief Introduction To The Problem #' #' [This is a conversion of the Mosel example 'Wagon Load Balancing'](https://www.fico.com/fico-xpress-optimization/docs/latest/examples/mosel/ApplBook/D_LoadCut/d1wagon2.mos). #' A brief introduction to this problem is given below, and to see the full #' mathematical modeling of this problem you may refer to section 9.1, page 122 of the #' book 'Applications of optimization with Xpress'. #' #' In this example, 16 boxes are required to be loaded into 3 railway wagons each with #' capacity of 100 quintals. This example is similar to the bin packing problem, but #' instead of minimizing the number of bins used, here we would like to minimize the #' heaviest wagon load. #' #' We #' formulate this example as a MIP problem and solve it using Xpress optimizer. #' #' The objective of the MIP problem is the maximum weight over all the wagon loads which #' we want to minimize. We define it as an integer variable 'maxweight'. Since the ideal #' case is that all wagons take the same load and the box weights are all integers, so #' we round the average of total box weights to the next larger integer and set this #' number as the lower bound of 'maxweight'. Without this lower bound, the optimum #' solution will be found easily but it will take a long time to prove its optimality. #' We also define binary variables for all possible pairs of boxes and wagons, and value #' '1' indicates that the box is assigned to the wagon and value '0' means no assignment. #' #' The constraints are clear, firstly, boxes are not fragmented and thus each box can #' only be assigned to one wagon. Then, we set 'maxweight' as upper bound of the weight #' of the boxes loaded into each wagon. By proceeding this way, in the optimal solution #' the minimization will force 'maxweight' to take the value corresponding to the weight #' of the heaviest wagon load. The mathematical formulations of these constraints are #' included in the guide book 'Applications of optimization with Xpress'. #' #' Before we solve the MIP, we create an initial solution using the following #' simple heuristic: #' In each iteration, we greedily choose the heaviest unloaded box and assign it #' to the wagon with the least load until all boxes are assigned. #' This does not yield a feasible solution for the data of the example, but can #' still be added to the MIP as an infeasible starting point to be repaired. #' #' #' For this example, we need packages 'xpress' and 'dplyr'. Besides, we use the #' function pretty_name to give the variables and constraints concise names. #' ## ----Load The Packages And The Function To Give Names------------------------- library(xpress) library(dplyr) pretty_name <- function(prefix, y) { "%s_%s" %>% sprintf(prefix, paste(lapply(names(y), function(name) paste(name, y[name], sep = "_")), collapse = "_")) } #' #' #' Add the values we need for this example. #' ## ----Data--------------------------------------------------------------------- # concrete values weight.df <- data.frame( Box = 1:16, Weight = c(34, 6, 8, 17, 16, 5, 13, 21, 25, 31, 14, 13, 33, 9, 25, 25) ) # carrying capacity of each wagon WMAX <- 100 # wagons Wagon <- 1:3 #' #' #' Create a function to get the heuristic solution. #' ## ----Heuristic---------------------------------------------------------------- solve_heur <- function(weight, wagon) { # sort the weight in descending order sorted <- sort( weight, decreasing = TRUE, method = "shell", index.return = TRUE ) weight.sorted <- sorted[[1]] box.sorted <- sorted[[2]] n = length(wagon) CurWeight <- rep(0, n) # current weight of wagon loads CurNumber <- rep(0, n) # current number of boxes per wagon Load <- vector("list", n) # boxes loaded onto the wagons l = length(weight.sorted) for (i in 1:l) { v = which.min(CurWeight) CurNumber[v] <- CurNumber[v] + 1 Load[[v]] = c(Load[[v]], box.sorted[i]) CurWeight[v] <- CurWeight[v] + weight.sorted[i] } results <- list(max(CurWeight), Load) names(results) <- c("maximum_weight", "Load") return(results) } #' #' #' Create a function to formulate this example as a MIP problem and solve it. #' ## ----Define And Solve the MIP Problem----------------------------------------- solve_MIP <- function(weight.df,Wagon){ # firstly, create a new problem prob <- createprob() # set the problem name setprobname(prob, "WagonLoadBalancing") # add columns # 1. variable 'maxweight' maxweight <- xprs_newcol( prob, lb = ceiling(sum(weight.df$Weight) / 3),
ub = Inf,
coltype = "I",
name = "maxweight",
objcoef = 1
)

# create a data frame to store the 'load' indices
idx = 1:(length(Wagon) * nrow(weight.df)),
Wagon = rep(Wagon, each = nrow(weight.df)),
Box = rep(weight.df$Box, length(Wagon)), Weight = rep(weight.df$Weight, length(Wagon))
)

loads_index$loads <- loads_index %>% apply(1, function(x) xprs_newcol( prob, lb = 0, ub = 1, coltype = "B", name = pretty_name("Load_", x["idx"]) )) # add rows # 1. each box can only be assigned to one wagon apply(weight.df, 1, function(x) xprs_newrow( prob, colind = (loads_index %>% filter(Box == x["Box"]))$loads,
rowcoef = rep(1, length(Wagon)),
rowtype = "E",
rhs = 1,
name = pretty_name("Box_Wagon", x["Box"])
))

# 2. set 'maxweight' as the upper bound of weight loaded into each wagon
group_map( ~ xprs_newrow(
prob,
colind = c(.x$loads, maxweight), rowcoef = c(.x$Weight, -1),
rowtype = 'L',
rhs = 0,
name = pretty_name("Wagon_Capacity", .y)
))

# solve the heuristic and get the solution
heuristic <- solve_heur(weight.df$Weight, Wagon) HeurMaxWeight <- heuristic$maximum_weight # the max weight of the heuristic solution
HeurLoad <- heuristic$Load # list of boxes loaded to each wagon # store the heuristic solution in a data frame Heur.df <- cbind(data.frame(c( rep(1, length(HeurLoad[[1]])), rep(2, length(HeurLoad[[2]])), rep(3, length(HeurLoad[[3]])) )), do.call(rbind, lapply(HeurLoad, data.frame))) names(Heur.df) = c("Wagon", "Box") # represent the heuristic solution in binary form loads_index$HeurSol <- rep(0, nrow(loads_index))
for (i in 1:nrow(Heur.df)) {
loads_index$HeurSol[which(loads_index$Wagon == Heur.df$Wagon[i] & loads_index$Box == Heur.df$Box[i])] <- 1 } # we create the heuristic solution in this order since we first add variable 'maxweight', # and then the 'loads' variables Heuristic <- c(HeurMaxWeight, loads_index$HeurSol)

prob,
solval = Heuristic,
colind = c(maxweight, loads_index$loads), name = "Heuristic_solution" ) # solve the MIP problem setoutput(prob) summary(xprs_optimize(prob)) solution <- xprs_getsolution(prob) # optimum objective value objsolval <- getdblattrib(prob, xpress:::MIPOBJVAL) # load solutions: # Note that the column indices are 0-based while the vector indices in R are 1-based, # so, to get the corresponding solutions, we add '1' to column indices loads_index$MIPSol <- solution[loads_index$loads + 1] wagon_loaded <- loads_index %>% filter(MIPSol == 1) %>% select(Wagon, Box, Weight) # the total weight each wagon loads: sum_weight <- c() for (i in Wagon) { sum_weight <- c(sum_weight, sum(wagon_loaded$Weight[which(wagon_loaded$Wagon == Wagon[i])])) } weight_loaded <- data.frame(wagon = Wagon, weight = sum_weight) results <- list( prob = prob, Heuristic_max_weight = HeurMaxWeight, Heuristic_solution = Heur.df, MIP_maximum_weight = objsolval, MIP_wagon_loaded = wagon_loaded, MIP_weight_loaded = weight_loaded ) return(results) } #' #' #' Now we can solve this problem. We formulate the problem as a MIP problem, #' add a heuristic solution to the MIP formulation and then solve it. #' ## ----Solve The Problem-------------------------------------------------------- MIP_solution <- solve_MIP(weight.df, Wagon) print(MIP_solution) if (MIP_solution$Heuristic_max_weight <= WMAX) {
print(
paste(
"Heuristic solution is feasible with weight ",
MIP_solution$Heuristic_max_weight ) ) } else { print( paste( "Heuristic solution weight exceeds maximum allowed weight of ", WMAX, " by ", MIP_solution$Heuristic_max_weight - WMAX
)
)
}

#'
#'