An Integer Programming based solution to a multi-Knapsack problem posted in StackOverflow. This solution uses R and in particular the lpSolveAPI package. The question can be found at: http://stackoverflow.com/questions/18036966/algorithm-for-combining-numbers-into-optimal-groups-which-do-not-exceed-a-vector

multi_knapsack.R

library(lpSolveAPI)
library(reshape2)
library(ggplot2)


#Problem definition
trim.lengths <- c(44, 57, 86, 86, 40, 88, 88, 41, 40, 40, 62, 62, 54, 54, 55, 55, 63, 63, 75, 75, 100, 100)
avail.lengths <- c(120, 103, rep(100, 9), 98, rep(97, 4), 95, rep(88, 3), 85, 65)
num.trims = length(trim_lengths) #22
num.avail = length(avail.lengths) #22
a.set<-1:num.avail
t.set <- 1:num.trims

#set rownames & colnames
setRowAndColNames<- function() {  
  a.and.t<- t(outer(a.set, t.set, paste, sep="_"))
  col.names <- paste("x_", a.and.t, sep="")
  s.vars <- paste("Surplus_", a.set, sep="")
  c.names <- c(s.vars, col.names)
  r.names <- setRowNames() 
  return( list(r.names, c.names)  )
}


setRowNames <- function() {
  cover.rows<- paste("CoverTrim", t.set, sep="_") #Cover constraints
  surplus.rows <- paste("DefineSurplus", a.set, sep="_") #Cover constraints
  r.names <- c(cover.rows, surplus.rows)
  return(r.names)
}


######## define constraint sets (called by populate.Amatrix)
populate.Amatrix <- function() {  
  df <- data.frame(matrix(0, n.rows, n.cols))  #create a dataframe shell  
  for (row in 1:num.trims) {
    skip = num.avail #for the surplus variables
    cols <- seq(0, (num.trims*num.avail)-1, by=num.avail)    
    df[row, skip+cols+row] <- 1
  }
  print("Done with cover Trims Constraints")
  st.row <- num.trims+1
  end.row<- st.row+num.avail-1
  for (row in st.row:end.row) {
    surplus.column <- row - num.trims
    df[row,surplus.column] <- 1
    current.a <- surplus.column
    acols <- num.avail + ((current.a-1)*num.trims) + (1:num.trims)
    df[row,acols] <- trim.lengths
    
  }
  return(df)
}


defineIP <- function(lp) {
  obj.vector <- c(rep(1, num.avail), rep(0, num.avail*num.trims))
  set.objfn(lp, obj.vector ) #minimize surplusses
  set.constr.type(lp, rep("=", n.rows), 1:n.rows) 
  rhs.vector <- c(rep(1, num.avail), avail.lengths)
  set.rhs(lp, b=rhs.vector, constraints=1:n.rows) # assign rhs values   
  
  #Set all the columns at once
  for (col in 1:n.cols) {
    set.column(lp, col, df[ ,col])  
  }
  
  for (col in (num.avail+1):n.cols) {
    print(col)
    set.type(lp, col, "binary")      
  }
  
  #assemble it all
  dimnames(lp) <- setRowAndColNames()
  write.lp(lp, "trimIP.lp", "lp")#write it out
}


# actual problem definition
n.rows <- num.trims + num.avail
n.cols <- (num.avail+1) * num.trims  #the extra one is for surplus variables
df <- populate.Amatrix()

lptrim <- make.lp(nrow=n.rows, ncol=n.cols)
defineIP(lptrim)
lptrim
solve(lptrim)
sol <- get.variables(lptrim)
sol <- c(sol, trim.lengths)
sol.df <- data.frame(matrix(sol, nrow=24, ncol=22 , byrow=T))
solution<- sol.df[-1, ]
solution$availlength <- c(avail.lengths, NA)
solution$surplus <- c(t(sol.df[1, ]), NA)
row.names(solution) <- c(paste("A", 1:22, sep=""), "trimlength")
solution$piece <- row.names(solution)
solution$first <- 0
solution$second <- 0

for (r in 1:22) {
  for (c in 1:22) {   
      if(solution[r,c]==1) { 
        if(solution$first[r] == 0) {
           solution$first[r] <- solution[23,c] 
          }
        else{
          solution$second[r] <- solution[23,c]           
        }
    }
    
  }
}

small.df <- solution[1:22, c(25, 24,26,27)] #just take the columns we want
df <- melt(small.df)
p <- ggplot(df, aes(x=piece)) + geom_bar(aes(y=value, fill=variable), stat="identity")
p #plot it.