pythonoptimizationor-toolscp-sat
Optimizing node placement in a 2D grid to match specific geodesic distances
I'm working on a problem where I need to arrange a set of nodes within a 2D grid such that the distance between pairs of nodes either approximate specific values as closely as possible or meet a minimum threshold.
In other words, for some node pairs, the distance should be approximately equal to a predefined value (≈). For other node pairs, the distance should be greater than or equal to a predefined threshold (≥).
Additional challenge: the grid is inscribed inside a concave polygon, so distances must be geodesic.
Question: Using OR-Tools, how can I efficiently approximate the location of the nodes given the constraints above-mentioned?
[EDIT] I've revised the example script, trying as best I can to apply Laurent's wise suggestions, but my poor understanding of OR-Tools (and its subtleties) still makes this a very difficult task.
This version:
- creates a simple conforming grid
- pre-computes geodesic distances between all pairs of cells within that grid
- indicates the number of nodes to place as well as their associated target pairwise distances
- each pairwise distance comes with an objective to match (either ≈ or ≥)
- declares the CP-SAT model and creates the main variables for the problem
- (new) ensures each node is assigned exactly to one position and each position can have at most one node
- (new) creates a Boolean variable checking if the distance constraint is met for each pair of nodes
- (new) use
AddImplicationto connect the Boolean variables with the node positions - (new) applies conditional penalties based on whether the distance condition is met and tries to minimize their sum
Unfortunately, I must be missing a few nuances because this implementation doesn't return any results even though the solution space not null.
from ortools.sat.python import cp_model
from itertools import combinations
import networkx as nx
# Dimensions of the original grid (width & height)
w, h = 7, 5
# Selection of grid-cell indices (conforming/concave grid)
cell_indices = list(sorted(set(range(w * h)) - set([0, 1, 2, 3, 7, 8, 9, 10, 28, 29])))
# Topology of the conforming/concave grid
T = nx.Graph()
for i in cell_indices:
if i >= w and (i - w) in cell_indices:
T.add_edge(i, i - w, weight=1)
if i < w * (h - 1) and (i + w) in cell_indices:
T.add_edge(i, i + w, weight=1)
if i % w != 0 and (i - 1) in cell_indices:
T.add_edge(i, i - 1, weight=1)
if (i + 1) % w != 0 and (i + 1) in cell_indices:
T.add_edge(i, i + 1, weight=1)
# Precompute geodesic distances using Dijkstra's algorithm
geodesic_distances = dict(nx.all_pairs_dijkstra_path_length(T))
# Get the largest geodesic distance
max_distance = float('-inf')
for i1 in geodesic_distances:
for i2 in geodesic_distances[i1]:
if i1 != i2 and i1 > i2:
distance = geodesic_distances[i1][i2]
if distance > max_distance: max_distance = distance
# Number of nodes to place
num_nodes = 5
# Target distances to match between each pair of nodes + type of objective (≈ or ≥)
objective_distances = {(0, 1): (3, '≈'),
(0, 2): (2, '≥'),
(0, 3): (2, '≥'),
(0, 4): (3, '≈'),
(1, 2): (3, '≥'),
(1, 3): (3, '≥'),
(1, 4): (4, '≥'),
(2, 3): (2, '≈'),
(2, 4): (4, '≥'),
(3, 4): (3, '≈')}
# Instantiate model
model = cp_model.CpModel()
# Ensure each position can have at most one node
node_at_position = {}
for index in cell_indices:
at_most_one = []
for node in range(num_nodes):
var = model.NewBoolVar(f'node_{node}_at_position_{index}')
node_at_position[node, index] = var
at_most_one.append(var)
# Apply at most one node per position constraint
model.AddAtMostOne(at_most_one)
# Ensure each node is assigned exactly to one position
for node in range(num_nodes):
model.AddExactlyOne(node_at_position[node, idx] for idx in cell_indices)
penalties = []
# For each pair of nodes:
for (node1, node2), (target_distance, constraint_type) in objective_distances.items():
# For each compatible pair of cells
for i1, i2 in combinations(cell_indices, 2):
# Get the corresponding geodesic distance
distance = geodesic_distances[i1][i2]
# Create a Boolean variable
is_compatible = model.NewBoolVar(f'compat_{node1}_{node2}_{i1}_{i2}')
# Create a penalty variable
penalty = model.NewIntVar(0, max_distance, f'penalty_{node1}_{node2}_{i1}_{i2}')
if constraint_type == '≈':
# Condition that `is_compatible` will be True if the distance approximates (deviation: -1/+1) the target distance
model.Add(is_compatible == (target_distance - 1 <= distance <= target_distance + 1))
elif constraint_type == '≥':
# Condition that `is_compatible` will be True if the distance is at least the target distance
model.Add(is_compatible == (distance >= target_distance))
# If 'is_compatible' is true -> implications to enforce node positions
model.AddImplication(is_compatible, node_at_position[node1, i1])
model.AddImplication(is_compatible, node_at_position[node2, i2])
# If it is not -> add a penalty
model.Add(penalty == abs(distance - target_distance)).OnlyEnforceIf(is_compatible.Not())
# Accumulate penalties
penalties.append(penalty)
# Objective to minimize total penalty
model.Minimize(sum(penalties))
# Solving the model
solver = cp_model.CpSolver()
status = solver.Solve(model)
print("Solver status:", solver.StatusName(status))
if status == cp_model.FEASIBLE or status == cp_model.OPTIMAL:
print("Solution found:")
for node in range(num_nodes):
for index in cell_indices:
if solver.Value(node_at_position[node, index]):
print(f'Node {node} is at position {index}')
else:
print("No solution found.")
Solution
If you're not beholden to ortools, here is a solution using pyomo and HiGHS solver. Minimally, this might give you some insights on the correct ortools syntax to implement some of these types of constraints, if this formulation works. I'm not familiar w/ the ortools syntax, but this aligns with the context of what @Laurent Perron mentioned above.
This is a pretty stout BIP (Binary Integer Program) and I've put an absolute gap in the solve such that if it gets within 1.0 of the theoretical limit, it will stop rather than continuing the solve. With that, it pops out an optimal answer in about 10 minutes on my machine, with the solution shown at bottom with a total penalty of 3 distance units for the example shown.
Code:
"""
placing nodes within a grid
Created on: 5/31/24
"""
from collections import defaultdict
import pyomo.environ as pyo
import highspy
import networkx as nx
# Dimensions of the original grid (width & height)
w, h = 7, 5
# Selection of grid-cell indices (conforming/concave grid)
cell_indices = list(sorted(set(range(w * h)) - {0, 1, 2, 3, 7, 8, 9, 10, 28, 29}))
# Topology of the conforming/concave grid
T = nx.Graph()
for i in cell_indices:
if i >= w and (i - w) in cell_indices:
T.add_edge(i, i - w, weight=1)
if i < w * (h - 1) and (i + w) in cell_indices:
T.add_edge(i, i + w, weight=1)
if i % w != 0 and (i - 1) in cell_indices:
T.add_edge(i, i - 1, weight=1)
if (i + 1) % w != 0 and (i + 1) in cell_indices:
T.add_edge(i, i + 1, weight=1)
# Precompute geodesic distances using Dijkstra's algorithm
geodesic_distances = dict(nx.all_pairs_dijkstra_path_length(T))
# Get the largest geodesic distance
max_distance = float('-inf')
for i1 in geodesic_distances:
for i2 in geodesic_distances[i1]:
if i1 != i2 and i1 > i2:
distance = geodesic_distances[i1][i2]
if distance > max_distance: max_distance = distance
# Number of nodes to place
num_nodes = 5
# Target distances to match between each pair of nodes + type of objective (≈ or ≥)
objective_distances = {(0, 1): (3, '≈'),
(0, 2): (2, '≥'),
(0, 3): (2, '≥'),
(0, 4): (3, '≈'),
(1, 2): (3, '≥'),
(1, 3): (3, '≥'),
(1, 4): (4, '≥'),
(2, 3): (2, '≈'),
(2, 4): (4, '≥'),
(3, 4): (3, '≈')
}
# print(geodesic_distances)
# let's chop up the pairs into needed subsets, and hook them up to penalty values
fudge_factor = 1.0 # the max allowable miss for ≈ constraints
cutoff_points = [t[0] for t in objective_distances.values()]
ge_pairs = defaultdict(dict) # GTE pairs at cutoff point
ae_pairs = defaultdict(dict) # Approx Equal pairs at cutoff point
for idx in cell_indices:
for other, dist in geodesic_distances[idx].items():
if other <= idx:
continue # we only need "sorted" pairs like the objective_distances
for c in cutoff_points:
if dist >= c:
ge_pairs[c][idx, other] = dist - c
if c - fudge_factor <= dist <= c + fudge_factor:
ae_pairs[c][idx, other] = abs(dist - c)
# now bin up the requirements as basis for the associated constraints
ge_reqts = {}
ae_reqts = {}
for (n1, n2), (c, c_type) in objective_distances.items():
if c_type == '≈':
ae_reqts[n1, n2] = c
elif c_type == '≥':
ge_reqts[n1, n2] = c
else:
raise ValueError()
# on to the model...
m = pyo.ConcreteModel()
# === SETS ===
m.G = pyo.Set(initialize=cell_indices, doc='Grid Points')
m.GG = pyo.Set(initialize=[(g1, g2) for g1 in m.G for g2 in m.G if g2 > g1], doc=('paired '
'assignment at '
'g_1, g_2)'))
m.N = pyo.Set(initialize=list(range(num_nodes)), doc='Nodes')
m.NN = pyo.Set(initialize=objective_distances.keys(), doc='assignment requirements')
def eligible_assignments(m, *nn):
# a helper function to determine (limit) assignments to eligible spots
if nn in ae_reqts:
cutoff = ae_reqts[nn]
eligible_locations = ae_pairs[cutoff]
return eligible_locations
if nn in ge_reqts:
cutoff = ge_reqts[nn]
eligible_locations = ge_pairs[cutoff]
return eligible_locations
m.eligible_locations = pyo.Set(m.NN, initialize=eligible_assignments)
m.NNGG = pyo.Set(initialize=[(nn, gg) for nn in m.NN for gg in m.eligible_locations[nn]])
# === VARS ===
m.place = pyo.Var(m.N, m.G, domain=pyo.Binary, doc='place node N at point G')
m.paired_assignment = pyo.Var(m.NNGG, domain=pyo.Binary, doc='node pair assigned to grid pair')
# === OBJ ===
# helper expressions
ge_penalties = sum(m.paired_assignment[nn, gg] * ge_pairs[ge_reqts[nn]][gg] for nn in ge_reqts for gg
in m.eligible_locations[nn])
ae_penalties = sum(m.paired_assignment[nn, gg] * ae_pairs[ae_reqts[nn]][gg] for nn in ae_reqts for gg
in m.eligible_locations[nn])
# the OBJECTIVE
m.obj = pyo.Objective(expr=ge_penalties + ae_penalties, sense=pyo.minimize)
# === CONSTRAINTS ===
@m.Constraint(m.N)
def assign_exactly_once(m, n):
return sum(m.place[n, g] for g in m.G) == 1
@m.Constraint(m.G)
def max_one_per_grid_point(m, g):
return sum(m.place[n, g] for n in m.N) <= 1
@m.Constraint(m.NNGG)
def link_assignment(m, n1, n2, *gg):
return m.paired_assignment[n1, n2, gg] <= sum(m.place[n1, g] + m.place[n2, g] for g in gg)/2
@m.Constraint(m.NN)
def assign_pair_once(m, *nn):
return sum(m.paired_assignment[nn, gg] for gg in m.eligible_locations[nn]) == 1
# === QA ===
# m.pprint()
# === SOLVE ===
opt = pyo.SolverFactory('appsi_highs')
res = opt.solve(m, options={'mip_abs_gap': 1.1}, tee=True)
# abs gap implies that if we get within 1.0 of theoretical limit, that's good enough to stop.
if pyo.check_optimal_termination(res):
print('Optimal termination')
print(res)
print(' --- placements ---')
for n in sorted(m.N):
for g in m.G:
if pyo.value(m.place[n, g]) > 0.5:
print(f'place {n} at {g}')
print(' --- pairing penalties ---')
for nngg in m.NNGG:
if pyo.value(m.paired_assignment[nngg]) > 0.5:
nn = nngg[:2]
gg = nngg[2:]
if nn in ge_reqts:
c = ge_reqts.get(nn, None)
penalty = ge_pairs[c][gg]
print(f'assign {nn} to (unordered) {gg} for penalty {penalty}', end='')
else:
c = ae_reqts.get(nn, None)
penalty = ae_pairs[c][gg]
print(f'assign {nn} to (unordered) {gg} for penalty {penalty}', end='')
print(f' (dist: {geodesic_distances[gg[0]][gg[1]]})')
else:
print('failed solve ... see log')
Output:
...
Solving report
Status Optimal
Primal bound 3
Dual bound 2
Gap 33.33% (tolerance: 36.67%)
Solution status feasible
3 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 782.32 (total)
0.28 (presolve)
0.00 (postsolve)
Nodes 170202
LP iterations 6791616 (total)
493334 (strong br.)
682673 (separation)
357533 (heuristics)
Optimal termination
Problem:
- Lower bound: 2.0
Upper bound: 3.0
Number of objectives: 1
Number of constraints: 0
Number of variables: 0
Sense: 1
Solver:
- Status: ok
Termination condition: optimal
Termination message: TerminationCondition.optimal
Solution:
- number of solutions: 0
number of solutions displayed: 0
--- placements ---
place 0 at 33
place 1 at 12
place 2 at 27
place 3 at 25
place 4 at 31
--- pairing penalties ---
assign (0, 1) to (unordered) (12, 33) for penalty 0 (dist: 3)
assign (0, 2) to (unordered) (27, 33) for penalty 0 (dist: 2)
assign (0, 3) to (unordered) (25, 33) for penalty 0 (dist: 2)
assign (0, 4) to (unordered) (31, 33) for penalty 1 (dist: 2)
assign (1, 2) to (unordered) (12, 27) for penalty 0 (dist: 3)
assign (1, 3) to (unordered) (12, 25) for penalty 0 (dist: 3)
assign (1, 4) to (unordered) (12, 31) for penalty 1 (dist: 5)
assign (2, 3) to (unordered) (25, 27) for penalty 0 (dist: 2)
assign (2, 4) to (unordered) (27, 31) for penalty 0 (dist: 4)
assign (3, 4) to (unordered) (25, 31) for penalty 1 (dist: 2)
Process finished with exit code 0
- Error: There were errors reading plugin package (please check also your plugin's metadata). Could not unzip file
- tensorflow: register numpy bfloat16 extension
- how to avoid row number in read_sql output
- How can I format Flask Werkzeug logging or move the request values to my logger format?
- pass both string and list to pandas .isin method
- ModuleNotFoundError: No module named 'taming'
- Baffled by Python math geometric output
- How to create a new text file using Python
- How do I get a list of tables, the schema, a dump, using the sqlite3 API?
- Leave incomplete line on screen when hitting Ctrl-C in ipython 5.0+
- How to understand the star quantifier (*) output
- Tk label not displaying on MacOS 13.0.1
- django-auditlog change list customization in django admin
- Issue installing matplotlib on Python 32-bit
- Discrete Real dimension spacing in scikit-optimize
- Annotate seaborn barplot with value from another column
- How to find the location of a specific "include" library in Python?
- Where are the python modules stored?
- Is there a way to destroy the tkinter messagebox without clicking OK?
- how to end flask app (e.g. control c) when i close/terminate a tkinter window
- pip install seems to be working, but when I try to run code I get thrown an error
- How can I format tickers for hh:mm in Python matplotlib over a fixed period
- How to run Uvicorn FastAPI server as a module from another Python file?
- Create random numbers with left skewed probability distribution
- How to set the logging level for the elasticsearch library differently to my own logging?
- Convert a python non-tail recursive function adding sub-folders info to a loop
- How to call a function that takes an argument in a Django template?
- Kivy app build successfully using buildozer but crashes when opened on android phone
- Django Template Not Displaying Vendor Contact and Description Fields
- Search for currency pairs