[Help-glpk] optimizing my ilp, stopping overlapping of tiles

[Paul Rubel]

Hi,

I have a tiling problem that I've formulated as an ILP. I even have
running mathprog for it. The problem is that my solution doesn't scale
very well. I can't scale past about 1000 tiles before it consumes too
much memory and crashes (sometimes after grinding for a week).

This was my first real exposure to LP and I suspect that the way I've
formulated the problem may not be as efficient as possible. I'm hoping
that there is some simple idea that an expert (or even a journeyman)
glpk user could help me with. 

The gist of the problem is that I have tiles on a grid. Each grid
square has an ideal value. I'd like to lay the tiles down to minimize
the difference between the tile values and the grid square values. So
far, no problem. The task is complicated by the fact that the tiles
come in pairs, like dominos, which is actually the domain. 

I need constraints that ensure that each square is covered by only one
tile. Currently I have a solution that tests each square more or less
exhaustively. For each grid square I have a constrain and check that
there is only a single tile in it:

var x {p in POSITIONS, d in TILES, o in ORIENTATION} >= 0, <= 1,integer;

/* In general a given non-edge square has 8 possibilities. It can go
out from that square or each of its neighbors can come in to it. For
example. 

      s.t. constraint1_2:sum {d in TILES} (0
         + x['P1_2',d,'UP']               % out
         + x['P1_2',d,'DOWN']             % out 
         + x['P1_2',d,'LEFT']             % out
         + x['P1_2',d,'RIGHT']            % out
         + x['P1_1',d,'UP']               % in
         + x['P1_3',d,'DOWN']             % in
         + x['P2_2',d,'LEFT']             % in
         + x['P0_2',d,'RIGHT']) = 1;      % in

For square 0,0, a special edge case, another tile can either be in
that square and either point up or right. Other tiles can come from
from the square above or left from the square to the right. */

s.t. constraint0_0:sum {d in TILES} (0
   + x['P0_0',d,'UP']
   + x['P0_0',d,'RIGHT']
   + x['P0_1',d,'DOWN']
   + x['P1_0',d,'LEFT']) = 1;


You can imagine that a solution like this gets large very quickly. I
keep thinking that there must be a better way to say this but I can't
figure it out. 

One idea was to split the squares into groups (of say 30 tiles) and
assign each square in the group a power of 2 value. Then I could sum
the groups up and check that it was what I expected. This would entail
a table indexed by group and position that gave out values. I'd them
sum for each group and make sure it added up to the expected value. I
couldn't find a way to say that a set was a union of a number of other
sets. Is this possible, does it sound reasonable? If I say that 
  set T := a b c z
  set V := w x y z
is the z in T and V the same if I use it as an index? This still seems
pretty hacky. Does anyone have a better idea? A lookup table of 
position-orientation x position-orientation is just too big as the number of 
squares increases. 

       Thanks for you help and the great software,
        Paul


A simple running example for a 2x2 grid: 

/* to run: glpsol -m small.mod -o small.sol */

/* sets */
set TILES;
set POSITIONS;
set ORIENTATION;


/* parameters */
param Cost {p in POSITIONS, d in TILES, o in ORIENTATION};



/* decision vars, is a tile in a position */
var x {p in POSITIONS, d in TILES, o in ORIENTATION} >= 0, <= 1,integer;


minimize layout: sum{p in POSITIONS, d in TILES, 
                     o in ORIENTATION} x[p,d,o] * Cost[p,d,o];

/* only use the tiles once */
s.t. onlyOnce00 : sum{p in POSITIONS, o in ORIENTATION} x[p,'D0_0', o] = 1;
s.t. onlyOnce01 : sum{p in POSITIONS, o in ORIENTATION} x[p,'D0_1', o] = 1;


/* fill each square */
s.t. constraint0_0:sum {d in TILES} (0
   + x['P0_0',d,'UP']
   + x['P0_0',d,'RIGHT']
   + x['P0_1',d,'DOWN']
   + x['P1_0',d,'LEFT']) = 1;

s.t. constraint0_1:sum {d in TILES} (0
   + x['P0_1',d,'DOWN']
   + x['P0_1',d,'RIGHT']
   + x['P0_0',d,'UP']
   + x['P1_1',d,'LEFT']) = 1;

s.t. constraint1_0:sum {d in TILES} (0
   + x['P1_0',d,'UP']
   + x['P1_0',d,'LEFT']
   + x['P1_1',d,'DOWN']
   + x['P0_0',d,'RIGHT']) = 1;

s.t. constraint1_1:sum {d in TILES} (0
   + x['P1_1',d,'DOWN']
   + x['P1_1',d,'LEFT']
   + x['P1_0',d,'UP']
   + x['P0_1',d,'RIGHT']) = 1;


data;

/* which way does the tile go? */ 
set ORIENTATION := "DOWN" "UP" "RIGHT" "LEFT";

set TILES       := D0_0 D0_1; 
set POSITIONS   := P0_0 P0_1 P1_0 P1_1; 

param Cost :=
   [*,*,'DOWN']:    D0_0   D0_1 :=
        P0_0          10000   10000
        P0_1            0      1
        P1_0          10000   10000
        P1_1            1      0
   [*,*,'UP']: D0_0   D0_1 :=
        P0_0            0      1
        P0_1          100000   100000
        P1_0            1      2
        P1_1          100000   100000
   [*,*,'RIGHT']:  D0_0   D0_1 :=
        P0_0            1      1
        P0_1            0      0
        P1_0          10000   10000
        P1_1          10000   10000
   [*,*,'LEFT']: D0_0    D0_1 :=
        P0_0          10000   10000
        P0_1          10000   10000
        P1_0            1      2
        P1_1            0      1;





end;