ILP Part 80 — Islanders – Random IT Utensils

This is the eightieth part of the ILP series. For your convenience you can find other parts in the table of contents in Part 1 – Boolean algebra

Today we are going to solve Islanders game.

We have different types of buildings. Each building has a size and a range. When building is constructed, it automatically adds points for its base gain and for all buildings in its range constructed earlier (which may be negative). We want to maximize the income.

Let’s start with model:

class Building{
	public string Name {get; set;}
	public int Range {get;set;}
	public int Size {get; set;}
	public int BaseScore {get;set;}
}

class Building{

public string Name {get; set;}

public int Range {get;set;}

public int Size {get; set;}

public int BaseScore {get;set;}

}

Each building has a name (for debugging only), range and size, and a base score (points added when the building is constructed).

Next, for each building we’ll create an ILP model:

class Placement{
	public Building Building {get;set;}
	// Center of building
	public IVariable X {get;set;}
	public IVariable Y {get;set;}
	public IVariable BuildOrder {get;set;}
	public int Identifier {get;set;}
}

class Placement{

public Building Building {get;set;}

// Center of building

public IVariable X {get;set;}

public IVariable Y {get;set;}

public IVariable BuildOrder {get;set;}

public int Identifier {get;set;}

}

This has two variables for position, and third variable for build order. It also has an identifier for debugging.

Now, the code:

// Model
var cityCenter = new Building(){
	Name = "City Center",
	Range = 5,
	Size = 1,
	BaseScore = 15
};
var house = new Building(){
	Name = "House",
	Range = 6,
	Size = 2,
	BaseScore = 1
};
var mansion = new Building(){
	Name = "Mansion",
	Range = 3,
	Size = 2,
	BaseScore = 1
};
var fountain = new Building(){
	Name = "Fountain",
	Range = 8,
	Size = 2,
	BaseScore = 0
};

var costs = new Dictionary<Tuple<Building, Building>, int>(){
	{ Tuple.Create(cityCenter, cityCenter), -40},
	{ Tuple.Create(cityCenter, house), 0},
	{ Tuple.Create(cityCenter, mansion), 0},
	{ Tuple.Create(cityCenter, fountain), 0},
	{ Tuple.Create(house, cityCenter), 6},
	{ Tuple.Create(house, house), 1},
	{ Tuple.Create(house, mansion), 0},
	{ Tuple.Create(house, fountain), 2},
	{ Tuple.Create(mansion, cityCenter), 8},
	{ Tuple.Create(mansion, house), 0},
	{ Tuple.Create(mansion, mansion), 1},
	{ Tuple.Create(mansion, fountain), 3},
	{ Tuple.Create(fountain, cityCenter), 7},
	{ Tuple.Create(fountain, house), 2},
	{ Tuple.Create(fountain, mansion), 3},
	{ Tuple.Create(fountain, fountain), -15}
};

Building[] buildings = new Dictionary<Building, int>(){
	{ cityCenter, 1},
	{ house, 4}, 
	{ mansion, 0},
	{ fountain, 0}
}.SelectMany(kv => Enumerable.Range(0, kv.Value).Select(i => kv.Key)).ToArray();

var boardSize = 20;

// ILP
Placement[] placements = buildings.Select((b, i) => new Placement(){
	Building = b,
	X = solver.CreateAnonymous(Domain.PositiveOrZeroInteger),
	Y = solver.CreateAnonymous(Domain.PositiveOrZeroInteger),
	BuildOrder = solver.CreateAnonymous(Domain.PositiveOrZeroInteger),
	Identifier = i + 1
}).ToArray();

// Make sure buildings are on board
foreach(var p in placements){
	p.X.Set<LessOrEqual>(solver.FromConstant(boardSize - p.Building.Size));
	p.Y.Set<LessOrEqual>(solver.FromConstant(boardSize - p.Building.Size));
}

// Make sure buildings do not overlap		
foreach(var p in placements){
	foreach(var p2 in placements){
		if(p == p2)continue;
		
		var overlapX = p.X.Operation<Subtraction>(p2.X).Operation<AbsoluteValue>().Operation<Multiplication>(solver.FromConstant(2)).Operation<IsLessThan>(solver.FromConstant(p.Building.Size + p2.Building.Size));
		
		var overlapY = p.Y.Operation<Subtraction>(p2.Y).Operation<AbsoluteValue>().Operation<Multiplication>(solver.FromConstant(2)).Operation<IsLessThan>(solver.FromConstant(p.Building.Size + p2.Building.Size));
		
		var overlap = overlapX.Operation<Disjunction>(overlapY);
		
		overlap.Set<Equal>(solver.FromConstant(0));
	}
}

