| |||||||||||
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.
wagonloadbalancing.R ##################################### # This file is part of the # # Xpress-R interface examples # # # # (c) 2022-2024 Fair Isaac Corporation # ##################################### #' --- #' title: "Wagon Load Balancing" #' 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 ) # 2. binary variables 'load' # create a data frame to store the 'load' indices loads_index <- data.frame( 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 loads_index %>% group_by(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) # add the heuristic solution addmipsol( 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 ) ) } #' #' | |||||||||||
© Copyright 2024 Fair Isaac Corporation. |