基于Gurobi Python的运输问题模型不可行原因排查与修正求助

模型不可行的核心原因 #

1. 未按花卉品类拆分变量与约束 #

原问题涉及玫瑰（R）、郁金香（T）、百合（L）3种花卉，但代码仅定义了单一运输变量x，未区分不同品类的运输流量。这导致供需约束无法对应到具体花卉，本质上是模型逻辑错误。

2. 数据读取与索引错误 #

• 节点数据读取不完整：循环for line in f.readlines()[1:-4]跳过了最后4行数据，导致仓库w7和所有零售商r8/r9/r10的信息未被加载，后续约束引用这些节点时会出现数据缺失。
• 供需约束索引错误：
  • 产地约束仅提取了百合（capL，索引5）的供应量，完全忽略了玫瑰（capR，索引3）和郁金香（capT，索引4）的供应。
  • 零售商约束仅提取了百合（demL，索引8）的需求量，未考虑玫瑰（demR，索引6）和郁金香（demT，索引7）的需求。
• 运输成本参数错误：代码尝试读取文件最后一行获取costFactor，但原数据文件无此内容；同时未区分三种花卉的单位公里运输成本（2/3/4），导致成本计算完全不符合要求。

3. 流量平衡约束逻辑不匹配 #

由于未拆分品类，节点流量平衡将不同花卉的流量混为一谈，即使总供需数值匹配，也无法满足单一品类的供需路径要求，直接导致模型不可行。

修正方案 #

步骤1：修复数据读取逻辑 #

改用split()分割数据，加载所有节点信息，并正确读取三种花卉的供需数据：

import gurobipy as gp
import numpy as np

# 1) 读取数据
instname = 'nodes'
mainFolderInput = r'C:\Users\karsp\Documents\Python Scripts\'
fileAddress = mainFolderInput + instname + '.txt'

f = open(fileAddress)
fDict = {}
rDict = {}
wDict = {}

# 读取表头后所有行，去掉[-4]的错误截断
for line in f.readlines()[1:]:
    asList = line.strip().split()
    lID = asList[0]
    # 转换数值列为整数
    asList[1:] = [int(num) for num in asList[1:]]
    if lID.startswith("f"): 
        fDict[lID] = asList
    elif lID.startswith("r"): 
        rDict[lID] = asList
    elif lID.startswith("w"): 
       wDict[lID] = asList

vDict = {}
vDict.update(fDict)
vDict.update(rDict)
vDict.update(wDict)

# 定义各花卉的单位公里成本
flower_costs = {
    'R': 2,
    'T': 3,
    'L': 4
}

步骤2：按品类定义运输变量 #

为每种花卉单独定义运输变量，变量维度为(弧, 花卉品类)：

# 初始化模型
TPModel = gp.Model('1CK100-TP')
TPModel.modelSense = gp.GRB.MINIMIZE

# 定义弧集与花卉品类
arcs = [('f1','w3'),('f1','w4'),('f2', 'w4'),('f2','w5'),('w3','w6'),('w4','w6'),('w4','r9'),('w4','w7'),('w5','w7'),('w5','r9'),('w6','r8'),('w6','r9'),('w6','r10'),('w7','r8'),('w7','r9'),('w7','r10')]
flowers = ['R', 'T', 'L']

# 按品类创建整数运输变量
x = TPModel.addVars(arcs, flowers, vtype=gp.GRB.INTEGER, name='x')

步骤3：修正运输成本计算 #

针对每种花卉，结合运输距离和对应单位成本计算总成本：

# 计算两点间距离（取整）
def calculate_distance(node1, node2):
    return np.math.trunc(np.hypot(node1[1] - node2[1], node1[2] - node2[2]))

