Finding all the Combinations - N Rectangles inside the Square

The following model finds solutions within a couple of seconds:

%  Chuffed: 1.6s
%  CPLEX:   3.9s
%  Gecode:  1.5s

int: noOfRectangles = 6;
int: squareLen = 10;
int: Empty = 0;

set of int: Coords = 1..squareLen;
set of int: Rectangles = 1..noOfRectangles;

%  decision variables:
%  The square matrix
%  Every tile is either empty or belongs to one of the rectangles
array[Coords, Coords] of var Empty .. noOfRectangles: s;

%  the edges of the rectangles
array[Rectangles] of var Coords: top;
array[Rectangles] of var Coords: bottom;
array[Rectangles] of var Coords: left;
array[Rectangles] of var Coords: right;

%  function
function var Coords: getCoord(Coords: row, Coords: col, Rectangles: r, Coords: coord, Coords: defCoord) =
  if s[row, col] == r then coord else defCoord endif;
    
%  ----------------------<  constraints  >-----------------------------
  
%  Determine rectangle limits as minima/maxima of the rows and columns for the rectangles.
%  Note: A non-existing rectangle would have top=squareLen, bottom=1, left=squareLen, right=1
%  This leads to a negative size and is thus ruled-out.
constraint forall(r in Rectangles) (
  top[r] == min([ getCoord(row, col, r, row, squareLen) | row in Coords, col in Coords])
);
constraint forall(r in Rectangles) (
  bottom[r] == max([ getCoord(row, col, r, row, 1) | row in Coords, col in Coords])
);
constraint forall(r in Rectangles) (
  left[r] == min([ getCoord(row, col, r, col, squareLen) | row in Coords, col in Coords])
);
constraint forall(r in Rectangles) (
  right[r] == max([ getCoord(row, col, r, col, 1) | row in Coords, col in Coords])
);

%  all tiles within the limits must belong to the rectangle
constraint forall(r in Rectangles) (
  forall(row in top[r]..bottom[r], col in left[r]..right[r]) 
    (s[row, col] == r)
);

%  enforce a minimum size per rectangle
constraint forall(r in Rectangles) (
  (bottom[r] - top[r] + 1) * (right[r] - left[r] + 1) in 2 .. 9 
);

%  symmetry breaking: 
%  order rectangles according to their top/left corners
constraint forall(r1 in Rectangles, r2 in Rectangles where r2 > r1) (
 (top[r1]*squareLen + left[r1]) < (top[r2]*squareLen + left[r2])
);


%  output solution

output [ if col == 1 then "\n" else "" endif ++ 
         if "\(s[row, col])" == "0" then "  " else "\(s[row, col]) " endif 
         | row in Coords, col in Coords];

The grid positions in the sqare can be empty or assume one of six values. The model determines the top and bottom rows of all rectangles. Together with the left and right columns, it makes sure that all tiles within these limits belong to the same rectangle.

To experiment, it is helpful to start with smaller square dimensions and/or smaller numbers of rectangles. It might also make sense to delimit the size of rectangles. Otherwise, the rectangles tend to become too small (1x1) or too big.

Symmetry breaking to enforce a certain ordering of rectangles, does speed-up the solving process.

Here's another solution using MiniZincs Geost constraint. This solution is heavily based on Patrick Trentins excellent answer here. Also make sure to see his explanation on the model.

I assume using the geost constraint speeds up the process a little. Symmetry breaking might further speed things up as Axel Kemper suggests.

include "geost.mzn";

int: k;
int: nObjects;
int: nRectangles;
int: nShapes; 

set of int: DIMENSIONS = 1..k;
set of int: OBJECTS    = 1..nObjects;
set of int: RECTANGLES = 1..nRectangles;
set of int: SHAPES     = 1..nShapes;

array[DIMENSIONS] of int:             l;
array[DIMENSIONS] of int:             u;
array[RECTANGLES,DIMENSIONS] of int:  rect_size;
array[RECTANGLES,DIMENSIONS] of int:  rect_offset;
array[SHAPES] of set of RECTANGLES:   shape;
array[OBJECTS,DIMENSIONS] of var int: x;
array[OBJECTS] of var SHAPES:         kind;


array[OBJECTS] of set of SHAPES: valid_shapes;

constraint forall (obj in OBJECTS) (
    kind[obj] in valid_shapes[obj]
);

