Оптимизация размещения узлов в 2D-сетке для соответствия определенным геодезическим расстояниям

Я бы не делал этого таким образом.

Создайте все упорядоченные пары точек. (i, j) и (j, i). Для каждой пары вычислите геодезическое расстояние.

Теперь посмотрите на ваши требования и ограничения: for (a, b) == k: создайте 1 bool var для каждой совместимой пары точек для (a, b) >= k создайте 1 bool var для каждой совместимой пары точек и добавьте термин bool_var * штраф этого выбора.

Теперь вам просто нужно добавить все ограничения at_most_one для обеспечения согласованности. Получить их все немного сложно, но что-то вроде:

для каждого узла и каждой позиции: создайте 1 bool_var. Добавьте ограничение «точно_один» для всех точек для 1 узла, добавьте 1 at_most_one для всех узлов для 1 точки. Затем добавьте импликации (пара a, b использует узел i, j), подразумевающую (a, i) == true и (b, j) == true.

Наконец, минимизируйте сумму штрафов.

Если вы не обязаны ortools, вот решение с использованием pyomo и HiGHS решателя. Как минимум, это может дать вам некоторое представление о правильном синтаксисе ortools для реализации некоторых из этих типов ограничений, если эта формулировка работает. Я не знаком с синтаксисом ortools, но это соответствует контексту того, что @Laurent Perron упомянул выше.

Это довольно толстая BIP (программа двоичных целых чисел), и я ввел в решение абсолютный пробел, так что, если он окажется в пределах 1,0 от теоретического предела, он остановится, а не продолжит решение. При этом на моей машине примерно за 10 минут выдается оптимальный ответ с решением, показанным внизу, с общим штрафом в 3 единицы расстояния для показанного примера.

Код:

"""
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')

Выход:

...
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