基于Gekko的设施选址问题建模优化及求解方案咨询

我用Python的Gekko构建类似设施选址的优化模型，目标是最小化选中设施与对应县的距离总和，约束条件如下：

• 每个县仅能映射到一个设施；
• 若某设施被选中，其服务的总人口需达到所有县总人口的至少10%；
• 不限制选中设施的数量，由算法自行决定最优选择。

目前求解时持续出现@error: Solution Not Found错误，想咨询：

1. 是否有更优的建模方式？
2. 若现有建模无法解决，推荐哪些其他算法/工具？

相关数据 #

注：以下为真实数据集的极小样本，仅作格式展示

import numpy as np

county_idx = [0,1,2]
selected_facilities = [0, 1, 2]
distances_mapped = np.array([[193, 85, 226],
                    [139, 112, 241],
                    [175, 110, 249]])

populations_mapped = [981447, 327286, 176622]

# 错误的截断值：实际应为总人口的10%，总人口=981447+327286+176622=1485355，10%为148535.5
cutoff_population = 153

现有代码 #

from gekko import GEKKO
import numpy as np

m = GEKKO(remote=False)
n_rows = len(county_idx)
n_cols = len(selected_facilities)

# 设施选择二进制变量
facility_assignment = m.Array(m.Var, (n_cols), lb=0, ub=1, integer=True)

# 县-设施映射二进制变量
assignment_matrix = m.Array(m.Var, (n_rows, n_cols), lb=0, ub=1, integer=True)

# 目标函数：最小化距离总和
m.Minimize(m.sum(assignment_matrix * distances_mapped))

# 约束：每个县仅分配给一个设施
for i in county_idx:
    m.Equation(m.sum([assignment_matrix[(i, j)] for j in selected_facilities]) == 1)

# 绑定设施选择与映射关系（现有错误约束）
for j in selected_facilities:
    m.Equation(facility_assignment[j] == m.sum([assignment_matrix[(i, j)] for i in county_idx]))

# 约束：选中的设施服务人口需达标（现有约束逻辑有瑕疵）
for j in selected_facilities:
    m.Equation(m.sum([assignment_matrix[(i, j)] * populations_mapped[i] for i in county_idx]) >= cutoff_population * facility_assignment[j])

# 使用APOPT整数规划求解器
m.options.SOLVER = 1
m.options.IMODE = 3

# 求解
m.solve(disp=True, debug=True)

优化建模方式 #

导致无解的核心问题在于设施选择与映射关系的约束错误，以及截断人口值的计算错误，优化方案如下：

1. 修正设施选择与映射的绑定约束 #

现有代码中用facility_assignment[j] == sum(assignment_matrix[i,j])是错误的：当多个县分配给同一个设施时，sum(assignment_matrix[i,j])会大于1，但facility_assignment[j]是0/1变量，等式无法满足，直接导致无解。

正确的绑定逻辑应该用大M法：

• 若有任何县分配给设施j，则facility_assignment[j]必须为1；
• 若设施j未被选中，则所有县都不能分配给它。

替换原有绑定约束为：

# 大M值取县的总数（最大可能的分配数量）
M = n_rows

for j in selected_facilities:
    # 约束1：未选中的设施不能分配任何县
    for i in county_idx:
        m.Equation(assignment_matrix[i,j] <= facility_assignment[j])
    # 约束2：有县分配的设施必须被标记为选中
    m.Equation(m.sum(assignment_matrix[i,j] for i in county_idx) <= M * facility_assignment[j])

2. 修正截断人口值的计算 #

你需要的是选中设施服务人口≥总人口的10%，所以正确计算cutoff_population：

total_pop = sum(populations_mapped)
cutoff_population = 0.1 * total_pop  # 示例中为148535.5

3. 调整人口约束的严谨性 #

原约束逻辑正确，但可以明确处理整数/浮点数的问题，同时确保当设施未被选中时，约束自动满足（左边为0，右边为0）：

for j in selected_facilities:
    assigned_pop = m.sum(assignment_matrix[i,j] * populations_mapped[i] for i in county_idx)
    m.Equation(assigned_pop >= cutoff_population * facility_assignment[j])

优化后的完整代码 #

from gekko import GEKKO
import numpy as np

county_idx = [0,1,2]
selected_facilities = [0, 1, 2]
distances_mapped = np.array([[193, 85, 226],
                    [139, 112, 241],
                    [175, 110, 249]])

populations_mapped = [981447, 327286, 176622]
total_pop = sum(populations_mapped)
cutoff_population = 0.1 * total_pop  # 正确的10%截断值

m = GEKKO(remote=False)
n_rows = len(county_idx)
n_cols = len(selected_facilities)

# 设施选择二进制变量
facility_assignment = m.Array(m.Var, (n_cols), lb=0, ub=1, integer=True)

# 县-设施映射二进制变量
assignment_matrix = m.Array(m.Var, (n_rows, n_cols), lb=0, ub=1, integer=True)

# 目标函数：最小化距离总和
m.Minimize(m.sum(assignment_matrix * distances_mapped))

# 约束：每个县仅分配给一个设施
for i in county_idx:
    m.Equation(m.sum(assignment_matrix[i,j] for j in selected_facilities) == 1)

# 大M法绑定设施选择与映射关系
M = n_rows
for j in selected_facilities:
    # 未选中的设施不能分配任何县
    for i in county_idx:
        m.Equation(assignment_matrix[i,j] <= facility_assignment[j])
    # 有县分配的设施必须被标记为选中
    m.Equation(m.sum(assignment_matrix[i,j] for i in county_idx) <= M * facility_assignment[j])

# 约束：选中的设施服务人口需≥总人口的10%
for j in selected_facilities:
    assigned_pop = m.sum(assignment_matrix[i,j] * populations_mapped[i] for i in county_idx)
    m.Equation(assigned_pop >= cutoff_population * facility_assignment[j])

m.options.SOLVER = 1  # APOPT求解器
m.options.IMODE = 3
m.solve(disp=True, debug=True)

# 输出结果
print("选中的设施：", [j for j in selected_facilities if facility_assignment[j].value[0] == 1])
print("县-设施映射：")
for i in county_idx:
    for j in selected_facilities:
        if assignment_matrix[i,j].value[0] == 1:
            print(f"县{i} → 设施{j}")

替代算法/工具推荐 #

如果Gekko的APOPT求解器仍无法解决（比如大规模数据集），可以选择以下专门的整数规划/组合优化工具：

1. PuLP + CBC/Gurobi #

• PuLP是Python的线性规划建模库，语法简洁，支持多种求解器；
• CBC是开源的整数规划求解器，适合中小规模问题；Gurobi是商业求解器，性能更强（需授权）。

2. OR-Tools #

• Google开源的组合优化工具包，内置MIP求解器、启发式算法，适合设施选址这类组合优化问题，性能优异，文档完善。

3. SCIP #

• 开源的高性能整数规划求解器，支持复杂的约束和大规模问题，可通过PySCIPOpt库在Python中调用。

这些工具都针对整数规划问题做了深度优化，相比Gekko更适合纯组合优化场景，尤其是大规模数据集下的求解效率更高。

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