使用Julia JuMP求解TSP时添加子回路消除约束后程序卡顿无响应

核心问题 #

你遇到的并非代码语法错误，而是中等规模TSP整数规划模型在GLPK求解器下的计算效率瓶颈：

• 当n=58时，MTZ子回路消除约束会生成58×57=3306个额外约束，加上原有的流量平衡约束，整个模型的整数变量与约束规模远超GLPK的处理能力，导致求解时间极长（表现为看似“冻结”）。
• GLPK作为开源求解器，在处理中等以上规模的整数规划问题时，性能远不如商业求解器（如Gurobi、CPLEX），甚至不如开源的Cbc求解器。

可行优化方案 #

1. 更换更高效的求解器

将GLPK替换为性能更优的开源求解器Cbc，或商业求解器，示例代码：

# 替换为Cbc求解器
using Cbc
m = Model(Cbc.Optimizer)
# 后续变量、约束定义保持不变

若使用商业求解器（以Gurobi为例）：

using Gurobi
m = Model(Gurobi.Optimizer)

这类求解器针对整数规划问题有更高效的分支定界与剪枝策略，能大幅缩短求解时间。

2. 优化MTZ约束的变量类型

MTZ约束中的p变量无需强制声明为整数，最优解中p会自动取整数值。将整数变量改为连续变量，可降低求解器的计算负担：

# 移除Int约束，改为连续变量
@variable(m, p[i in N])
@constraint(m, [i in N], 1 <= p[i] <= n)
@constraint(m, p[1] == 1)
# 其余约束保持不变

3. 改用动态子回路消除（分支定界+割平面）

直接加入所有MTZ约束会引入大量冗余约束，更高效的方式是先求解无回路约束的模型，检测到子回路后再动态添加消除约束，直到无剩余子回路为止。示例思路：

function solve_tsp_with_subtour_elimination(D, n)
    N = 1:n
    m = Model(Cbc.Optimizer)
    @variable(m, x[i in N, j in N], Bin)
    @objective(m, Min, sum(x[i,j] * D[i,j] for i in N, j in N))
    @constraint(m, [i in N], sum(x[i,j] for j in N) == 1)
    @constraint(m, [j in N], sum(x[i,j] for i in N) == 1)
    @constraint(m, [i in N], x[i,i] == 0)
    
    while true
        optimize!(m)
        # 获取当前解的x变量值
        x_val = value.(x)
        # 检测子回路
        visited = falses(n)
        subtours = []
        for i in N
            if !visited[i]
                subtour = []
                current = i
                while !visited[current]
                    visited[current] = true
                    push!(subtour, current)
                    # 找到当前节点的下一个节点
                    current = findfirst(j -> x_val[current,j] ≈ 1, N)
                end
                if length(subtour) < n
                    push!(subtours, subtour)
                end
            end
        end
        if isempty(subtours)
            break
        end
        # 添加子回路消除约束
        for S in subtours
            @constraint(m, sum(x[i,j] for i in S, j in S) <= length(S) - 1)
        end
    end
    return objective_value(m), value.(x)
end

这种方式仅在需要时添加必要约束，避免了一开始就引入大量冗余约束，能显著提升求解效率。

补充说明 #

若不确定Julia是真冻结还是在计算，可开启求解器日志输出实时查看进度（以Cbc为例）：

m = Model(Cbc.Optimizer)
set_optimizer_attribute(m, "logLevel", 1) # 开启日志，跟踪求解过程

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