# Land Lot Optimization

by Philippe Olivier
Created 2023-10-30

A couple months back I was asked about land lot optimization. How does a contractor maximize the yield of a lot? It turns out that the interplay between land use regulations (at the municipal/state/federal levels), incentive programs, and environmental laws (just to name a few) make answering this question fairly difficult.

The problem can be more specifically framed as follows: Given a lot of a certain size, and considering a set of laws and regulations, what is the maximum yield (living space) that can be achieved? We will illustrate this with a toy problem:

We have a lot of size 60×40. We want to place up to 5 residential buildings (blue), up to 2 parking lots (grey), and 1 park (green). There are two 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 residential building. The combined area of the parking lots must be at least 10% of the combined area of the residential buildings. We want to generate a 2D plan that maximizes the yield of this lot, which is the combined area of the residential buildings.

At first, it might seem like modeling this would be complicated, especially since the expected output of the model is a 2D plan integrating the items mentioned above, and following the stated rules. However, constraint programming can make use of so-called interval variables that allow us to directly model spacial constraints. We will use CP-SAT, from the OR-Tools optimization suite, to model this problem.

First off, an interval variable is a special kind of variable that is made of three "sub-variables": A starting value `start`, a size `size`, and an ending value `end`. The interval variable enforces that `start + size == end`. All of `start`, `size`, and `end` are (or rather, can be) variables themselves, so it is easy to see that interval variables are powerful and expressive constructs.

We start by importing the necessary packages for optimization and for displaying the solution of the problem.

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

model = cp_model.CpModel()
```

We now define some constants based on the problem definition.

```SIZE_X = 60
SIZE_Y = 40
NUM_BUILDINGS = 5
NUM_PARKING_LOTS = 2

FLOODABLE_X_START = 10
FLOODABLE_X_SIZE = 7
FLOODABLE_X_END = FLOODABLE_X_START + FLOODABLE_X_SIZE
FLOODABLE_Y_START = 20
FLOODABLE_Y_SIZE = 12
FLOODABLE_Y_END = FLOODABLE_Y_START + FLOODABLE_Y_SIZE

UTILITY_X_START = 40
UTILITY_X_SIZE = 5
UTILITY_X_END = UTILITY_X_START + UTILITY_X_SIZE
UTILITY_Y_START = 30
UTILITY_Y_SIZE = 5
UTILITY_Y_END = UTILITY_Y_START + UTILITY_Y_SIZE
```

Before moving on to defining the model, we create a function `entity_2d()` to generate 2D entities. An entity encapsulates a bunch of useful features: its X and Y positions, sizes, and intervals, as well as its area.

```def entity_2d(model, name):
"""Return a 2D entity.

{"X": {"Start": [X starting position],
"Size": [length of X],
"End": [X ending position],
"Interval": [interval variable associated with Start, Size, and End]},
"Y": {"Start": [Y starting position],
"Size": [length of Y],
"End": [Y ending position],
"Interval": [interval variable associated with Start, Size, and End]},
"Area": [area occupied by the entity]}
"""
entity = {"X": {"Start": model.NewIntVar(0, SIZE_X, f"{name}_x_start"),
"Size": model.NewIntVar(0, SIZE_X, f"{name}_x_size"),
"End": model.NewIntVar(0, SIZE_X, f"{name}_x_end")},
"Y": {"Start": model.NewIntVar(0, SIZE_Y, f"{name}_y_start"),
"Size": model.NewIntVar(0, SIZE_Y, f"{name}_y_duration"),
"End": model.NewIntVar(0, SIZE_Y, f"{name}_y_end")}}

entity["X"]["Interval"] = model.NewIntervalVar(entity["X"]["Start"],
entity["X"]["Size"],
entity["X"]["End"],
f"{name}_x_interval")

entity["Y"]["Interval"] = model.NewIntervalVar(entity["Y"]["Start"],
entity["Y"]["Size"],
entity["Y"]["End"],
f"{name}_y_interval")

entity["Area"] = model.NewIntVar(0, SIZE_X*SIZE_Y, f"{name}_area")

# Enforce the size of the area
entity["Y"]["Size"]])

return entity
```

Now we can start defining the model. First, we create entities for the residential buildings, the parking lots, and the park.

```# Building variables
buildings = {i: entity_2d(model, f"building_{i}")
for i in range(NUM_BUILDINGS)}

# Symmetry breaking for buildings
for i in range(NUM_BUILDINGS-1):

# Parking lots variables
parking_lots = {i: entity_2d(model, f"parking_lots_{i}")
for i in range(NUM_PARKING_LOTS)}

# Symmetry breaking for parking lots

