Constraint-guided 2D spatial generator.

lot-generator.py

"""Constraint-guided 2D spatial generator.

We have a lot of size 60x40. We want to place up to 5 residential buildings (blue), up to 2 parking
lots (grey), and 1 park (green). There are 2 problematic areas on the lot. The first one (red) is a
floodable area, and nothing can be built there. The second one (orange) is a utility pole: we cannot
build a residential building there nor have it in the park, but we are allowed to build a parking
lot around it. The size of the park must be at least as large as the largest of the residential
buildings. The combined area of the parkings must be at least 10% of the combined area of the
residential buildings. We want to maximize the yield of the lot, which is the combined area of the
residential buildings.
"""


import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from ortools.sat.python import cp_model


SIZE_X = 60
SIZE_Y = 40
NUM_BUILDINGS = 5
NUM_PARKINGS = 2

FLOODABLE_X_S = 10
FLOODABLE_X_D = 7
FLOODABLE_X_E = FLOODABLE_X_S + FLOODABLE_X_D
FLOODABLE_Y_S = 20
FLOODABLE_Y_D = 12
FLOODABLE_Y_E = FLOODABLE_Y_S + FLOODABLE_Y_D

UTILITY_X_S = 40
UTILITY_X_D = 5
UTILITY_X_E = UTILITY_X_S + UTILITY_X_D
UTILITY_Y_S = 30
UTILITY_Y_D = 5
UTILITY_Y_E = UTILITY_Y_S + UTILITY_Y_D

model = cp_model.CpModel()

# Building variables
buildings = {i: {"x": {"s": model.NewIntVar(0, SIZE_X, f"building_{i}_x_start"),
                       "d": model.NewIntVar(0, SIZE_X, f"building_{i}_x_duration"),
                       "e": model.NewIntVar(0, SIZE_X, f"building_{i}_x_end")},
                 "y": {"s": model.NewIntVar(0, SIZE_Y, f"building_{i}_y_start"),
                       "d": model.NewIntVar(0, SIZE_Y, f"building_{i}_y_duration"),
                       "e": model.NewIntVar(0, SIZE_Y, f"building_{i}_y_end")}}
             for i in range(NUM_BUILDINGS)}

buildings_intervals = {i: {"x": model.NewIntervalVar(buildings[i]["x"]["s"],
                                                     buildings[i]["x"]["d"],
                                                     buildings[i]["x"]["e"],
                                                     f"building_{i}_interval_x"),
                           "y": model.NewIntervalVar(buildings[i]["y"]["s"],
                                                     buildings[i]["y"]["d"],
                                                     buildings[i]["y"]["e"],
                                                     f"building_{i}_interval_y")}
                       for i in range(NUM_BUILDINGS)}

building_areas = {i: model.NewIntVar(0, SIZE_X*SIZE_Y, f"building_{i}_area")
                  for i in range(NUM_BUILDINGS)}

for i in range(NUM_BUILDINGS):
    model.AddMultiplicationEquality(building_areas[i], [buildings[i]["x"]["d"],
                                                        buildings[i]["y"]["d"]])

# Parking variables
parkings = {i: {"x": {"s": model.NewIntVar(0, SIZE_X, f"parking_{i}_x_start"),
                      "d": model.NewIntVar(0, SIZE_X, f"parking_{i}_x_duration"),
                      "e": model.NewIntVar(0, SIZE_X, f"parking_{i}_x_end")},
                "y": {"s": model.NewIntVar(0, SIZE_Y, f"parking_{i}_y_start"),
                      "d": model.NewIntVar(0, SIZE_Y, f"parking_{i}_y_duration"),
                      "e": model.NewIntVar(0, SIZE_Y, f"parking_{i}_y_end")}}
            for i in range(NUM_PARKINGS)}

parkings_intervals = {i: {"x": model.NewIntervalVar(parkings[i]["x"]["s"],
                                                    parkings[i]["x"]["d"],
                                                    parkings[i]["x"]["e"],
                                                    f"parking_{i}_interval_x"),
                          "y": model.NewIntervalVar(parkings[i]["y"]["s"],
                                                    parkings[i]["y"]["d"],
                                                    parkings[i]["y"]["e"],
                                                    f"parking_{i}_interval_y")}
                      for i in range(NUM_PARKINGS)}

parking_areas = {i: model.NewIntVar(0, SIZE_X*SIZE_Y, f"parking_{i}_area")
                 for i in range(NUM_PARKINGS)}

for i in range(NUM_PARKINGS):
    model.AddMultiplicationEquality(parking_areas[i], [parkings[i]["x"]["d"],
                                                       parkings[i]["y"]["d"]])

# Park variables
park = {"x": {"s": model.NewIntVar(0, SIZE_X, "park_x_start"),
              "d": model.NewIntVar(0, SIZE_X, "park_x_duration"),
              "e": model.NewIntVar(0, SIZE_X, "park_x_end")},
        "y": {"s": model.NewIntVar(0, SIZE_Y, "park_y_start"),
              "d": model.NewIntVar(0, SIZE_Y, "park_y_duration"),
              "e": model.NewIntVar(0, SIZE_Y, "park_y_end")}}

