##################################### # 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 ) ) } #' #'