MILP formulation of Islanders game logic

milp_pulp.py

#!python
# -*- coding: utf-8 -*-
"""
Module milp_pulp.py: MILP solution of Islanders building placing problem

Created on 2019/04/23

author: George Zogopoulos
"""

import pulp
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

prob = pulp.LpProblem('Islanders', pulp.LpMaximize)

# =============================================================================
# Helper structures
# =============================================================================


class Building:
    def __init__(self, name, x_c, y_c, dim, b_range):
        self.name = name
        self.x_coord = x_c
        self.y_coord = y_c
        self.dimension = dim
        self.b_range = b_range
        self.order = None

    def set_order(self, order):
        self.order = order

    def plot(self, axh):
        # Draw building range
#        aoe = mpl.patches.Circle((self.x_coord, self.y_coord),
#                                 radius=self.b_range, alpha=0.3, color='green')
        diag_range = self.b_range*1.4142
        aoe = mpl.patches.Rectangle((self.x_coord, self.y_coord-self.b_range),
                                    diag_range, diag_range,
                                    facecolor='none', alpha=10,
                                    edgecolor='black', linewidth=1,
                                    linestyle='--', angle=45)
        axh.add_artist(aoe)
        # Draw building
        rect = mpl.patches.Rectangle((self.x_coord-self.dimension/2, self.y_coord-self.dimension/2),
                                     self.dimension, self.dimension,
                                     facecolor='cyan', edgecolor='black', linewidth=3)
        axh.add_artist(rect)
        # Write building name
        axh.annotate(f'{self.name}:{self.order}', (self.x_coord, self.y_coord),
                     horizontalalignment='center', verticalalignment='center')


# =============================================================================
# Input data
# =============================================================================

# Building base score array; City Center, House, Mansion, Fountain
base_score = np.array([15, 1, 1, 0])

# Building proximity score matrix
prox_score = np.array([
    [-40, 0, 0, 0],
    [6, 1, 0, 2],
    [8, 0, 1, 3],
    [7, 2, 3, -15],
    ])

# Building dimension
dim = np.array([1, 2, 2, 2])

# Building range
build_range = np.array([5, 6, 3, 8])

# =============================================================================
# Define number of buildings
# =============================================================================

#num_buildings = np.array([1, 20, 8, 1])
num_buildings = np.array([1, 4, 0, 0])

# =============================================================================
# Build other parameters
# =============================================================================

T = num_buildings.sum()  # Maximum time
L = 10  # Maximum grid dimension
eps = 0.01  # Epsilon value

# =============================================================================
# Build decision variables
# =============================================================================

city_center_list = ['C'+str(i) for i in range(num_buildings[0])]
houses_list = ['H'+str(i) for i in range(num_buildings[1])]
mansions_list = ['M'+str(i) for i in range(num_buildings[2])]
fountains_list = ['F'+str(i) for i in range(num_buildings[3])]
buildings_list = city_center_list + houses_list + mansions_list + fountains_list
build_range_list = [build_range[0]]*len(city_center_list) + \
                    [build_range[1]]*len(houses_list) + \
                    [build_range[2]]*len(mansions_list) + \
                    [build_range[3]]*len(fountains_list)

dim_list = [dim[0]]*len(city_center_list) + \
            [dim[1]]*len(houses_list) + \
            [dim[2]]*len(mansions_list) + \
            [dim[3]]*len(fountains_list)

ranges_dict = dict(zip(buildings_list, build_range_list))
dimension_dict = dict(zip(buildings_list, dim_list))

# Build primary variables
N = num_buildings.sum()
px = pulp.LpVariable.dicts('x_coord', (buildings_list), -L, L, pulp.LpContinuous)  # North positive
py = pulp.LpVariable.dicts('y_coord', (buildings_list), -L, L, pulp.LpContinuous)  # East positive
t = pulp.LpVariable.dicts('build_time', (buildings_list), 0, T-1, pulp.LpInteger)

# Build auxiliary variables
placed_before = pulp.LpVariable.dicts("placed_before", (buildings_list, buildings_list),
                                      0, 1, pulp.LpInteger)
east_of = pulp.LpVariable.dicts("east_of", (buildings_list, buildings_list),
                                      0, 1, pulp.LpInteger)
north_of = pulp.LpVariable.dicts("north_of", (buildings_list, buildings_list),
                                      0, 1, pulp.LpInteger)
covers = pulp.LpVariable.dicts("covers", (buildings_list, buildings_list),
                                 0, 1, pulp.LpInteger)  # b1 covers b2 when placed
x_overlaps = pulp.LpVariable.dicts("x_overlaps", (buildings_list, buildings_list),
                                 0, 1, pulp.LpInteger)
y_overlaps = pulp.LpVariable.dicts("y_overlaps", (buildings_list, buildings_list),
                                 0, 1, pulp.LpInteger)
gives_score = pulp.LpVariable.dicts("gives_score", (buildings_list, buildings_list),
                                    0, 1, pulp.LpInteger)  # b1 gives score to b2 (b1 placed before b2))

