Maximize enclosed area of given figures on 2d grid

Question

I need to solve an optimization problem for a given set of polyominoes, for example the five Tetrominoes known from Tetris. The goal is to place each one of the figures on the 2d grid, so the area they enclose is maximal.

Staying with our Tetromino example, the optimal solution would be:

Bruteforcing the problem becomes unfeaseable very quickly, especially since my check with bfs search alone that an area is even enclosed at all must be performed countless times. I tried to implement systems to recognize when it won't be possible to close the circle with the remaining polyominoes, but didn't achieve a significant speed-up with them.

I also researched into strategies for the related polyomino tiling problem to see if I can adapt them, but had no success so far. The most promising approach I found was to translate it into Mixed-Integer Linear Programming like here or in this variation. But to apply this approach to my problem I need to both add the restriction of the enclosed area and then formulate its size as a linear objective function, whereby especially the encoding of the latter may not be possible at all.

Are there other algorithms to find the global maximum for this type of optimization problem that I've overlooked? Or are there data structures or strategies for branch elimination to drastically speed up the search of all possibilities?

Answer 1

Here is one possible approach. Introduce boolean variables $x_\ell,y_\ell,z_{\ell,p,r}$, where $\ell$ ranges over locations in the grid, $p$ ranges over polyominoes, and $r$ over rotation amounts (0, 90, 180, 270), with the following meaning:

• $x_\ell$ is true means that $\ell$ is part of the interior region that is enclosed (not including the boundary/perimeter)

• $y_\ell$ is true means that $\ell$ is part of the boundary/perimeter (it is covered by a polyomino)

• $z_{\ell,p,r}$ means that polyomino $p$ is rotated by $r$ degrees and then placed at location $\ell$

Now you can write down some boolean constraints regarding the relationships between these. For instance:

• $x_\ell \land \neg y_{\ell'} \implies x_{\ell'}$, for all pairs of adjacent locations $\ell,\ell'$

• $z_{\ell,p,r} \implies y_{\ell+t}$ for all offsets $t$ that are part of $p$ after being rotated by $r$ degrees

• $y_{\ell} \implies \bigvee_{p,r,t} z_{\ell-t,p,r}$ where the logical-OR ranges over all $p,r$ and all offsets $t$ that are part of $p$ after being rotated by $r$ degrees

