Monday, October 11, 2021

2d knapsack problem

This post is about what sometimes is called a 2d knapsack problem. It can also be considered as a bin-packing problem. These are famously difficult problems, so let's have a look at a relatively small instance. 

So, we have one 2d container of a given size. Furthermore, we have a collection of items: rectangles of a given size, with a value. The goal is to pack items in our container without overlapping each other such that the total value is maximized.

To make things concrete, let's have a look at a small data set (randomly generated):


----     18 PARAMETER data  problem data

                width      height   available       value  value/size

k1             20.000       4.000       2.000     338.984       4.237
k2             12.000      17.000       6.000     849.246       4.163
k3             20.000      12.000       2.000     524.022       2.183
k4             16.000       7.000       9.000     263.303       2.351
k5              3.000       6.000       3.000     113.436       6.302
k6             13.000       5.000       3.000     551.072       8.478
k7              4.000       7.000       6.000      86.166       3.077
k8              6.000      18.000       8.000     755.094       6.992
k9             14.000       2.000       7.000     223.516       7.983
k10             9.000      11.000       5.000     369.560       3.733
container      30.000      20.000


Our container is \(30\times 20\). We have 10 different types of items. The column available indicates the number of items of the current type. This means we have actually 51 items. Instead of value, it may be better to look at the value/size column. This is the relative value, normalized for size (i.e. height times width). The four most valuable items are k6, k9, k8, and k5. Hopefully, we will see them back in the solution.

An interesting wrinkle is: the items can be rotated. There are different ways of handling this. We can duplicate all items and add a rotated version. However, we need to make sure that the availability limit is applied to both. In the first MIP model below, I decided to add a rotation index \(r\) to indicate whether an item is rotated. This is almost the same as duplicating items. In the second model, I introduce an explicit binary variable to indicate whether the item is rotated.

Our solution will consist of a subset of the items plus a 2d no-overlapping configuration of those items. There is no way to just ask for the optimal selection of items: we always need to figure out the actual layout.

In this post, I'll discuss three different models:

  • A covering MIP model that works with integer widths, heights, and locations,
  • a continuous size MIP model where the sizes, and locations, can be fractional,
  • and a CP-SAT model where all coefficients need to be integer-valued.

These three models are very different from each other. They share almost nothing!


A covering model


For this model, we assume that the container is a grid of cells, and each item occupies a number of cells. In other words, I assume the width and height are integer-valued. The main decision variable is: \[\color{darkred}x_{k,r,i,j} = \begin{cases} 1 & \text{if an item of type $k$ and direction $r$ is placed at cell \((i,j)\)} \\ 0 & \text{otherwise}\end{cases}\] The index \(r\) has two elements: \({r,nr}\) indicating if the corresponding item \(k\) is rotated or not. My convention is: if \(\color{darkred}x_{k,r,i,j}=1\) for a non-rotated item, then cells \((i',j')\), with \(i'=i,..,i+\color{darkblue}{\mathit{width}}_k-1\) and \(j'=j,..,j+\color{darkblue}{\mathit{height}}_k-1\) are occupied by this item. We need some rule to precisely state what it means if an item \(k\) is placed at \((i,j)\).

The first thing I did was creating two sets: \(\color{darkblue}{\mathit{ok}}\) and \(\color{darkblue}{\mathit{cover}}\).

The first set is \(\color{darkblue}{\mathit{ok}}_{k,r,i,j}\) indicating the allowed cells \((i,j)\) where we can place an item of type \(k\) (rotated or not). Below we see the allowed assignments of item k3:


ok(k,r,i,j)

The second set is \(\color{darkblue}{\mathit{cover}}_{k,r,i,j,i',j'}\) indicating which cells \((i',j')\) are covered when placing item \(k\) on cell \((i,j)\). This is a large 6 dimensional set. It has 282,944 elements.


cover(k,r,i,j,ii,jj)


Here we see the coverage when we place item k1 at (r4,c6) (not rotated).


With this we can write our covering model:


Covering Model
\[\begin{align}\max&\>\color{darkred}{\mathit{totalValue}} = \sum_{k,r,i,j|\color{darkblue}{\mathit{ok}}(k,r,i,j)} \color{darkblue}{\mathit{value}}_k \cdot \color{darkred}x_{k,r,i,j} && &&\text{Maximize total value}\\ & \sum_{k,r,i,j|\color{darkblue}{\mathit{cover}}(k,r,i,j,i',j')}\color{darkred}x_{k,r,i,j} \le 1 && \forall i',j' && \text{No overlap} \\ & \sum_{r,i,j|\color{darkblue}{\mathit{ok}}(k,r,i,j)} \color{darkred}x_{k,r,i,j} \le \color{darkblue}{\mathit{avail}}_k &&\forall k && \text{Availability} \\ & \color{darkred}x_{k,r,i,j} \in \{0,1\} \end{align}\]


