Islanders.cs

// 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();
}










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