# 构建分品类的运输成本字典
transport_costs = {}
for flower in flowers:
    cost_per_km = flower_costs[flower]
    # 产地到仓库
    for f_id in fDict:
        for w_id in wDict:
            if (f_id, w_id) in arcs:
                dist = calculate_distance(fDict[f_id], wDict[w_id])
                transport_costs[(f_id, w_id, flower)] = cost_per_km * dist
    # 仓库到零售商
    for w_id in wDict:
        for r_id in rDict:
            if (w_id, r_id) in arcs:
                dist = calculate_distance(wDict[w_id], rDict[r_id])
                transport_costs[(w_id, r_id, flower)] = cost_per_km * dist
    # 仓库到仓库
    for w1_id in wDict:
        for w2_id in wDict:
            if (w1_id, w2_id) in arcs and w1_id != w2_id:
                dist = calculate_distance(wDict[w1_id], wDict[w2_id])
                transport_costs[(w1_id, w2_id, flower)] = cost_per_km * dist

步骤4：修正目标函数 #

按品类和弧汇总运输成本：

TPModel.setObjective(
    gp.quicksum(transport_costs[i, j, f] * x[i, j, f] for (i,j) in arcs for f in flowers),
    sense=gp.GRB.MINIMIZE
)

步骤5：修正流量平衡约束 #

为每个节点、每个品类单独设置流量平衡约束：

v = vDict.keys()

# 产地约束：每个品类的流出量等于供应量
for f_id in fDict:
    # 玫瑰（R）
    TPModel.addConstr(
        gp.quicksum(x[f_id, j, 'R'] for j in v if (f_id, j) in arcs) == fDict[f_id][3],
        name=f'balance_f_{f_id}_R'
    )
    # 郁金香（T）
    TPModel.addConstr(
        gp.quicksum(x[f_id, j, 'T'] for j in v if (f_id, j) in arcs) == fDict[f_id][4],
        name=f'balance_f_{f_id}_T'
    )
    # 百合（L）
    TPModel.addConstr(
        gp.quicksum(x[f_id, j, 'L'] for j in v if (f_id, j) in arcs) == fDict[f_id][5],
        name=f'balance_f_{f_id}_L'
    )

# 零售商约束：每个品类的流入量等于需求量
for r_id in rDict:
    # 玫瑰（R）
    TPModel.addConstr(
        gp.quicksum(x[i, r_id, 'R'] for i in v if (i, r_id) in arcs) == rDict[r_id][6],
        name=f'balance_r_{r_id}_R'
    )
    # 郁金香（T）
    TPModel.addConstr(
        gp.quicksum(x[i, r_id, 'T'] for i in v if (i, r_id) in arcs) == rDict[r_id][7],
        name=f'balance_r_{r_id}_T'
    )
    # 百合（L）
    TPModel.addConstr(
        gp.quicksum(x[i, r_id, 'L'] for i in v if (i, r_id) in arcs) == rDict[r_id][8],
        name=f'balance_r_{r_id}_L'
    )

# 仓库约束：每个品类的流入量等于流出量
for w_id in wDict:
    for flower in flowers:
        TPModel.addConstr(
            gp.quicksum(x[w_id, j, flower] for j in v if (w_id, j) in arcs) 
            == gp.quicksum(x[i, w_id, flower] for i in v if (i, w_id) in arcs),
            name=f'balance_w_{w_id}_{flower}'
        )

# 非负约束（整数变量默认非负，可省略）
TPModel.addConstrs((x[i,j,f] >= 0 for (i,j) in arcs for f in flowers), name='nonnegativity')

步骤6：求解并输出结果 #

TPModel.optimize()

# 输出最优解或排查不可行原因
if TPModel.status == gp.GRB.OPTIMAL:
    print("最优运输方案：")
    for (i,j,f) in x.keys():
        if x[i,j,f].X > 0:
            print(f"从{i}运输{x[i,j,f].X}箱{f}到{j}")
    print(f"总运输成本：{TPModel.ObjVal}")
elif TPModel.status == gp.GRB.INFEASIBLE:
    print("模型仍不可行，生成不可行约束集文件：")
    TPModel.computeIIS()
    TPModel.write("model.ilp")

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