I considered adding the constraint: \[\sum_{k,r,i,j|\color{darkblue}{\mathit{ok}}(k,r,i,j)} \color{darkblue}{\mathit{area}}_k \cdot \color{darkred}x_{k,r,i,j} \le \color{darkblue}{\mathit{containerSize}}\] This did not seem to be much of a help.

The model is not very small:


MODEL STATISTICS

BLOCKS OF EQUATIONS           3     SINGLE EQUATIONS          611
BLOCKS OF VARIABLES           2     SINGLE VARIABLES        4,273
NON ZERO ELEMENTS       291,489     DISCRETE VARIABLES      4,272

However, this model solves quickly: just 63 nodes were needed. The results are:
 
----     83 VARIABLE x.L  item k is placed at (i,j)

             r1.c1      r1.c14      r1.c28      r3.c14      r5.c14       r6.c1       r6.c6      r7.c11      r7.c29

k5 .nr                                   1
k6 .nr           1
k6 .r                                                                        1           1
k8 .r                                                                                                1
k9 .nr                       1                       1           1
k9 .r                                                                                                            1

     +     r13.c11      r19.c1     r19.c15

k8 .r            1
k9 .nr                       1           1


----     83 VARIABLE totalValue.L          =     4617.938  objective (maximized)


This is a bit difficult to interpret. A plot can help here:

 

Notice that all four valuable items (k6, k9, k8 and k5) are present in the solution!

This picture makes it look like a very small problem. But remember that we had to pick rectangles from a 51-item collection. Because of the possibility of rotation, it is even better described as picking from a 102-item inventory.

A continuous size model


Instead of using a grid, here we use continuous widths and heights. For this, we use pairwise no-overlap constraints.  Let's show the idea behind this.

Intermezzo: pairwise no-overlap constraints. These constraints are: \[\begin{align} \color{darkred}x_a  &\ge \color{darkred}x_b + \color{darkblue}w_b \\ {\bf or}\>\color{darkred}x_b &\ge \color{darkred}x_a + \color{darkblue}w_a \\ {\bf or}\>\color{darkred}y_a &\ge \color{darkred}y_b + \color{darkblue}h_b \\ {\bf or}\>\color{darkred}y_b &\ge \color{darkred}y_a + \color{darkblue}h_a \end{align}\] The or conditions can be implemented with additional binary variable and big-M constraints: \[\begin{align} & \color{darkred}x_a  \ge \color{darkred}x_b + \color{darkblue}w_b  - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,1})\\ & \color{darkred}x_b \ge \color{darkred}x_a + \color{darkblue}w_a - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,2})\\ &\color{darkred}y_a \ge \color{darkred}y_b + \color{darkblue}h_b - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,3})\\ & \color{darkred}y_b \ge \color{darkred}y_a + \color{darkblue}h_a - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,4}) \\ & \sum_k \color{darkred}\delta_{a,b,k} \ge 1 \end{align}\] We need to add one more detail: we want to have these constraints only active if both items a and b are selected to be in our container. So our constraints become:  \[\begin{align} & \color{darkred}x_a  \ge \color{darkred}x_b + \color{darkblue}w_b  - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,1}) - \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_a) - \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_b)\\ & \color{darkred}x_b \ge \color{darkred}x_a + \color{darkblue}w_a - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,2}) - \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_a) - \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_b)\\ &\color{darkred}y_a \ge \color{darkred}y_b + \color{darkblue}h_b - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,3})- \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_a) - \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_b) \\ & \color{darkred}y_b \ge \color{darkred}y_a + \color{darkblue}h_a - \color{darkblue}M \cdot (1-\color{darkred} \delta_{a,b,4}) - \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_a) - \color{darkblue}M\cdot (1-\color{darkred}{\mathit{selected}}_b)\\ & \sum_k \color{darkred}\delta_{a,b,k} \ge 1 \end{align}\] The final complication is that we can rotate items. So the width can become the hight and vice versa.


We need to form pairs of items we may need to compare. I started with enumerating all possible items:

----     56 SET kn  (k,n) combinations

             n1          n2          n3          n4          n5          n6          n7          n8          n9

k1          YES         YES
k2          YES         YES         YES         YES         YES         YES
k3          YES         YES
k4          YES         YES         YES         YES         YES         YES         YES         YES         YES
k5          YES         YES         YES
k6          YES         YES         YES
k7          YES         YES         YES         YES         YES         YES
k8          YES         YES         YES         YES         YES         YES         YES         YES
k9          YES         YES         YES         YES         YES         YES         YES
k10         YES         YES         YES         YES         YES


----     56 PARAMETER num  sequence numbers for (k,n)

             n1          n2          n3          n4          n5          n6          n7          n8          n9

k1            1           2
k2            3           4           5           6           7           8
k3            9          10
k4           11          12          13          14          15          16          17          18          19
k5           20          21          22
k6           23          24          25
k7           26          27          28          29          30          31
k8           32          33          34          35          36          37          38          39
k9           40          41          42          43          44          45          46
k10          47          48          49          50          51
  