# Park variables
park = entity_2d(model, "park")
```

Notice the symmetry breaking constraints. They reduce the size of the solution space by ensuring that similar entities are (mostly) not interchangeable. We define the floodable area and the area occupied by the utility pole manually since their positions are already known.

```# Floodable interval variable
floodable_interval = {"X": model.NewIntervalVar(FLOODABLE_X_START,
FLOODABLE_X_SIZE,
FLOODABLE_X_END,
"floodable_interval_x"),
"Y": model.NewIntervalVar(FLOODABLE_Y_START,
FLOODABLE_Y_SIZE,
FLOODABLE_Y_END,
"floodable_interval_y")}

# Utility pole interval variable
utility_interval = {"X": model.NewIntervalVar(UTILITY_X_START,
UTILITY_X_SIZE,
UTILITY_X_END,
"utility_interval_x"),
"Y": model.NewIntervalVar(UTILITY_Y_START,
UTILITY_Y_SIZE,
UTILITY_Y_END,
"utility_interval_y")}
```

We want to ensure that the entities do not overlap because, well, we can't build a building inside another building. CP-SAT has a very useful constraint that can easily take care of this: AddNoOverlap2D().

```# The buildings, the parking lots, the park, and the floodable area cannot overlap
model.AddNoOverlap2D([buildings[i]["X"]["Interval"] for i in range(NUM_BUILDINGS)] +
[parking_lots[i]["X"]["Interval"] for i in range(NUM_PARKING_LOTS)] +
[park["X"]["Interval"]] +
[floodable_interval["X"]],
[buildings[i]["Y"]["Interval"] for i in range(NUM_BUILDINGS)] +
[parking_lots[i]["Y"]["Interval"] for i in range(NUM_PARKING_LOTS)] +
[park["Y"]["Interval"]] +
[floodable_interval["Y"]])

# The utility pole cannot overlap with the buildings
for i in range(NUM_BUILDINGS):
[buildings[i]["Y"]["Interval"], utility_interval["Y"]])

# The utility pole cannot overlap with the park
[park["Y"]["Interval"], utility_interval["Y"]])
```

We now enforce the park constraints. The park must be at least 10% of the combined building areas and must be at least as large as the largest building.

```# The combined areas of the parkings must be at least 10% the combined areas of the buildings
model.Add(cp_model.LinearExpr.Sum([parking_lots[i]["Area"] for i in range(NUM_PARKING_LOTS)]) * 10 >=
cp_model.LinearExpr.Sum([buildings[i]["Area"] for i in range(NUM_BUILDINGS)]))

# 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, [buildings[i]["Area"] for i in range(NUM_BUILDINGS)])
```

Finally, our objective is to maximize the yield of the lot.

```# The objective is to maximize the yield (building area)
lot_yield = model.NewIntVar(0, SIZE_X*SIZE_Y, "lot_yield")
model.Add(lot_yield == cp_model.LinearExpr.Sum([buildings[i]["Area"] for i in range(NUM_BUILDINGS)]))
model.Maximize(lot_yield)
```

We give the solver only 5 seconds to solve the problem. As we will see later, this is largely enough. We display the results with `matplotlib`.

```# 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")

print(f"Yield: {solver.Value(lot_yield)}")

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

# Residential buildings
for i in range(NUM_BUILDINGS):
solver.Value(buildings[i]["Y"]["Start"])),
solver.Value(buildings[i]["X"]["Size"]),
solver.Value(buildings[i]["Y"]["Size"]),
edgecolor='black',
facecolor='royalblue',
fill=True,
lw=3))

# Parking lots
for i in range(NUM_PARKING_LOTS):
solver.Value(parking_lots[i]["Y"]["Start"])),
solver.Value(parking_lots[i]["X"]["Size"]),
solver.Value(parking_lots[i]["Y"]["Size"]),
edgecolor='black',
facecolor='grey',
fill=True,
lw=3))

# Park
solver.Value(park["Y"]["Start"])),
solver.Value(park["X"]["Size"]),
solver.Value(park["Y"]["Size"]),
edgecolor='black',
facecolor='limegreen',
fill=True,
lw=3))

# Floodable area
FLOODABLE_Y_START),
FLOODABLE_X_SIZE,
FLOODABLE_Y_SIZE,
edgecolor='black',
facecolor='red',
fill=True,
lw=3))

# Utility pole
UTILITY_Y_START),
UTILITY_X_SIZE,
UTILITY_Y_SIZE,
edgecolor='black',
facecolor='orange',
fill=True,
lw=3))

plt.show()
```

Below are some sample results of the model.

The yields of the above lots hover around 1740. The theoretical maximum yield for this problem (disregarding any spatial constraints) is 1781, so our solutions are pretty good. The complete source code of the model can be found here.

This website is generated by Emacs.