基于OR-Tools的凹多边形内2D网格节点测地线距离约束布局优化

我要解决一项节点布局任务：在凹多边形内嵌的二维网格中布置一组节点，要求节点对间的测地线距离要么尽可能接近指定值（允许±1偏差），要么满足最小阈值要求。我用OR-Tools的CP-SAT模型编写了实现脚本，包含合规网格构建、测地线距离预计算、节点唯一性约束、约束兼容性判断及penalty最小化目标，但脚本始终无法返回结果，而我确认求解空间是有解的，恳请提供优化指导。

from ortools.sat.python import cp_model
from itertools import combinations
import networkx as nx

# 原始网格尺寸（宽、高）
w, h = 7, 5

# 合规网格单元格索引（凹多边形区域）
cell_indices = list(sorted(set(range(w * h)) - set([0, 1, 2, 3, 7, 8, 9, 10, 28, 29])))

# 构建合规网格拓扑图
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)

# 预计算所有单元格对的测地线距离（Dijkstra算法）
geodesic_distances = dict(nx.all_pairs_dijkstra_path_length(T))

# 获取最大测地线距离
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

# 待布置节点数量
num_nodes = 5

# 节点对的目标距离及约束类型：(目标距离, 约束类型≈/≥)
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, '≈')}

# 初始化CP-SAT模型
model = cp_model.CpModel()

# 定义节点位置布尔变量：node_at_position[node, index]表示节点node是否在位置index
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)
    # 约束每个位置最多放一个节点
    model.AddAtMostOne(at_most_one)

# 约束每个节点必须恰好放在一个位置
for node in range(num_nodes):
    model.AddExactlyOne(node_at_position[node, idx] for idx in cell_indices)

penalties = []

# 遍历所有节点对的约束要求
for (node1, node2), (target_distance, constraint_type) in objective_distances.items():
    # 遍历所有单元格对
    for i1, i2 in combinations(cell_indices, 2):
        distance = geodesic_distances[i1][i2]
        # 创建兼容性布尔变量
        is_compatible = model.NewBoolVar(f'compat_{node1}_{node2}_{i1}_{i2}')
        # 创建惩罚变量
        penalty = model.NewIntVar(0, max_distance, f'penalty_{node1}_{node2}_{i1}_{i2}')

        if constraint_type == '≈':
            # 距离在目标值±1范围内则兼容
            model.Add(is_compatible == (target_distance - 1 <= distance <= target_distance + 1))
        elif constraint_type == '≥':
            # 距离≥目标值则兼容
            model.Add(is_compatible == (distance >= target_distance))
        
        # 兼容时强制节点位置对应
        model.AddImplication(is_compatible, node_at_position[node1, i1])
        model.AddImplication(is_compatible, node_at_position[node2, i2])
        # 不兼容时添加惩罚
        model.Add(penalty == abs(distance - target_distance)).OnlyEnforceIf(is_compatible.Not())
        
        penalties.append(penalty)

# 目标：最小化总惩罚
model.Minimize(sum(penalties))

# 求解模型
solver = cp_model.CpSolver()
status = solver.Solve(model)

print("求解器状态:", solver.StatusName(status))    
if status == cp_model.FEASIBLE or status == cp_model.OPTIMAL:
    print("找到解决方案:")
    for node in range(num_nodes):
        for index in cell_indices:
            if solver.Value(node_at_position[node, index]):
                print(f'节点 {node} 位于位置 {index}')
else:
    print("未找到解决方案。")

优化指导 #

1. 减少变量爆炸问题 #

原代码为每对节点+每对单元格创建is_compatible和penalty变量，导致变量数量极度过剩（比如5个节点对应10对节点，若网格有25个单元格，会生成3000组冗余变量）。建议改用节点位置整数变量替代布尔变量：

# 定义节点位置整数变量，直接取值为单元格索引
pos = [model.NewIntVar(min(cell_indices), max(cell_indices), f'pos_{node}') for node in range(num_nodes)]
# 约束位置必须在合规单元格内
for node in range(num_nodes):
    model.AddAllowedAssignments([pos[node]], [[idx] for idx in cell_indices])
# 约束所有节点位置互不重复
model.AddAllDifferent(pos)

2. 优化距离约束与惩罚计算 #

利用AddMapDomain预查距离，结合条件约束计算惩罚，避免大量布尔变量：

total_penalty = model.NewIntVar(0, len(objective_distances)*max_distance, 'total_penalty')
penalty_terms = []

for (node1, node2), (target, type_) in objective_distances.items():
    # 定义节点对的测地线距离变量
    d = model.NewIntVar(0, max_distance, f'd_{node1}_{node2}')
    # 预构建位置对到距离的映射表
    distance_map = {}
    for i1 in cell_indices:
        for i2 in cell_indices:
            distance_map[(i1, i2)] = geodesic_distances[i1][i2]
    # 通过映射表关联位置变量与距离变量
    model.AddMapDomain([pos[node1], pos[node2]], d, distance_map)
    
    penalty = model.NewIntVar(0, max_distance, f'penalty_{node1}_{node2}')
    if type_ == '≈':
        # 距离在±1范围内惩罚为0，否则惩罚为距离与目标的差值绝对值
        is_close = model.NewBoolVar(f'close_{node1}_{node2}')
        model.Add(d >= target - 1).OnlyEnforceIf(is_close)
        model.Add(d <= target + 1).OnlyEnforceIf(is_close)
        model.Add(d < target - 1).OnlyEnforceIf(is_close.Not())
        model.Add(d > target + 1).OnlyEnforceIf(is_close.Not())
        
        model.Add(penalty == 0).OnlyEnforceIf(is_close)
        model.Add(penalty == cp_model.Abs(d - target)).OnlyEnforceIf(is_close.Not())
    elif type_ == '≥':
        # 满足最小距离要求惩罚为0，否则惩罚为目标与实际距离的差值
        meets_min = model.NewBoolVar(f'meets_{node1}_{node2}')
        model.Add(d >= target).OnlyEnforceIf(meets_min)
        model.Add(d < target).OnlyEnforceIf(meets_min.Not())
        
        model.Add(penalty == 0).OnlyEnforceIf(meets_min)
        model.Add(penalty == target - d).OnlyEnforceIf(meets_min.Not())
    
    penalty_terms.append(penalty)

model.Add(total_penalty == sum(penalty_terms))
model.Minimize(total_penalty)

3. 求解器参数调优 #

添加求解器参数缩短等待时间、开启调试日志：

solver = cp_model.CpSolver()
# 设置最大求解时间（秒）
solver.parameters.max_time_in_seconds = 60
# 开启搜索进度日志
solver.parameters.log_search_progress = True
# 使用组合搜索启发式
solver.parameters.search_branching = cp_model.PORTFOLIO_SEARCH
status = solver.Solve(model)

4. 提前过滤无效约束 #

预检查目标距离是否在网格的最小/最大测地线距离范围内，对于≥类型约束，若某对单元格的最大可能距离小于目标值，直接排除该位置组合，减少无效搜索空间。

内容的提问来源于stack exchange，提问作者solub