// Make sure building order is different
solver.Set<AllDifferent>(solver.FromConstant(0), placements.Select(p => p.BuildOrder).ToArray());

// Calculate cost
IVariable baseInput = solver.Operation<Addition>(placements.Select(p => p.Building.BaseScore).Select(s => solver.FromConstant(s)).ToArray());
IVariable[] rangeInput = placements.SelectMany(p => {
	return placements.Where(p2 => p != p2).Select(p2 => {
		var isAfter = p.BuildOrder.Operation<IsGreaterOrEqual>(p2.BuildOrder);
		var isInRange = p.X.Operation<Subtraction>(p2.X).Operation<AbsoluteValue>().Operation<Addition>(p.Y.Operation<Subtraction>(p2.Y).Operation<AbsoluteValue>()).Operation<IsLessOrEqual>(solver.FromConstant(p.Building.Range));
		
		return solver.Operation<Condition>(
			isAfter.Operation<Conjunction>(isInRange),
			solver.FromConstant(costs[Tuple.Create(p.Building, p2.Building)]),
			solver.FromConstant(0)
		);
	});
}).ToArray();

var goal = baseInput.Operation<Addition>(rangeInput);

solver.AddGoal("Goal", goal);

solver.Solve();

Console.WriteLine("Result: " + solver.GetValue(goal));