• no two polyominoes overlap, which can be expressed in terms of $\neg (z_{\ell,p,r} \land z_{\ell',p',r'})$ where $p,p',r,r'$ range over all possibilities and then $\ell,\ell'$ have an appropriate relationship given $p,p',r,r'$

• For each $p$, exactly one of $z_{\ell,p,r}$ is true (ranging over all $\ell,r$) - this is a one-out-of-$n$ constraint, which can be encoded in standard ways

• $\neg x_\ell$, for each location $\ell$ along the four outermost edges of the grid of cells - this forces them to be on on the outside of the shape

• $\neg y_\ell$, for each location $\ell$ along the four outermost edges of the grid of cells - this forces them to not be covered by any polyomino

Finally the goal is to find a satisfying assignment that maximizes the number of $x_\ell$ that are true. This could be solved using ILP (each variable can be treated as a 0-or-1 integer variable, each boolean constraint can be expressed as one or more linear inequalities, and the objective function is $\sum_\ell x_\ell$) or using SAT (use binary search on the number of cells the region encloses, add a pseudo-boolean constraint that requires at least that many of the $x_\ell$ are true). I don't know whether this will be efficient enough to work well in practice, but you could give it a try. That would let you benefit from all of the heuristics and algorithms that have been developed over the years for ILP and/or SAT.

Answer 2

In case anyone is interested, this is the solution implemented with https://github.com/afish/MilpManager

It's inspired by https://cs.stackexchange.com/a/153823/79175 (including the discussions we had in the comments). Source code:

var integerWidth = 15;
var epsilon = 0.001;
var storeDebugExpressions = false;
var cacheConstants = false;
var storeDebugConstraints = false;
var solver = new CplexMilpSolver(new CplexMilpSolverSettings
{
    IntegerWidth = integerWidth,
    Epsilon = epsilon,
    StoreDebugExpressions = storeDebugExpressions,
    CacheConstants = cacheConstants,
    StoreDebugConstraints = storeDebugConstraints
});
solver.Cplex.SetParam(ILOG.CPLEX.Cplex.Param.MIP.Tolerances.Integrality, 0.000000000001);
        
var pieces = new []{
    new []{
        "XX",
        "X ",
        "X "
    },
    new[]{
        " XX",
        "XX "
    },
    new[]{
        "XX ",
        " XX"
    },
    new[]{
        "X ",
        "XX",
        "X "
    },
    new[]{
        "XX",
        "XX"
    },
    new[]{
        "XXXX"
    },
    new[]{
        "X ",
        "X ",
        "XX"
    }
};

var totalLength = pieces.Select(p => (int)Math.Max(p.Length, p[0].Length)).Sum();
var maxLength = pieces.Max(p => (int)Math.Max(p.Length, p[0].Length));
var width = 11; // Or totalLength + 2 to make sure it's big enough for any shapes
var height = 9; // Or totalLength + 2 to make sure it's big enough for any shapes
var piecesCount = pieces.Length;
// Right, down, left, up
var directions = 4; 

var heads = Enumerable.Range(0, height).Select(h =>
    Enumerable.Range(0, width).Select(w => 
        Enumerable.Range(0, piecesCount).Select(p => 
            Enumerable.Range(0, directions).Select(d =>
                solver.CreateAnonymous(Domain.BinaryInteger)
            ).ToArray()
        ).ToArray()
    ).ToArray()
).ToArray();

var isInterior = Enumerable.Range(0, height).Select(h =>
    Enumerable.Range(0, width).Select(w => 
        solver.CreateAnonymous(Domain.BinaryInteger)
    ).ToArray()
).ToArray();

var isBoundary = Enumerable.Range(0, height).Select(h =>
    Enumerable.Range(0, width).Select(w => 
        solver.CreateAnonymous(Domain.BinaryInteger)
    ).ToArray()
).ToArray();

Func<string[], int, string[]> rotatePiece = (piece, direction) => {
    if(direction == 0){
        return piece;
    }else if(direction == 2){
        return piece.Select(l => string.Join("", l.Reverse())).Reverse().ToArray();
    }else if(direction == 1){
        return Enumerable.Range(0, piece[0].Length)
            .Select(c => string.Join("", Enumerable.Range(0, piece.Length).Select(w => piece[piece.Length - 1 - w][c])))
            .ToArray();
    }else {
        return Enumerable.Range(0, piece[0].Length)
            .Select(c => string.Join("", Enumerable.Range(0, piece.Length).Select(w => piece[w][piece[0].Length - 1 - c])))
            .ToArray();
    }
};

// Put pieces on the board
for(int p=0;p<piecesCount;++p){
    solver.Set<Equal>(
        solver.FromConstant(1),
        solver.Operation<Addition>(
            Enumerable.Range(0, height).SelectMany(h =>
                Enumerable.Range(0, width).SelectMany(w =>
                    Enumerable.Range(0, directions).Select(d =>
                        heads[h][w][p][d]
                    )
                )
            ).ToArray()
        )
    );
}

// Make sure pieces don't fall off the board
for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        for(int p=0;p<piecesCount;++p){
            for(int d = 0; d < directions;++d){
                var piece = rotatePiece(pieces[p], d);
                if(d == 0){
                    if(w + piece[0].Length - 1 >= width){
                        heads[h][w][p][d].Set<Equal>(solver.FromConstant(0));
                    }
                }else if(d == 1){
                    if(h + piece.Length - 1 >= height){
                        heads[h][w][p][d].Set<Equal>(solver.FromConstant(0));
                    }
                }else if(d == 2){
                    if(w - piece[0].Length + 1 < 0){
                        heads[h][w][p][d].Set<Equal>(solver.FromConstant(0));
                    }
                }else if(d == 3){
                    if(h - piece.Length + 1 < 0){
                        heads[h][w][p][d].Set<Equal>(solver.FromConstant(0));
                    }
                }
            }
        }
    }
}