We see we have 51 items to consider. The sequence number can be used to form pairs we need to compare. The big trick is to compare item a with item b only once. I created a set for this:

Items to compare

An important advantage of creating these sets and parameters is that our model equations become much simpler. Also, this allows us to develop and debug a lot of complex logic in isolation and before the model is even working.

The model has the following decision variables: \[\color{darkred}x_{k,n} \in \{0,1\}\] indicates if item \((k,n)\) is placed in our container. Another binary variable \[\color{darkred}{\mathit{rot}}_{k,n} \in \{0,1\}\] tells us if item \((k,n)\) is rotated. In addition, we have a bunch of binary \(\color{darkred}\delta\) variables for the no-overlap constraints. Continuous variables are used to represent the location \[\color{darkred}{\mathit{pos}}_{k,n,xy} \ge 0\] and the end-point \[\color{darkred}{\mathit{pos2}}_{k,n,xy} \in [0,\color{darkblue}{\mathit{WH}}_{xy}]\] Here we use the set \(\color{darkblue}xy=\{x,y\}\) and the size of the container \(\color{darkblue}{\mathit{WH}}_{xy}\). Let's try to formulate a model:

Continuous size model
\[\begin{align} \max\> & \color{darkred}{\mathit{totalValue}} = \sum_{\color{darkblue}{\mathit{kn}}(k,n)} \color{darkblue}{\mathit{value}}_k \cdot \color{darkred}x_{k,n} && && \text{Maximize total value}\\ & \color{darkred}{\mathit{pos2}}_{k,n,xy}  = \color{darkred}{\mathit{pos}}_{k,n,xy} + \color{darkblue}{\mathit{wh}}_{k,R,xy}\cdot \color{darkred}{\mathit{rot}}_{k,n} + \color{darkblue}{\mathit{wh}}_{k,{\mathit{NR}},xy}\cdot (1-\color{darkred}{\mathit{rot}}_{k,n})  &&\forall \color{darkblue}{\mathit{kn}}_{k,n},xy&&\text{End-point of item} \\ & \color{darkred}{\mathit{pos}}_{k,n,xy} \ge \color{darkred}{\mathit{pos2}}_{k',n',xy} - \color{darkblue}M \cdot (3-\color{darkred}\delta_{k,n,k',n',xy,1}-\color{darkred}x_{k,n}-\color{darkred}x_{k',n'})&&\forall \color{darkblue}{\mathit{compare}}_{n,k,n',k'}, xy && \text{No overlap}\\ &\color{darkred}{\mathit{pos}}_{k',n',xy} \ge \color{darkred}{\mathit{pos2}}_{k,n,xy} - \color{darkblue}M \cdot (3-\color{darkred}\delta_{k,n,k',n',xy,2}-\color{darkred}x_{k,n}-\color{darkred}x_{k',n'})&&\forall \color{darkblue}{\mathit{compare}}_{n,k,n',k'},xy && \text{No overlap} \\ & \sum_{xy} \left( \color{darkred}\delta_{k,n,k',n',xy,1}+\color{darkred}\delta_{k,n,k',n',xy,2}\right)\ge 1 && \forall \color{darkblue}{\mathit{compare}}_{n,k,n',k'} && \text{No overlap}\\ &\color{darkred}x_{k,n} \in \{0,1\}\\ &\color{darkred}{\mathit{rot}}_{k,n} \in \{0,1\}\\ &\color{darkred}{\mathit{pos}}_{k,n,xy} \ge 0 \\ &\color{darkred}{\mathit{pos2}}_{k,n,xy} \in [0,\color{darkblue}{\mathit{WH}}_{xy}] \end{align}\]


We need good values for the big-M constants. I just used the width or height of the container. Note that it looks as if we have only 2 no-overlap big-M constraints for each pair of items. But if you look carefully, you can see that the set \(xy\) is used. So each constraint actually generates two conditions for each pair. 

I further added three constraints. They are not essential but may help performance.
  1. We can order the variable \(\color{darkred}x_{k,n}\) by \(n\): \[\color{darkred}x_{k,n-1}\ge\color{darkred}x_{k,n}\] This makes the solution easier to read.
  2. We can order variables \(\color{darkred}{\mathit{pos}}\) of the same type by their \(x\)-coordinate: \[\color{darkred}{\mathit{pos}}_{k,n-1,x}\le \color{darkred}{\mathit{pos}}_{k,n,x}\]
  3. We know that the total area of items cannot exceed the area of the container: \[\sum_{\color{darkblue}{\mathit{kn}}(k,n)}  \color{darkblue}{\mathit{area}}_k \cdot \color{darkred}x_{k,n} \le \color{darkblue}{\mathit{containerArea}}\]


The results of this model are:

----    127 VARIABLE x.L  item (k,n) is placed (rotated/non-rotated)

             n1          n2          n3          n4          n5          n6

k5            1
k6            1           1           1
k8            1           1
k9            1           1           1           1           1           1


