约束编程：OR-Tools中用变量作矩阵索引的优化问题求解

你遇到的是标准二次分配问题（QAP），两个报错/卡壳的原因非常明确：

1. 第一种写法错误：建模阶段求解器还没有对变量赋值，你用np.where查询变量的取值，本质是在操作未求解的变量对象，根本拿不到0/1的实际结果，必然触发索引错误。
2. 第二种写法思路完全正确，你担心的双布尔变量乘积“非线性”问题，在OR-Tools的CP-SAT求解器里不需要额外处理——两个布尔变量的乘积仍然是布尔值，求解器原生支持这类约束，只需要用辅助变量显式绑定乘积关系即可，不需要换求解库。

首先统一变量定义避免下标混乱：

• 点位索引i,j对应原点位编号'1'、'2'、'3'
• 名称索引k,l对应原名称标识'a'、'b'、'c'
• 布尔变量x[i,k]取1表示点位i分配名称k，取0表示不分配
• dist[i,j]为点位i到j的距离，flow[k,l]为名称k到l的通行次数
• 优化目标和你推导的完全一致：最小化所有通行需求的总距离，即$\sum_{k,l}\sum_{i,j} flow[k,l] * dist[i,j] * x[i,k] *x[j,l]$

from ortools.sat.python import cp_model
import numpy as np

# 输入矩阵定义
dist = np.array([
    [0, 10, 20],
    [10, 0, 30],
    [20, 30, 0]
])
flow = np.array([
    [0, 1, 3],
    [1, 0, 2],
    [3, 2, 0]
])
n = len(dist)

model = cp_model.CpModel()

# 定义分配布尔变量
x = {}
for i in range(n):
    for k in range(n):
        x[i, k] = model.NewBoolVar(f'x_{i}_{k}')

# 添加置换约束：每个点位唯一名称，每个名称唯一对应点位
for i in range(n):
    model.AddExactlyOne(x[i, k] for k in range(n))
for k in range(n):
    model.AddExactlyOne(x[i, k] for i in range(n))

# 构造目标函数
cost_terms = []
for k in range(n):
    for l in range(n):
        current_flow = flow[k, l]
        if current_flow == 0:
            continue  # 跳过0流量项减少计算量
        for i in range(n):
            for j in range(n):
                current_dist = dist[i, j]
                if current_dist == 0:
                    continue  # 跳过0距离项（同点位）
                # 定义辅助变量存储两个布尔变量的乘积
                prod = model.NewBoolVar('')
                # 绑定乘积约束：prod=1 当且仅当x[i,k]和x[j,l]同时为1
                model.Add(prod == x[i, k] * x[j, l])
                cost_terms.append(current_flow * current_dist * prod)

model.Minimize(sum(cost_terms))

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

# 结果输出
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
    print(f"最小总通行距离：{solver.ObjectiveValue()}")
    point_id = {0: '1', 1: '2', 2: '3'}
    name_id = {0: 'a', 1: 'b', 2: 'c'}
    for i in range(n):
        for k in range(n):
            if solver.Value(x[i, k]) == 1:
                print(f"点位{point_id[i]} -> 名称{name_id[k]}")
else:
    print("无可行解")

针对你给出的示例矩阵，求解得到最优值为200，对应分配方案为：

• 点位'1' -> 名称'c'
• 点位'2' -> 名称'a'
• 点位'3' -> 名称'b'

• CP-SAT处理双布尔变量乘积的本质是自动添加了三个线性约束：prod ≤ x[i,k]、prod ≤ x[j,l]、prod ≥ x[i,k] + x[j,l] -1，完全不需要你手动处理非线性问题。
• 当点位/名称数量小于20时，上述代码的求解效率足够；如果规模更大，可以考虑针对QAP设计的专用启发式算法，常规场景用OR-Tools完全可以满足需求。

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