// Set impacts
IList<Tuple<int, int, int, int>>[][] allImpacts = Enumerable.Range(0, height).Select(h =>
    Enumerable.Range(0, width).Select(w => 
        new List<Tuple<int, int, int, int>>()
    ).ToArray()
).ToArray();

for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        for(int p=0;p<piecesCount;++p){
            for(int d = 0; d < directions;++d){
                var piece = rotatePiece(pieces[p], d);
                var pieceWidth = piece[0].Length;
                var pieceHeight = piece.Length;
                for(int y=0;y<pieceHeight;++y){
                    for(int x=0;x<pieceWidth;++x){
                        var isPieceX = piece[y][x] == 'X';
                        if(!isPieceX) continue;
                        
                        var finalY = d == 0 || d == 1 ? h + y : y - pieceHeight + 1 + y;
                        var finalX = d == 0 || d == 1 ? w + x : x - pieceWidth + 1 + x;
                        
                        if(finalY >= 0 && finalY < height && finalX >=0 && finalX < width){
                            allImpacts[finalY][finalX].Add(Tuple.Create(h, w, p, d));
                        }
                    }
                }
            }
        }
    }
}

// Set the boundary
for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        isBoundary[h][w].Set<Equal>(
            solver.Operation<Addition>(
                allImpacts[h][w].Select(t => heads[t.Item1][t.Item2][t.Item3][t.Item4]).ToArray()
            )
        );
    }
}

// Make sure pieces do not overlap
for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        solver.Set<LessOrEqual>(
            solver.Operation<Addition>(
                allImpacts[h][w].Select(t => heads[t.Item1][t.Item2][t.Item3][t.Item4]).ToArray()
            ),
            solver.FromConstant(1)
        );
    }
}


// Extend interior
for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        for(int y=-1;y<=1;++y){
            for(int x=-1;x<=1;++x){
                if(y == 0 && x == 0){
                    continue;
                }
                
                if(h + y < 0 || h + y >= height){
                    continue;
                }
                
                if(w + x < 0 || w + x >= width){
                    continue;
                }
                
                solver.Set<Equal>(
                    solver.FromConstant(1),
                    solver.Operation<MaterialImplication>(
                        solver.Operation<Conjunction>(
                            isInterior[h][w],
                            isBoundary[h + y][w + x].Operation<BinaryNegation>()
                        ),
                        isInterior[h + y][w + x]
                    )
                );
            }
        }
    }
}

// Make sure there is some exterior
for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        if(h == 0 || h == height - 1 || w == 0 || w == width - 1){
            solver.Set<Equal>(
                solver.FromConstant(0),
                isInterior[h][w]
            );
            solver.Set<Equal>(
                solver.FromConstant(0),
                isBoundary[h][w]
            );
        }
    }
}

// Make sure a field cannot be both exterior and boundary
for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        solver.Set<LessOrEqual>(
            solver.Operation<Addition>(
                isInterior[h][w],
                isBoundary[h][w]
            ),
            solver.FromConstant(1)
        );
    }
}

var goal = solver.Operation<Addition>(isInterior.SelectMany(x => x).ToArray());
solver.AddGoal("goal", goal);
solver.Solve();

Console.WriteLine(solver.GetValue(goal));

var board = Enumerable.Range(0, height).Select(h => new String(' ', width).ToCharArray()).ToArray();

for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        for(int p=0;p<piecesCount;++p){
            for(int d = 0; d < directions;++d){
                var piece = rotatePiece(pieces[p], d);
                var pieceWidth = piece[0].Length;
                var pieceHeight = piece.Length;
                if(solver.GetValue(heads[h][w][p][d]) > 0.5){
                    for(int y=0;y<pieceHeight;++y){
                        for(int x=0;x<pieceWidth;++x){
                            var finalY = d == 0 || d == 1 ? h + y : y - pieceHeight + 1 + y;
                            var finalX = d == 0 || d == 1 ? w + x : x - pieceWidth + 1 + x;
                            
                            if(finalY >= 0 && finalY < height && finalX >=0 && finalX < width){
                                board[finalY][finalX] = piece[y][x];
                            }else {
                                Console.WriteLine("Error in " + h + " " + w + " " + p  + " " + d  + " " + y + " " + x);
                            }
                        }
                    }
                }
            }
        }
    }
}