----    127 PARAMETER positions  rectangles (x1,y1)-(x2,y2)

               x1          y1          x2          y2     rotated

k5.n1                               3.000       6.000
k6.n1      17.000                  30.000       5.000
k6.n2      20.000       5.000      25.000      18.000       1.000
k6.n3      25.000       5.000      30.000      18.000       1.000
k8.n1       2.000      12.000      20.000      18.000       1.000
k8.n2       2.000       6.000      20.000      12.000       1.000
k9.n1                   6.000       2.000      20.000       1.000
k9.n2       2.000      18.000      16.000      20.000
k9.n3       3.000                  17.000       2.000
k9.n4       3.000       4.000      17.000       6.000
k9.n5       3.000       2.000      17.000       4.000
k9.n6      16.000      18.000      30.000      20.000


----    127 VARIABLE totalValue.L          =     4617.938  objective (maximized)


We see the optimal objective is identical to the results of the covering model. Implementing two very different models is a good way to certify results. This model takes a little bit more than 200 seconds.

 
Solution by continuous size model


Notes


  • Slightly different data can cause very different performance.
  • The covering model is a bit faster than the continuous size model.
  • The problem is very simple. The MIP models not so much.
  • I did not try to physically duplicate data items to model rotation. Just have two items: one not rotated and one rotated. Impose a constraint that says we can only pick one. That may be a bit simpler than what I tried, especially in the second model. The Python code in the next section is exactly doing this.

Covering modelContinuous size model
Equations6116,561
Variables4,2735,407
Discrete Variables4,2725,202
Objective4617.93784617.9378
StatusOptimalOptimal
Time (sec.)40232
Iterations26,0669,598,476
Nodes63157,747
OptionsDefaultsDefaults


OR-Tools CP-SAT implementation 


In this section, I'll discuss my OR-Tools CP-SAT implementation of the optimization model. The main idea is to use the high-level NoOverlap2D constraint, but make the width and height of an item zero when it is not used.

The first thing I did was to physically create individual items (51 of them) and add explicitly rotated versions. So I have 102 items: the first 51 unrotated and the second 51 rotated. For the NoOverlap2D constraint, we need to create so-called interval variables for both the \(x\) and \(y\) directions. The interval construct has three components: a start, a length, and an end. In my approach, the start is the location of an item. The length is the width or height, but it can also be zero (so it is a variable). The end is the end-point of an item. 

A summary of the model can look like: 

CP-SAT Model
\[\begin{align}\max & \sum_i {\bf int}(1000 \cdot\color{darkblue}{\mathit{value}}_i)\cdot\color{darkred}u_i && && \text{integer coefficients}\\ & \color{darkred}w'_i = \color{darkblue}w_i \cdot \color{darkred}u_i && \forall i && \text{zero or }\color{darkblue}w_i\\ & \color{darkred}h'_i = \color{darkblue}h_i \cdot \color{darkred}u_i && \forall i && \text{zero or }\color{darkblue}h_i \\ & \color{darkred}u_i + \color{darkred}u_{i+51} \le 1 && i \in 1,\dots,51 && \text{rotated or not} \\ & {\bf NoOverlap2D}([\color{darkred}{x1},\color{darkred}w',\color{darkred}{x2}],[\color{darkred}{y1},\color{darkred}h',\color{darkred}{y2}]) && && \text{uses interval variables} \\ & \sum_i \color{darkblue}h_i \cdot \color{darkblue}w_i \cdot \color{darkred}u_i \le \color{darkblue}H \cdot \color{darkblue}W && && \text{extra} \\ & \color{darkred}{x1}_i,\color{darkred}{x2}_i \in \{0,\dots,\color{darkblue}H\} \\ & \color{darkred}{y1}_i,\color{darkred}{y2}_i \in \{0,\dots,\color{darkblue}W\} \\ &\color{darkred}w'_i \in \{0,\dots,\color{darkblue}w_i\} \\ &\color{darkred}h'_i \in \{0,\dots,\color{darkblue}h_i\}\\ & \color{darkred}u_i \in \{0,1\} \end{align}\]


The constants \(\color{darkblue}W\) and \(\color{darkblue}H\) are the width and height of the container. The interval variables make sure that \(\color{darkred}{x1}_i+\color{darkred}w'_i=\color{darkred}{x2}_i\) and  \(\color{darkred}{y1}_i+\color{darkred}h'_i=\color{darkred}{y2}_i\). The extra constraint was important for the performance. We see:


return code:4
status:OPTIMAL

CpSolverResponse summary:
status: OPTIMAL
objective: 4617934
best_bound: 4617934
booleans: 432
conflicts: 87347
branches: 125277
propagations: 362526
integer_propagations: 1855211
restarts: 269
lp_iterations: 209669
walltime: 142.573
usertime: 142.573
deterministic_time: 51.1133
primal_integral: 93.1079


     x   y   w   h  x2  y2       v