constraint geost_bb(k, rect_size, rect_offset, shape, x, kind, l, u);

And the corresponding data:

k = 2;                 % Number of dimensions
nObjects = 6;          % Number of objects
nRectangles = 4;       % Number of rectangles
nShapes = 4;           % Number of shapes

l = [0, 0];            % Lower bound of our bounding box
u = [10, 10];          % Upper bound of our bounding box

rect_size = [|
     2, 3|
     3, 2|

     3, 5|
     5, 3|];
     

rect_offset = [|
     0, 0|
     0, 0|
     
     0, 0|
     0, 0|];
     
shape = [{1}, {2}, {3}, {4}];
    
    
valid_shapes = [{1, 2}, {1, 2}, {1, 2}, {1, 2}, {1, 2}, {3, 4}];

The output reads a little different. Take this example:

x = array2d(1..6, 1..2, [7, 0, 2, 5, 5, 0, 0, 5, 3, 0, 0, 0]);
kind = array1d(1..6, [1, 1, 1, 1, 1, 3]);

This means rectangle one is placed at [7, 0] and takes the shape [2,3] as seen in this picture:

Building on the answer of @Phonolog, one way to obtain the wanted output format is to use a 2d-array m which is mapped to x through constraints (here size is the bounding box size):

% mapping to a 2d-array output format
set of int: SIDE     = 0..size-1;

array[SIDE, SIDE] of var 0..nObjects: m;
    
constraint forall (i, j in SIDE) ( m[i,j] = sum(o in OBJECTS)(o * 
    (i >= x[o,1] /\ 
     i <= x[o,1] + rect_size[kind[o],1]-1 /\ 
     j >= x[o,2] /\ 
     j <= x[o,2] + rect_size[kind[o],2]-1)) );
    
% symmetry breaking between equal objects
array[OBJECTS] of var int: pos = [ size*x[o,1] + x[o,2] | o in OBJECTS ];

constraint increasing([pos[o] | o in 1..nObjects-1]);
    
solve satisfy;

output ["kind=\(kind)\n"] ++
["x=\(x)\n"] ++
["m="] ++ [show2d(m)]

Edit: Here is the complete code:

include "globals.mzn";

int: k = 2;
int: nObjects = 6;
int: nRectangles = 4;
int: nShapes = 4;
int: size = 10;

set of int: DIMENSIONS = 1..k;
set of int: OBJECTS    = 1..nObjects;
set of int: RECTANGLES = 1..nRectangles;
set of int: SHAPES     = 1..nShapes;

array[DIMENSIONS] of int:             l = [0, 0];
array[DIMENSIONS] of int:             u = [size, size];
array[OBJECTS,DIMENSIONS] of var int: x;
array[OBJECTS] of var SHAPES:         kind;

array[RECTANGLES,DIMENSIONS] of int: rect_size = [|
    3, 2|
    2, 3|
    5, 3|
    3, 5|];
array[RECTANGLES,DIMENSIONS] of int: rect_offset = [|
    0, 0|
    0, 0|
    0, 0|
    0, 0|];
array[SHAPES] of set of SHAPES: shape = [
    {1}, {2}, {3}, {4}];

array[OBJECTS] of set of SHAPES: valid_shapes = 
    [{1, 2}, {1, 2}, {1, 2}, {1, 2}, {1, 2}, {3, 4}];
    
constraint forall (obj in OBJECTS) (
    kind[obj] in valid_shapes[obj]
);

constraint
    geost_bb(
        k,
        rect_size,
        rect_offset,        
        shape,
        x,
        kind,
        l,
        u
    );
    
% mapping to a 2d-array output format
set of int: SIDE     = 0..size-1;

array[SIDE, SIDE] of var 0..nObjects: m;

constraint forall (i, j in SIDE) ( m[i,j] = sum(o in OBJECTS)(o * 
    (i >= x[o,1] /\ 
     i <= x[o,1] + rect_size[kind[o],1]-1 /\ 
     j >= x[o,2] /\ 
     j <= x[o,2] + rect_size[kind[o],2]-1)) );


% symmetry breaking between equal objects
array[OBJECTS] of var int: pos = [ size*x[o,1] + x[o,2] | o in OBJECTS ];

constraint increasing([pos[o] | o in 1..nObjects-1]);

solve satisfy;

output ["kind=\(kind)\n"] ++
["x=\(x)\n"] ++
["m="] ++ [show2d(m)]