foreach(var l in board){
    Console.WriteLine(string.Join("", l));
}
Console.WriteLine();


for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        Console.Write(Math.Round(solver.GetValue(isInterior[h][w])));
    }
    Console.WriteLine();
}
Console.WriteLine();


for(int h = 0; h < height; ++h){
    for(int w = 0; w < width; ++ w){
        Console.Write(Math.Round(solver.GetValue(isBoundary[h][w])));
    }
    Console.WriteLine();
}
Console.WriteLine();

It should generalize to any polyominoes you have (just replace the array at the beginning). This can be easily translated to SAT if you need, so you have plenty of solvers out there to play with it.

CPLEX output for the sample above:

Tried aggregator 2 times.
MIP Presolve eliminated 4659 rows and 4747 columns.
MIP Presolve modified 688 coefficients.
Aggregator did 1432 substitutions.
Reduced MIP has 809 rows, 847 columns, and 4647 nonzeros.
Reduced MIP has 847 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.03 sec. (10.92 ticks)
Found incumbent of value 0.000000 after 0.06 sec. (33.44 ticks)
Probing time = 0.00 sec. (2.70 ticks)
Tried aggregator 1 time.
Reduced MIP has 809 rows, 847 columns, and 4647 nonzeros.
Reduced MIP has 847 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.02 sec. (2.75 ticks)
Probing time = 0.02 sec. (2.69 ticks)
Clique table members: 1918.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 8 threads.
Root relaxation solution time = 0.02 sec. (17.35 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                            0.0000       35.0000      565     ---
      0     0       29.7678   243        0.0000       29.7678      565     ---
*     0+    0                            2.0000       29.7678      657     ---
      0     0       29.2337   242        2.0000      Cuts: 36      657     ---
*     0+    0                            3.0000       29.2337      657  874.46%
      0     0       28.7133   242        3.0000      Cuts: 99      816  857.11%
      0     0       28.5944   248        3.0000      Cuts: 26      881  853.15%
      0     0       28.3575   243        3.0000      Cuts: 21      958  845.25%
      0     0       28.1942   256        3.0000      Cuts: 27     1036  839.81%
      0     0       28.1076   244        3.0000      Cuts: 25     1120  836.92%
*     0+    0                            4.0000       28.1076     1120  602.69%
      0     0       28.1076   245        4.0000      Cuts: 12     1130  602.69%
      0     0       28.0851   255        4.0000       Cuts: 9     1150  597.76%
      0     0       27.9794   239        4.0000      Cuts: 28     1229  597.76%
      0     0       27.8878   280        4.0000      Cuts: 25     1330  597.20%
      0     0       27.8093   275        4.0000      Cuts: 23     1402  595.23%
      0     0       27.6921   275        4.0000      Cuts: 40     1520  592.30%
      0     0       27.6921   276        4.0000       Cuts: 7     1525  592.30%
*     0+    0                           20.0000       27.6921     1525   38.46%
      0     2       27.6921   276       20.0000       27.6060     1525   38.03%
Elapsed time = 0.97 sec. (335.52 ticks, tree = 0.01 MB, solutions = 5)
*     5+    5                           21.0000       27.5336     1944   31.11%
*    56    40      integral     0       22.0000       27.4817     5646   24.92%
*    63+   38                           25.0000       27.4817     5959    9.93%

Clique cuts applied:  15
Cover cuts applied:  19
Implied bound cuts applied:  63
Mixed integer rounding cuts applied:  35
Zero-half cuts applied:  53
Gomory fractional cuts applied:  1

Root node processing (before b&c):
  Real time             =    0.97 sec. (335.26 ticks)
Parallel b&c, 8 threads:
  Real time             =    0.52 sec. (250.13 ticks)
  Sync time (average)   =    0.26 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    1.49 sec. (585.39 ticks)
25

   XXXXXX
  XX    X
 XX     X
 X      XX
 X      X
 XX    XX
  XXXXXXX


00000000000
00000000000
00001111000
00011111000
00111111000
00111111000
00011110000
00000000000
00000000000

00000000000
00011111100
00110000100
01100000100
01000000110
01000000100
01100001100
00111111100
00000000000