park_interval = {"x": model.NewIntervalVar(park["x"]["s"],
                                           park["x"]["d"],
                                           park["x"]["e"],
                                           "park_interval_x"),
                 "y": model.NewIntervalVar(park["y"]["s"],
                                           park["y"]["d"],
                                           park["y"]["e"],
                                           "park_interval_y")}

park_area = model.NewIntVar(0, SIZE_X*SIZE_Y, f"park_area")

model.AddMultiplicationEquality(park_area, [park["x"]["d"],
                                            park["y"]["d"]])

# Floodable interval
floodable_interval = {"x": model.NewIntervalVar(FLOODABLE_X_S,
                                                FLOODABLE_X_D,
                                                FLOODABLE_X_E,
                                                "floodable_interval_x"),
                      "y": model.NewIntervalVar(FLOODABLE_Y_S,
                                                FLOODABLE_Y_D,
                                                FLOODABLE_Y_E,
                                                "floodable_interval_y")}

# Utility pole interval
utility_interval = {"x": model.NewIntervalVar(UTILITY_X_S,
                                              UTILITY_X_D,
                                              UTILITY_X_E,
                                              "utility_interval_x"),
                    "y": model.NewIntervalVar(UTILITY_Y_S,
                                              UTILITY_Y_D,
                                              UTILITY_Y_E,
                                              "utility_interval_y")}

# The buildings, the parking lots, the park, and the floodable area cannot overlap
model.AddNoOverlap2D([buildings_intervals[i]["x"] for i in range(NUM_BUILDINGS)] +
                     [parkings_intervals[i]["x"] for i in range(NUM_PARKINGS)] +
                     [park_interval["x"]] +
                     [floodable_interval["x"]],
                     [buildings_intervals[i]["y"] for i in range(NUM_BUILDINGS)] +
                     [parkings_intervals[i]["y"] for i in range(NUM_PARKINGS)] +
                     [park_interval["y"]] +
                     [floodable_interval["y"]])

# The utility pole cannot overlap with the buildings and the park
for i in range(NUM_BUILDINGS):
    model.AddNoOverlap2D([buildings_intervals[i]["x"], utility_interval["x"]],
                         [buildings_intervals[i]["y"], utility_interval["y"]])

model.AddNoOverlap2D([park_interval["x"], utility_interval["x"]],
                     [park_interval["y"], utility_interval["y"]])

# The combined areas of the parkings must be at least 10% the combined areas of the buildings
model.Add(cp_model.LinearExpr.Sum(parking_areas.values()) * 10 >=
          cp_model.LinearExpr.Sum(building_areas.values()))

# The area of the park must be at least as large as the area of the largest building
largest_building = model.NewIntVar(0, SIZE_X*SIZE_Y, "largest_building")
model.AddMaxEquality(largest_building, building_areas.values())
model.Add(park_area >= largest_building)

# The objective is to maximize the yield (building area)
model.Maximize(cp_model.LinearExpr.Sum(building_areas.values()))

# Solve the problem with a time limit
solver = cp_model.CpSolver()
solver.parameters.max_time_in_seconds = 5
status = solver.Solve(model)
if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
    print("Model neither optimal nor feasible")

# Display results
fig, ax = plt.subplots()
plt.xlim([0, SIZE_X])
plt.ylim([0, SIZE_Y])

# Residential buildings
for i in range(NUM_BUILDINGS):
    ax.add_patch(Rectangle((solver.Value(buildings[i]["x"]["s"]),
                            solver.Value(buildings[i]["y"]["s"])),
                           solver.Value(buildings[i]["x"]["d"]),
                           solver.Value(buildings[i]["y"]["d"]),
                           edgecolor='black',
                           facecolor='royalblue',
                           fill=True,
                           lw=3))

# Parkings
for i in range(NUM_PARKINGS):
    ax.add_patch(Rectangle((solver.Value(parkings[i]["x"]["s"]),
                            solver.Value(parkings[i]["y"]["s"])),
                           solver.Value(parkings[i]["x"]["d"]),
                           solver.Value(parkings[i]["y"]["d"]),
                           edgecolor='black',
                           facecolor='grey',
                           fill=True,
                           lw=3))

# Park
ax.add_patch(Rectangle((solver.Value(park["x"]["s"]),
                        solver.Value(park["y"]["s"])),
                       solver.Value(park["x"]["d"]),
                       solver.Value(park["y"]["d"]),
                       edgecolor='black',
                       facecolor='limegreen',
                       fill=True,
                       lw=3))

# Floodable area
ax.add_patch(Rectangle((FLOODABLE_X_S,
                        FLOODABLE_Y_S),
                       FLOODABLE_X_D,
                       FLOODABLE_Y_D,
                       edgecolor='black',
                       facecolor='red',
                       fill=True,
                       lw=3))

# Utility pole
ax.add_patch(Rectangle((UTILITY_X_S,
                        UTILITY_Y_S),
                       UTILITY_X_D,
                       UTILITY_Y_D,
                       edgecolor='black',
                       facecolor='orange',
                       fill=True,
                       lw=3))

# plt.savefig("lot.png")
plt.show()