0    0  14  14   2  14  16  223516
1    0  16  14   2  14  18  223516
2    0  18  14   2  14  20  223516
3    2   0  14   2  16   2  223516
4   16   0  14   2  30   2  223516
5   20   2   5  13  25  15  551072
6   25   2   5  13  30  15  551072
7    2   2  18   6  20   8  755093
8    2   8  18   6  20  14  755093
9   14  14   3   6  17  20  113436
10  17  15  13   5  30  20  551072
11   0   0   2  14   2  14  223516


Note that we need to divide the objective by 1,000. We multiplied by a thousand to make the effect of the truncation to an integer value smaller. The time was 142 seconds, falling somewhere in between our MIP models. However, it is noted that (1) this was on a different machine (a small 2 core VM on colab.google.com) and (2) both performance (and the final solution) can vary a lot from one run to the next. Or-tools is less deterministic than MIP solvers like Cplex and Gurobi.
 

CP-SAT solution


This model is interesting as it is very different from the previous MIP models. There is very little reuse of the equations we developed earlier.  

Appendix 1: GAMS covering model


$ontext

  
2d knapsack problem

  
Fill a container with most valuable rectangles.

$offtext

set
  dummy
'for display' /xmin,xmax,ymin,ymax,width,height,available,value/

  i
'rows in container' /r1*r20/
  j
'columns in container' /c1*c30/
  k
'rectangular items to place in container' /k1*k10/
  r
'rotate' /nr,r/
;

*---------------------------------------------------------------------
* data
*---------------------------------------------------------------------

parameter data(*,*) 'problem data';
data(
'container','height') = card(i);
data(
'container','width') = card(j);
data(k,
'height') = uniformint(1,20);
data(k,
'width') = uniformint(1,20);
data(k,
'value/size') = uniform(1,10);
data(k,
'value') = data(k,'value/size')*data(k,'height')*data(k,'width');
data(k,
'available') = uniformint(1,10);
display data;

*---------------------------------------------------------------------
* derived data
*---------------------------------------------------------------------

alias (i,ii),(j,jj);

parameters
   w(k,r)
'item width'
   h(k,r)
'item height'
;

w(k,r) = data(k,
'width')$sameas(r,'nr') + data(k,'height')$sameas(r,'r');
h(k,r) = data(k,
'height')$sameas(r,'nr') + data(k,'width')$sameas(r,'r');

sets
  ok(k,r,i,j)
'item (k,r) can be placed at (i,j)'
  cover(k,r,i,j,ii,jj)
'(ii,jj) is covered when item (k,r) is placed at (i,j)'
;
ok(k,r,i,j) =
ord(i)+h(k,r)-1 <= card(i)
             
and ord(j)+w(k,r)-1 <= card(j);

cover(ok(k,r,i,j),ii,jj) =
          
ord(ii)>=ord(i) and ord(ii)<ord(i)+h(k,r) and
          
ord(jj)>=ord(j) and ord(jj)<ord(j)+w(k,r);


*-----------------------------------------------------------------------
* model
*-----------------------------------------------------------------------

binary variable x(k,r,i,j) 'item k is placed at (i,j)';
variable totalValue 'objective (maximized)';

equations
   noOverlap(i,j)
'only one item can cover (i,j)'
   count(k)      
'limit number of each item'
   objective     
'maximize total value'
;

noOverlap(ii,jj)..
sum(cover(k,r,i,j,ii,jj),x(k,r,i,j)) =l= 1;

count(k)..
sum(ok(k,r,i,j), x(k,r,i,j)) =l= data(k,'available');

objective.. totalValue =e=
sum(ok(k,r,i,j),x(k,r,i,j)*data(k,'value'));

model m /all/;
option optcr=0,threads=8;
solve m maximizing totalValue using mip;

*-----------------------------------------------------------------------
* reporting
*-----------------------------------------------------------------------

option x:0:2:2;
display x.l,totalValue.l;

parameter plotdata(k,r,i,j,*);
loop((k,r,i,j)$(x.l(k,r,i,j)>0.5),
    plotdata(k,r,i,j,
'xmin') = ord(j) - 1 + EPS;
    plotdata(k,r,i,j,
'xmax') = plotdata(k,r,i,j,'xmin') + w(k,r);
    plotdata(k,r,i,j,
'ymin') = ord(i) - 1 + EPS;
    plotdata(k,r,i,j,
'ymax') = plotdata(k,r,i,j,'ymin') + h(k,r);
    plotdata(k,r,i,j,
'value') = data(k,'value');
);
options plotdata:2:4:1;
display plotdata;




Appendix 2: GAMS continuous size model


$ontext

  
2d knapsack problem

  
Fill a container with most valuable rectangles.

  
Continuous size formulation: pairwise no-overlap constraints.

$offtext

*---------------------------------------------------------------
* data
*---------------------------------------------------------------

set
  dummy
'for display' /xmin,xmax,ymin,ymax,width,height,available,value/
  k
'rectangular items to place in container' /k1*k10/
  n
'item# for each k' /n1*n10/
  r
'rotate' /nr,r/
  xy