# =============================================================================
# Build Constraints
# =============================================================================

# Build order
for idx_1 in range(N):
    for idx_2 in np.arange(idx_1+1, N):
        b1 = buildings_list[idx_1]
        b2 = buildings_list[idx_2]
        prob += t[b2] - t[b1] + 1 <= placed_before[b1][b2]*T, ""
        prob += t[b1] - t[b2] + 1 <= placed_before[b2][b1]*T, ""
        prob += placed_before[b1][b2] + placed_before[b2][b1] == 1, ""

# Build relative direction auxiliary variables
for idx_1 in range(N):
    for idx_2 in np.arange(idx_1+1, N):
        b1 = buildings_list[idx_1]
        b2 = buildings_list[idx_2]
        # Find if b1 is east of b2
        prob += px[b2] - px[b1] >= -east_of[b1][b2]*L, ""
        prob += px[b1] - px[b2] >= -east_of[b2][b1]*L, ""
        prob += east_of[b2][b1] + east_of[b1][b2] == 1, ""
        # Find if b1 is north of b2
        prob += py[b2] - py[b1] >= -north_of[b1][b2]*L, ""
        prob += py[b1] - py[b2] >= -north_of[b2][b1]*L, ""
        prob += north_of[b2][b1] + north_of[b1][b2] == 1, ""

# Build overlapping auxiliary variables
for idx_1 in range(N):
    for idx_2 in np.arange(idx_1+1, N):
        b1 = buildings_list[idx_1]
        b2 = buildings_list[idx_2]
        # East-West direction
        # Case: b2 is east of b1
        prob += px[b2] - px[b1] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) + eps - (x_overlaps[b1][b2]+east_of[b1][b2])*L, ""  # case is true, x_overlaps not needed
        prob += px[b2] - px[b1] <= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - eps + (1-x_overlaps[b1][b2]+east_of[b1][b2])*L, ""  # case is false, x_overlaps needed 0
        # Case: b2 is west of b1
        prob += px[b1] - px[b2] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) + eps - (x_overlaps[b1][b2]+east_of[b2][b1])*L, ""  # case is true, x_overlaps not needed
        prob += px[b1] - px[b2] <= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - eps + (1-x_overlaps[b1][b2]+east_of[b2][b1])*L, ""  # case is false, x_overlaps needed 0
        prob += x_overlaps[b1][b2] == x_overlaps[b2][b1], ""
        # North-South direction
        # Case: b2 is north of b1
        prob += py[b2] - py[b1] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) + eps - (y_overlaps[b1][b2]+north_of[b1][b2])*L, ""  # case is true, y_overlaps not needed
        prob += py[b2] - py[b1] <= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - eps + (1-y_overlaps[b1][b2]+north_of[b1][b2])*L, ""  # case is false, y_overlaps needed 0
        # Case: b2 is south of b1
        prob += py[b1] - py[b2] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) + eps - (y_overlaps[b1][b2]+north_of[b2][b1])*L, ""  # case is true, y_overlaps not needed
        prob += py[b1] - py[b2] <= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - eps + (1-y_overlaps[b1][b2]+north_of[b2][b1])*L, ""  # case is false, y_overlaps needed 0
        prob += y_overlaps[b1][b2] == y_overlaps[b2][b1], ""

# Build covers auxiliary variable
# Using the rhombus shape (Manhattan distance) to calculate coverage
for idx_1 in range(N):
    for idx_2 in np.arange(idx_1+1, N):
        # Two passes are needed because convers[b1][b2] is independent of convers[b2][b1]
        # First pass to build covers[b1][b2]
        b1 = buildings_list[idx_1]
        b2 = buildings_list[idx_2]
        # b1 is east and north of b2
        prob += px[b1]-px[b2]+py[b1]-py[b2] <= \
            ranges_dict[b1] + (east_of[b2][b1]+north_of[b2][b1]+1-covers[b1][b2])*L, ""
        # b1 is west and north of b2
        prob += px[b2]-px[b1]+py[b1]-py[b2] <= \
            ranges_dict[b1] + (east_of[b1][b2]+north_of[b2][b1]+1-covers[b1][b2])*L, ""
        # b1 is east and south of b2
        prob += px[b1]-px[b2]+py[b2]-py[b1] <= \
            ranges_dict[b1] + (east_of[b2][b1]+north_of[b1][b2]+1-covers[b1][b2])*L, ""
        # b1 is west and south of b2
        prob += px[b2]-px[b1]+py[b2]-py[b1] <= \
            ranges_dict[b1] + (east_of[b1][b2]+north_of[b1][b2]+1-covers[b1][b2])*L, ""
        # Second pass to build covers[b2][b1]
        b1 = buildings_list[idx_2]
        b2 = buildings_list[idx_1]
        # b1 is east and north of b2
        prob += px[b1]-px[b2]+py[b1]-py[b2] <= \
            ranges_dict[b1] + (east_of[b2][b1]+north_of[b2][b1]+1-covers[b1][b2])*L, ""
        # b1 is west and north of b2
        prob += px[b2]-px[b1]+py[b1]-py[b2] <= \
            ranges_dict[b1] + (east_of[b1][b2]+north_of[b2][b1]+1-covers[b1][b2])*L, ""
        # b1 is east and south of b2
        prob += px[b1]-px[b2]+py[b2]-py[b1] <= \
            ranges_dict[b1] + (east_of[b2][b1]+north_of[b1][b2]+1-covers[b1][b2])*L, ""
        # b1 is west and south of b2
        prob += px[b2]-px[b1]+py[b2]-py[b1] <= \
            ranges_dict[b1] + (east_of[b1][b2]+north_of[b1][b2]+1-covers[b1][b2])*L, ""