var board = Enumerable.Range(0, boardSize).Select(i => new int[boardSize]).ToArray();
for(int i=0;i<boardSize;++i){
	for(int j=0;j<boardSize;++j){
		board[i][j] = 0;
	}
}
foreach(var p in placements){
	int x = (int)Math.Round(solver.GetValue(p.X));
	int y = (int)Math.Round(solver.GetValue(p.Y));
	for(int i=0;i<p.Building.Size;++i){
		for(int j=0;j<p.Building.Size;++j){
			board[x+i][y+j] = p.Identifier;
		}
	}
}
for(int i=0;i<boardSize;++i){
	for(int j=0;j<boardSize;++j){
		if(board[i][j] != 0){
			Console.Write(board[i][j]);
		}else{
			Console.Write(" ");
		}
	}
	Console.WriteLine();
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

// Model

var cityCenter = new Building(){

Name = "City Center",

Range = 5,

Size = 1,

BaseScore = 15

};

var house = new Building(){

Name = "House",

Range = 6,

Size = 2,

BaseScore = 1

};

var mansion = new Building(){

Name = "Mansion",

Range = 3,

Size = 2,

BaseScore = 1

};

var fountain = new Building(){

Name = "Fountain",

Range = 8,

Size = 2,

BaseScore = 0

};

var costs = new Dictionary<Tuple<Building, Building>, int>(){

{ Tuple.Create(cityCenter, cityCenter), -40},

{ Tuple.Create(cityCenter, house), 0},

{ Tuple.Create(cityCenter, mansion), 0},

{ Tuple.Create(cityCenter, fountain), 0},

{ Tuple.Create(house, cityCenter), 6},

{ Tuple.Create(house, house), 1},

{ Tuple.Create(house, mansion), 0},

{ Tuple.Create(house, fountain), 2},

{ Tuple.Create(mansion, cityCenter), 8},

{ Tuple.Create(mansion, house), 0},

{ Tuple.Create(mansion, mansion), 1},

{ Tuple.Create(mansion, fountain), 3},

{ Tuple.Create(fountain, cityCenter), 7},

{ Tuple.Create(fountain, house), 2},

{ Tuple.Create(fountain, mansion), 3},

{ Tuple.Create(fountain, fountain), -15}

};

Building[] buildings = new Dictionary<Building, int>(){

{ cityCenter, 1},

{ house, 4},

{ mansion, 0},

{ fountain, 0}

}.SelectMany(kv => Enumerable.Range(0, kv.Value).Select(i => kv.Key)).ToArray();

var boardSize = 20;

// ILP

Placement[] placements = buildings.Select((b, i) => new Placement(){

Building = b,

X = solver.CreateAnonymous(Domain.PositiveOrZeroInteger),

Y = solver.CreateAnonymous(Domain.PositiveOrZeroInteger),

BuildOrder = solver.CreateAnonymous(Domain.PositiveOrZeroInteger),

Identifier = i + 1

}).ToArray();

// Make sure buildings are on board

foreach(var p in placements){

p.X.Set<LessOrEqual>(solver.FromConstant(boardSize - p.Building.Size));

p.Y.Set<LessOrEqual>(solver.FromConstant(boardSize - p.Building.Size));

}

// Make sure buildings do not overlap

foreach(var p in placements){

foreach(var p2 in placements){

if(p == p2)continue;

var overlapX = p.X.Operation<Subtraction>(p2.X).Operation<AbsoluteValue>().Operation<Multiplication>(solver.FromConstant(2)).Operation<IsLessThan>(solver.FromConstant(p.Building.Size + p2.Building.Size));

var overlapY = p.Y.Operation<Subtraction>(p2.Y).Operation<AbsoluteValue>().Operation<Multiplication>(solver.FromConstant(2)).Operation<IsLessThan>(solver.FromConstant(p.Building.Size + p2.Building.Size));

var overlap = overlapX.Operation<Disjunction>(overlapY);

overlap.Set<Equal>(solver.FromConstant(0));

}

// Make sure building order is different

solver.Set<AllDifferent>(solver.FromConstant(0), placements.Select(p => p.BuildOrder).ToArray());

// Calculate cost

IVariable baseInput = solver.Operation<Addition>(placements.Select(p => p.Building.BaseScore).Select(s => solver.FromConstant(s)).ToArray());

IVariable[] rangeInput = placements.SelectMany(p => {

return placements.Where(p2 => p != p2).Select(p2 => {

var isAfter = p.BuildOrder.Operation<IsGreaterOrEqual>(p2.BuildOrder);

var isInRange = p.X.Operation<Subtraction>(p2.X).Operation<AbsoluteValue>().Operation<Addition>(p.Y.Operation<Subtraction>(p2.Y).Operation<AbsoluteValue>()).Operation<IsLessOrEqual>(solver.FromConstant(p.Building.Range));

return solver.Operation<Condition>(

isAfter.Operation<Conjunction>(isInRange),

solver.FromConstant(costs[Tuple.Create(p.Building, p2.Building)]),

solver.FromConstant(0)

);

});

}).ToArray();

var goal = baseInput.Operation<Addition>(rangeInput);

solver.AddGoal("Goal", goal);

solver.Solve();

Console.WriteLine("Result: " + solver.GetValue(goal));

var board = Enumerable.Range(0, boardSize).Select(i => new int[boardSize]).ToArray();

for(int i=0;i<boardSize;++i){

for(int j=0;j<boardSize;++j){

board[i][j] = 0;

}

foreach(var p in placements){

int x = (int)Math.Round(solver.GetValue(p.X));

int y = (int)Math.Round(solver.GetValue(p.Y));

for(int i=0;i<p.Building.Size;++i){

for(int j=0;j<p.Building.Size;++j){

board[x+i][y+j] = p.Identifier;

}

for(int i=0;i<boardSize;++i){

for(int j=0;j<boardSize;++j){

if(board[i][j] != 0){

Console.Write(board[i][j]);

}else{

Console.Write(" ");

}

Console.WriteLine();

}

In lines 1-25 we define buildings. Next, in lines 27-44 we specify points for constructing buildings in ranges. Finally, we ask to build one city center and for houses, just to keep model simple. We also specify that our board side length is 20.

We initialize variables in lines 55-62. Then, we start adding constraints.

First, we make sure that buildings are built on the board (lines 64-68).

Next, we want to make sure that buildings do not overlap. Building size is set to one means that it spans half a field each direction (so it forms a square of size one). To avoid floating point errors we multiply differences by two and not divide. This is in lines 70-83.

Next, we make sure that buildings are constructed in some order (line 86).

Finally, we calculate the cost in lines 88-101

Sample output:

Tried aggregator 2 times.
MIP Presolve eliminated 545 rows and 105 columns.
MIP Presolve modified 1807 coefficients.
Aggregator did 456 substitutions.
Reduced MIP has 1050 rows, 655 columns, and 3540 nonzeros.
Reduced MIP has 560 binaries, 95 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.02 sec. (7.63 ticks)
Probing fixed 160 vars, tightened 80 bounds.
Probing changed sense of 40 constraints.
Probing time = 0.00 sec. (6.62 ticks)
Tried aggregator 2 times.
MIP Presolve eliminated 320 rows and 160 columns.
MIP Presolve modified 376 coefficients.
Aggregator did 200 substitutions.
Reduced MIP has 530 rows, 295 columns, and 1860 nonzeros.
Reduced MIP has 200 binaries, 95 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.01 sec. (2.63 ticks)
Probing time = 0.00 sec. (0.53 ticks)
Tried aggregator 1 time.
Reduced MIP has 530 rows, 295 columns, and 1860 nonzeros.
Reduced MIP has 200 binaries, 95 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.85 ticks)
Probing changed sense of 10 constraints.
Probing time = 0.00 sec. (0.60 ticks)
Clique table members: 310.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 8 threads.
Root relaxation solution time = 0.00 sec. (1.22 ticks)

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

      0     0       55.0000    52                     55.0000       59
      0     0       49.0000    80                    Cuts: 88      170
      0     0       49.0000    64                   Cuts: 104      217
      0     0       49.0000    64                  MIRcuts: 1      218
*     0+    0                           39.0000       49.0000      218   25.64%
      0     2       49.0000    64       39.0000       49.0000      218   25.64%
Elapsed time = 0.17 sec. (59.33 ticks, tree = 0.01 MB, solutions = 1)
*  1155   462      integral     0       40.0000       49.0000    16771   22.50%
*  4782+ 1764                           43.0000       49.0000    73723   13.95%
   5076  1836       48.8125    37       43.0000       49.0000    77941   13.95%
*  5190+  648                           45.0000       49.0000    80331    8.89%
   5623   126       47.9286    40       45.0000       49.0000    90495    8.89%

Clique cuts applied:  10
Cover cuts applied:  22
Implied bound cuts applied:  10
Flow cuts applied:  3
Mixed integer rounding cuts applied:  63
Gomory fractional cuts applied:  12

Root node processing (before b&c):
  Real time             =    0.16 sec. (59.06 ticks)
Parallel b&c, 8 threads:
  Real time             =    2.26 sec. (751.26 ticks)
  Sync time (average)   =    0.33 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    2.42 sec. (810.32 ticks)
Result: 45
      55
      55
        33
        33
    1

22
22
  44
  44

Tried aggregator 2 times.

MIP Presolve eliminated 545 rows and 105 columns.

MIP Presolve modified 1807 coefficients.

Aggregator did 456 substitutions.

Reduced MIP has 1050 rows, 655 columns, and 3540 nonzeros.

Reduced MIP has 560 binaries, 95 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.02 sec. (7.63 ticks)

Probing fixed 160 vars, tightened 80 bounds.

Probing changed sense of 40 constraints.

Probing time = 0.00 sec. (6.62 ticks)

Tried aggregator 2 times.

MIP Presolve eliminated 320 rows and 160 columns.

MIP Presolve modified 376 coefficients.

Aggregator did 200 substitutions.

Reduced MIP has 530 rows, 295 columns, and 1860 nonzeros.

Reduced MIP has 200 binaries, 95 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.01 sec. (2.63 ticks)

Probing time = 0.00 sec. (0.53 ticks)

Tried aggregator 1 time.

Reduced MIP has 530 rows, 295 columns, and 1860 nonzeros.

Reduced MIP has 200 binaries, 95 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.00 sec. (0.85 ticks)

Probing changed sense of 10 constraints.

Probing time = 0.00 sec. (0.60 ticks)

Clique table members: 310.

MIP emphasis: balance optimality and feasibility.

MIP search method: dynamic search.

Parallel mode: deterministic, using up to 8 threads.

Root relaxation solution time = 0.00 sec. (1.22 ticks)

Nodes Cuts/

Node Left Objective IInf Best Integer Best Bound ItCnt Gap

0 0 55.0000 52 55.0000 59

0 0 49.0000 80 Cuts: 88 170

0 0 49.0000 64 Cuts: 104 217

0 0 49.0000 64 MIRcuts: 1 218

* 0+ 0 39.0000 49.0000 218 25.64%

0 2 49.0000 64 39.0000 49.0000 218 25.64%

Elapsed time = 0.17 sec. (59.33 ticks, tree = 0.01 MB, solutions = 1)

* 1155 462 integral 0 40.0000 49.0000 16771 22.50%

* 4782+ 1764 43.0000 49.0000 73723 13.95%

5076 1836 48.8125 37 43.0000 49.0000 77941 13.95%

* 5190+ 648 45.0000 49.0000 80331 8.89%

5623 126 47.9286 40 45.0000 49.0000 90495 8.89%

Clique cuts applied: 10

Cover cuts applied: 22

Implied bound cuts applied: 10

Flow cuts applied: 3

Mixed integer rounding cuts applied: 63

Gomory fractional cuts applied: 12

Root node processing (before b&c):

Real time = 0.16 sec. (59.06 ticks)

Parallel b&c, 8 threads:

Real time = 2.26 sec. (751.26 ticks)

Sync time (average) = 0.33 sec.

Wait time (average) = 0.00 sec.

------------

Total (root+branch&cut) = 2.42 sec. (810.32 ticks)

Result: 45