Question: 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 AddImplication to 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.") 

Answer: 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