'coordinates' /x,y/
;

parameter data(*,*);
data(k,
'height') = uniformint(1,20);
data(k,
'width') = uniformint(1,20);
data(k,
'value/area') = uniform(1,10);
data(k,
'available') = uniformint(1,10);
data(
'container','height') = 20;
data(
'container','width') = 30;

data(k,
'area') = data(k,'height')*data(k,'width');
data(
'container','area') = data('container','height')*data('container','width');
data(k,
'value') = data(k,'value/area')*data(k,'area');
display data;

*---------------------------------------------------------------
* derived data
*---------------------------------------------------------------

parameter
   wh(k,r,xy)
'width or height'
   whn(k,n,r,xy)
'auxiliary'
   num(k,n)
'sequence numbers for (k,n)'
   cWH(xy)
'width or height of container'

;
wh(k,r,
'x') = data(k,'width')$sameas(r,'nr') + data(k,'height')$sameas(r,'r');
wh(k,r,
'y') = data(k,'height')$sameas(r,'nr') + data(k,'width')$sameas(r,'r');
whn(k,n,r,xy) = wh(k,r,xy);

cWH(xy) = data(
'container','width')$sameas(xy,'x') +
             data(
'container','height')$sameas(xy,'y');

set kn(k,n) '(k,n) combinations';
kn(k,n)$(
ord(n)<=data(k,'available')) = yes;
alias(kn,kn2);

scalar cnt /0/;
loop(kn,
   cnt = cnt+1;
   num(kn) = cnt;
);

option num:0;
display kn,num;

alias (k,kk),(n,nn),(r,rr);

set compare(k,n,kk,nn) 'we need no-overlap between these';
compare(kn,kn2) = num(kn) < num(kn2);


*---------------------------------------------------------------
* model
*---------------------------------------------------------------

binary variables
   x(k,n)
'item (k,n) is placed (rotated/non-rotated)'
   rot(k,n)
'item (k,n) is rotated'
   delta
'used in no-overlap constraints'
;
positive variables
    pos(k,n,xy)
'position'
    pos2(k,n,xy)
'pos+width/height'
;
variable totalValue 'objective (maximized)';

pos.up(kn(k,n),xy) = cWH(xy) -
smin(r,wh(k,r,xy));
pos2.up(kn(k,n),xy) = cWH(xy);

equations
   epos2(k,n,xy)          
'pos+width/height'
   noOverlap1(k,n,k,n,xy) 
'no-overlap between items'
   noOverlap2(k,n,k,n,xy) 
'no-overlap between items'
   noOverlap3(k,n,k,n)    
'no-overlap between items'
   objective

* these are extras. They may help performance
   order
   order2
   area
;

epos2(kn(k,n),xy)..
   pos2(kn,xy) =e= pos(kn,xy) + rot(kn)*wh(k,
'r',xy) + (1-rot(kn))*wh(k,'nr',xy);

noOverlap1(compare(kn,kn2),xy)..
   pos(kn,xy) =g= pos2(kn2,xy)
                  - cWH(xy)*(3-delta(kn,kn2,xy,
'1')-x(kn)-x(kn2));
noOverlap2(compare(kn,kn2),xy)..
   pos(kn2,xy) =g= pos2(kn,xy)
                   - cWH(xy)*(3-delta(kn,kn2,xy,
'2')-x(kn)-x(kn2));
noOverlap3(compare(kn,kn2))..
   
sum(xy, delta(kn,kn2,xy,'1')+delta(kn,kn2,xy,'2')) =g= 1;

objective.. totalValue =e=
sum(kn(k,n),x(kn)*data(k,'value'));

order(k,n)$(kn(k,n)
and kn(k,n-1)).. x(k,n-1) =g= x(k,n);
order2(k,n)$(kn(k,n)
and kn(k,n-1)).. pos(k,n-1,'x') =l= pos(k,n,'x');

area..
sum(kn(k,n), x(kn)*data(k,'area')) =l= prod(xy,cWH(xy));

model m /all/;
option optcr=0, threads=8;
solve m maximizing totalValue using mip;


*---------------------------------------------------------------
* reporting
*---------------------------------------------------------------

parameter positions(k,n,*) 'rectangles (x1,y1)-(x2,y2)';
positions(kn,
'x1')$(x.l(kn)>0.5) = pos.l(kn,'x');
positions(kn,
'y1')$(x.l(kn)>0.5) = pos.l(kn,'y');
positions(kn,
'x2')$(x.l(kn)>0.5) = pos2.l(kn,'x');
positions(kn,
'y2')$(x.l(kn)>0.5) = pos2.l(kn,'y');
positions(kn,
'rotated')$(x.l(kn)>0.5) = rot.l(kn);

option x:0;
display x.l,positions,totalValue.l;

