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 model	Continuous size model
Equations	611	6,561
Variables	4,273	5,407
Discrete Variables	4,272	5,202
Objective	4617.9378	4617.9378
Status	Optimal	Optimal
Time (sec.)	40	232
Iterations	26,066	9,598,476
Nodes	63	157,747
Options	Defaults	Defaults

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)