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 approximate the location of the nodes given the constraints above-mentioned?

Here’s an incomplete minimal example script that:

• 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

Having absolutely no experience of OR-tools I’m unfortunately having a hard time implementing the constraints and correctly defining the objectives to be optimized.

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 = 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 smallest and largest geodesic distances 
max_distance = float('-inf')
min_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
            if distance < min_distance: min_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): (2, '≈'),
                       (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()

# Possible locations for a node (indices of grid points)
indices = [model.NewIntVarFromDomain(cp_model.Domain.FromValues(list(cell_indices)), f'p{i}') for i in range(num_nodes)]

# Ensure all nodes are placed at different grid locations
model.AddAllDifferent(indices)

# Objective variable to minimize (probably needed for distance approximation -> ≈)
total_deviation = model.NewIntVar(0, 1000, 'total_deviation')

# Accumulate the absolute differences between geodesic distances and target distances
deviations = []

# For each pair of nodes
for i1, i2 in combinations(range(num_nodes), 2):

    # Get the target distance that should ideally separate the pair
    target_dist, obj_type = objective_distances[(i1, i2)]
    
    # Set a variable for the actual distance
    actual_dist = model.NewIntVar(min_distance, max_distance, f'dist_{i1}_{i2}')
    
    # !! Don't know what to do from there !!
    
    # if obj_type == '≈':
    #
    #    - how to set a deviation from the target distance based on the set of all possible geodesic distances
    #    - how to minimize that deviation
    #
    # elif obj_type == '≥'
    #
    #    - how to make sure that 'actual_dist' is greater than or equal to the 'target_dist'