# Enforce dimension separation if other dimension is close
for idx_1 in range(N):
    for idx_2 in np.arange(idx_1+1, N):
        b1 = buildings_list[idx_1]
        b2 = buildings_list[idx_2]
        # b1 is east of b2
        prob += px[b1]-px[b2] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - (y_overlaps[b1][b2]+east_of[b2][b1])*L, ""
        # b2 is east of b1
        prob += px[b2]-px[b1] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - (y_overlaps[b1][b2]+east_of[b1][b2])*L, ""
        # b1 is north of b2
        prob += py[b1]-py[b2] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - (x_overlaps[b1][b2]+north_of[b2][b1])*L, ""
        # b2 is north of b1
        prob += py[b2]-py[b1] >= 0.5*(dimension_dict[b1]+dimension_dict[b2]) - (x_overlaps[b1][b2]+north_of[b1][b2])*L, ""
        prob += x_overlaps[b1][b2] + y_overlaps[b1][b2] <= 1, ""  # Simultaneous overlap restricted
        prob += x_overlaps[b1][b2] == x_overlaps[b2][b1], ""
        prob += y_overlaps[b1][b2] == y_overlaps[b2][b1], ""

# Score counting
for idx_1 in range(N):
    for idx_2 in np.arange(idx_1+1, N):
        b1 = buildings_list[idx_1]
        b2 = buildings_list[idx_2]
        # b1 is placed befor b2
        prob += gives_score[b1][b2] <= placed_before[b1][b2], ""
        prob += gives_score[b1][b2] <= covers[b2][b1], ""
        # b2 is placed before b1
        prob += gives_score[b2][b1] <= placed_before[b2][b1], ""
        prob += gives_score[b2][b1] <= covers[b1][b2], ""

# =============================================================================
# Cost function
# =============================================================================
#prob += pulp.lpSum(covers)
prob += pulp.lpSum(gives_score)
# =============================================================================
# Run solver
# =============================================================================
prob.writeLP("Islanders.lp")
prob.solve()
print("Status:", pulp.LpStatus[prob.status])

# =============================================================================
# Display results
# =============================================================================

for b1 in buildings_list:
    print(f'Building name: {b1}')
    print(f'Building coordinates: ({px[b1].value()}, {py[b1].value()})')
    print(f'Build order: {t[b1].value()}')
    for b2 in buildings_list:
        if (b1 != b2):
            if east_of[b1][b2].value():
                print(f'{b1} is east of {b2}')
            else:
                print(f'{b1} is west of {b2}')
            if north_of[b1][b2].value():
                print(f'{b1} is north of {b2}')
            else:
                print(f'{b1} is south of {b2}')
            print('---Build Order---:')
            print(f'{b1} is placed before {b2}: {bool(placed_before[b1][b2].value())}')
            print(f'{b1} is not placed before {b2}: {bool(placed_before[b2][b1].value())}')
            print('---Ranges---:')
            if covers[b1][b2].value():
                print(f'{b1} covers {b2}')
            else:
                print(f'{b1} does not cover {b2}')
            print('---Scoring---:')
            print(f'{b1} gives score to {b2}: {bool(gives_score[b1][b2].value())}')
            print()

# Check constraints for violation
for c_name, c in prob.constraints.items():
    value = c.value()
    c_type = c.sense
    if (c_type==0) and (value!=0):
        print(f'Constraint {c_name}:{c} is {value} and not zero')
    if (c_type==1) and (value<0):
        print(f'Constraint {c_name}:{c} is {value} and not >=0')
    if (c_type==-1) and (value>0):
        print(f'Constraint {c_name}:{c} is {value} and not <=0')


buildings = []
for b in buildings_list:
    building = Building(b, px[b].value(), py[b].value(),
                        dimension_dict[b], ranges_dict[b])
    building.set_order(t[b].value())
    buildings.append(building)

# Create a figure
fig = plt.figure(figsize=(10,10))
fig.suptitle('Building Plan')
axh = fig.add_axes([0.1, 0.1, 0.8, 0.8], aspect='equal')
axh.set_xlim([-L-5, L+5])
axh.set_ylim([-L-5, L+5])
axh.set_xticks(np.linspace(-L, L, 11))
axh.set_yticks(np.linspace(-L, L, 11))
axh.grid()

for b in buildings:
    b.plot(axh)

plt.show()