Need to create $3-4$ different box sizes and to minimize material waste for a set of $n$ objects that need to fit into these boxes Excuse the non-mathematical way I've phrased the question.
I have the following problem:


*

*I have $N$ square paper documents with side lengths between $150$mm and $860$mm. I know each document side's length.

*I need to create $3-4$ differently sized boxes to fit all the documents, e.g. Three box types: Box $1$ side $l_1=300$mm, Box $2$ side $l_2=600$mm, Box $3$: $l_3=860$mm.

*There are as many boxes as documents, i.e. each document goes into its own separate box (of the smallest possible size so as to minimize waste of cardboard).

*What is the best way to decide on the size of the boxes, so as to minimize the total amount of (surface area) of cardboard used?
I'm not necessarily looking for the analytical solution to this problem. 
Two ideas I've had:
a) Pick $l_1$ and $l_2$ values at random and calculate the total surface area of cardboard. Guess the values again, see if the total surface area is smaller, and so on and on. 
b) A more analytical approach where I'm computing $l_1$ and $l_2$ value in say $1$mm increments and I calculate the total surface area for each combination of box lengths between say ($150$mm, $151$nmm,$860$mm) and ($858$mm,$859$mm,$860$mm). 
What would you suggest is the most practical way of going about solving this?
BTW, I'm great with Excel, less so with Mathlab, etc. I can program well in Ruby if that helps in any way. 
 A: As Erwin points out in his blog post, you can model this as a network. I would take that approach, in part because it requires no specialized software. Per Erwin's post, you have 384 distinct paper sizes. Create one node for each, and let $s_i$ be the paper size for node $i$ and $n_i$ the count for size $i$. For each pair of nodes $i < j$, draw an arc from node $i$ to node $j$ whose cost is $s_j^2 \sum_{k=i+1}^j n_k$. This arc represents the cost (surface area) of putting all pages with sizes between $s_{i+1}$ and $s_j$ into boxes of size $s_j$.
You can now iterate over the graph using either two or three nested loops (since you limited yourself to three or four box sizes). Start at node 1 and look at each possible successor node (outer loop), each possible successor to that node (inner loop), each possible successor to that node (nested loop if you are allowing four sizes), recognizing that you must take the arc from the node in the innermost loop to node 384. You sum the lengths of the selected arcs, then compare the sum to the best solution so far. If it's shorter, update the best solution. Finally, note that you can break off any inner loop if its cumulative sum equals or exceeds the best some so far, since adding more (positive) arc costs cannot reduce the sum.
I would not advocate brute force in general, but with a maximum of four box sizes, and given the speed of a contemporary PC, this should be rather doable (and, again, requires no special software, other than a compiler/interpreter for some programming language).
A: A Dynamic Programming algorithm in R can look like:
data <- c(
 156,162,168,178,178,180,185,185,190,192,193,194,195,195,
 . . .
 805,806,820,823,827,827,855,855,864)

# get unique values and counts
t<-table(data)
count <- as.vector(t)
size <- as.numeric(rownames(t))

cumulative <- cumsum(count)

# number of box sizes, number of item sizes
NB <- 4
NI <- length(size)

# allocate matrix NI rows, NB columns (initialize with NAs)
# f[ni,nb] = cost when we have ni items and nb blocks
F <- matrix(NA,NI,NB)
S <- matrix("",NI,NB)

# initialize for nb=1
F[1:NI,1] <- cumulative * size[1:NI]^2
S[1:NI,1] <- paste("(",1,"-",1:NI,")",sep="")

# dyn programming loop
for (nb in 2:NB) {
  for (ni in nb:NI) {
     k <- (nb-1):(ni-1)
     v <- F[k,nb-1] + (cumulative[ni]-cumulative[k])*size[ni]^2
     F[ni,nb] <- min(v)

     # create path (string)
     mink <- which.min(v) + nb - 2
     s <- paste("(",mink+1,"-",ni,")",sep="")
     S[ni,nb] <- paste(S[mink,nb-1],s,sep=",")
  }
}

for (nb in 1:NB) {
   cat(sprintf("%s boxes: %s, total area = %g\n",nb,S[NI,nb],F[NI,nb]))
}

The solution looks like:
1 boxes: (1-384), total area = 8.70414e+08
2 boxes: (1-268),(269-384), total area = 4.59274e+08
3 boxes: (1-155),(156-284),(285-384), total area = 3.66837e+08
4 boxes: (1-155),(156-268),(269-351),(352-384), total area = 3.29199e+08

We found the same solution as with the other methods.