parameter plotdata(k,n,*);
loop(kn(k,n)$(x.l(kn)>0.5),
   plotdata(kn,
'xmin') = pos.l(kn,'x')+EPS;
   plotdata(kn,
'xmax') = pos2.l(kn,'x');
   plotdata(kn,
'ymin') = pos.l(kn,'y')+EPS;
   plotdata(kn,
'ymax') = pos2.l(kn,'y');
   plotdata(kn,
'value') = data(k,'value');
);
option plotdata:2:2:1;
display plotdata;



Appendix 3: Python code for Or-tools CP-SAT model


from ortools.sat.python import cp_model
from io import StringIO
import pandas as pd
import numpy as np

#---------------------------------------------------
# data
#---------------------------------------------------

data =
'''
item         width    height available    value    color
k1             20       4       2        338.984   blue
k2             12      17       6        849.246   orange
k3             20      12       2        524.022   green
k4             16       7       9        263.303   red
k5              3       6       3        113.436   purple
k6             13       5       3        551.072   brown
k7              4       7       6         86.166   pink
k8              6      18       8        755.094   grey
k9             14       2       7        223.516   olive
k10             9      11       5        369.560   cyan
'''

df = pd.read_table(StringIO(data),sep=
'\s+')
print("Input data")
display(df)

H =
20
W =
30

print(f"Container width:{W} height:{H}")

#---------------------------------------------------
# derived data
# expand to individual items
#---------------------------------------------------

w0 = df[
'width'].to_numpy()
h0 = df[
'height'].to_numpy()
a0 = df[
'available'].to_numpy()
v0 = df[
'value'].to_numpy()
indx0 = np.arange(np.size(w0))

w = np.repeat(w0,a0)
h = np.repeat(h0,a0)
v = np.repeat(v0,a0)
indx = np.repeat(indx0,a0)

n = len(w)
print(f"Number of individual items: {n}")

#---------------------------------------------------
# create rotated items
# by duplicating
#---------------------------------------------------

wr = np.concatenate((w,h))
hr = np.concatenate((h,w))
vr = np.concatenate((v,v))
vr = (vr*
1000.0).astype('int'# scale values to become integers
indxr = np.concatenate((indx,indx))

nr =
2*n
print(f"Number of individual items (after adding rotations): {nr}")

#---------------------------------------------------
# or-tools model
#---------------------------------------------------

model = cp_model.CpModel()

#
# variables
#

# u[i] : item i is used
u = [ model.NewBoolVar(f
"u{i}") for i in range(nr) ]

# x[i],y[i] : location of item i
x = [ model.NewIntVar(
0,W,f"x{i}") for i in range(nr) ]
y = [ model.NewIntVar(
0,H,f"y{i}") for i in range(nr) ]

# xw[i],yh[i] : width/height item i
xw = [ model.NewIntVar(
0,int(wr[i]),f"xw{i}") for i in range(nr) ]
yh = [ model.NewIntVar(
0,int(hr[i]),f"yh{i}") for i in range(nr) ]


# x2[i],y2[i] : upper limit of interval variable
x2 = [ model.NewIntVar(
0,W,f"x2{i}") for i in range(nr) ]
y2 = [ model.NewIntVar(
0,H,f"y2{i}") for i in range(nr) ]

# interval variables
xival = [ model.NewIntervalVar(x[i],xw[i],x2[i],f
"xival{i}") for i in range(nr) ]
yival = [ model.NewIntervalVar(y[i],yh[i],y2[i],f
"yival{i}") for i in range(nr) ]

#
# constraints
#
# The main idea here is: if an item is not selected, make its
# height and width zero.

# u[i] = 0  ==> xw[i]=yh[i]=0
# u[i] = 1  ==> xw[i]=hr[i],yh[i]=wr[i]
for i in range(nr):
  model.Add(xw[i]==wr[i]*u[i])
  model.Add(yh[i]==hr[i]*u[i])

# only one of non-rotated/rotated pair can be used
for i in range(n):
    model.Add(u[i]+u[i+n]<=
1)

# no overlap.
model.AddNoOverlap2D(xival,yival)

# extra: this constraint helps performance enormously
model.Add(sum([wr[i]*hr[i]*u[i]
for i in range(nr)])<=W*H)

# objective
model.Maximize(sum([u[i]*vr[i]
for i in range(nr)]))


#
# solve model
#
solver = cp_model.CpSolver()
# log does not work inside a Jupyter notebook
solver.parameters.log_search_progress = True
solver.parameters.num_search_workers =
8
rc = solver.Solve(model)
print(f"return code:{rc}")
print(f"status:{solver.StatusName()}")
print()
print(solver.ResponseStats())
print()

#
# report solution
#
if rc == 4:
    used = {i
for i in range(nr) if solver.Value(u[i]) == 1 }
    df = pd.DataFrame({
       
'x'   : [solver.Value(x[i]) for i in used],
       
'y'   : [solver.Value(y[i]) for i in used],
       
'w'   : [wr[i] for i in used],
       
'h'   : [hr[i] for i in used],
       
'x2'  : [solver.Value(x2[i]) for i in used],
       
'y2'  : [solver.Value(y2[i]) for i in used],
       
'v'   : [vr[i] for i in used] })
   
print(df)

13 comments:

  1. Erwin, I looked at a similar model years ago. Despite the larger size, the covering MIP far outperformed the "continuous size" model. I believe this is because the continuous size model depends on disjunctive constraints; in the LP relaxation, disjunctive constraints tend to spread the values to 1/n, so the LP relaxation gives poor information. If necessary, I believe you could improve the performance of the covering MIP via symmetry reduction and/or column generation.

    ReplyDelete
    Replies
    1. That is certainly confirmed for these models. Sometimes the covering formulation can lead to really, really big models.

      Delete
  2. Hi,

    I tweaked a little the model (not sure it helps)

    ```for i in range(nr):
    model.Add(xw[i] == wr[i]).OnlyEnforceIf(u[i])
    model.Add(xw[i] == 0).OnlyEnforceIf(u[i].Not())
    model.Add(yh[i] == hr[i]).OnlyEnforceIf(u[i])
    model.Add(yh[i] == 0).OnlyEnforceIf(u[i].Not())
    ```


    On my machine, it solves quite fast

    CpSolverResponse summary:
    status: OPTIMAL
    objective: 4617936
    best_bound: 4617936
    booleans: 592
    conflicts: 30049
    branches: 51905
    propagations: 627691
    integer_propagations: 670265
    restarts: 237
    lp_iterations: 24505
    walltime: 6.20536
    usertime: 6.20536
    deterministic_time: 13.1242
    primal_integral: 106.057

    return code:4
    status:OPTIMAL

    CpSolverResponse summary:
    status: OPTIMAL
    objective: 4617936
    best_bound: 4617936
    booleans: 592
    conflicts: 30049
    branches: 51905
    propagations: 627691
    integer_propagations: 670265
    restarts: 237
    lp_iterations: 24505
    walltime: 6.20536
    usertime: 6.20536
    deterministic_time: 13.1242
    primal_integral: 106.057


    x y w h x2 y2 v
    0 0 0 14 2 14 2 223516
    1 0 2 14 2 14 4 223516
    2 0 4 14 2 14 6 223516
    3 14 18 14 2 28 20 223516
    4 18 5 5 13 23 18 551072
    5 23 5 5 13 28 18 551072
    6 0 18 14 2 14 20 223516
    7 0 6 18 6 18 12 755094
    8 14 0 3 6 17 6 113436
    9 17 0 13 5 30 5 551072
    10 0 12 18 6 18 18 755094
    11 28 5 2 14 30 19 223516

    ReplyDelete
    Replies
    1. Hi Laurent: that was my first formulation. I replaced it with a linear expression just because it was a little bit more compact. I ran it on colab.google.com. That may be a less "stable" environment (a small VM) than running it locally on a real machine.

      Delete
    2. In the morning it seems a bit slower:

      return code:4
      status:OPTIMAL

      CpSolverResponse summary:
      status: OPTIMAL
      objective: 4617934
      best_bound: 4617934
      booleans: 444
      conflicts: 1608380
      branches: 1653210
      propagations: 11072083
      integer_propagations: 38767100
      restarts: 1714
      lp_iterations: 0
      walltime: 331.133
      usertime: 331.133
      deterministic_time: 191.74
      primal_integral: 347.965

      Delete
  3. Using the old formulation does not change performance. Which version of or-tools did you use ?

    ReplyDelete
    Replies
    1. I did not print the version, but I see it installed: ortools-9.1.9490-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl

      Delete
  4. So, the jupyter backend seems to be mono-core, or close to.

    ReplyDelete
  5. I have rewritten the example, added a second model where the dimensions of the item depends on its state(not_selected, no_rotation, rotated).
    In the meantime, I have added a linear relaxation to the no_overlap_2d constraint that makes your equation obsolete.
    I have also implemented detection of possible linearization of the product of two affine expression, which is necessary to rebuild the linear term in the energetic redundant equation.

    This currently requires using the master branch of or-tools, while waiting for 9.2 to be released (hopefully this month).
    With this version, I solve the problem much faster than the previous code.

    The example is available at: https://github.com/google/or-tools/blob/master/examples/python/packing_2d_sat.py

    ReplyDelete
    Replies
    1. Last update. I have changed the name of the example to knapsack_2d_sat.py. And we have implemented support for floating point coefficients and offsets in the objective. See https://github.com/google/or-tools/blob/master/examples/python/knapsack_2d_sat.py

      Delete
  6. Problem: a timber warehouse has a bundle of assorted widths and lengths, thickness is uniform. A customer places an order for a number of boards of varying widths and lengths, thickness is uniform. Is there a mathematical solution or software which identifies the most efficient way of meeting the order from the available timber?

    ReplyDelete
    Replies
    1. This is a two-dimensional cutting stock problem. There is a body of literature about this problem.

      Delete
  7. can you provide the code to display the output in container form for the 3rd case(using cp-sat solver)?

    ReplyDelete