Python中PuLP无法返回整数线性规划模型解的求助

用户使用PuLP求解整数线性规划模型，确认问题存在整数解，排查变量和系数未发现异常，怀疑问题出在Qmax常量上。以下是完整代码：

from pulp import *
import networkx as nx


def solver(k,G,safe,Qmax):

    n = len(G)
    mu = [1.0/n]*n #Temporary

    unsafe = [i for i in range(n) if not i in safe]

    maet = LpProblem("ResolveMAET", LpMinimize)

    q = {i:(str(i)) for i in G}
    p = {(i,j):(str(i)+","+str(j)) for i in unsafe for j in G.neighbors(i)}

    q_var = LpVariable.dicts("q",[qi for qi in q.values()],lowBound = 0, upBound = Qmax, cat = 'Continuous')
    p_var = LpVariable.dicts("p",[pij for pij in p.values()],lowBound = 0, upBound= 1, cat = 'Binary')

    ### Objective function
    maet += lpSum([ mu[i] * q_var[q[i]] for i in q])

    ### Constraints

    #Markov Chain (Constraint n°1)
    for i in unsafe:
        Mij = 1.0 / G.degree(i) #Temporary

        maet += q_var[q[i]] - lpSum([Mij*q_var[q[j]] for j in G.neighbors(i)]) + Qmax * lpSum([p_var[p[i,j]] for j in G.neighbors(i)]) >= lpSum([ Mij*G[i][j]["weight"] for j in G.neighbors(i) ])

    #Sign (Constraint n°2)
    for i,j in p:

        maet += q_var[q[i]] - q_var[q[j]] - Qmax*p_var[p[i,j]] >= G[i][j]["weight"] - Qmax

    #At Most One Sign Per Node (Constraint n°3)
    for i in unsafe:
        maet += lpSum([p_var[p[i,j]] for j in G.neighbors(i)]) <= 1

    #At Most k Signs (Constraint n°4)
    maet += lpSum([p_var[p[i,j]] for i,j in p]) <= k

    #Shelter (Constraint n°5)
    for i in safe:
        maet += q_var[q[i]] <= 0 #=0 car q_var[q[i]] >= 0

    ### Résolution
    maet.solve()

    return [(i,j) for i,j in p if p_var[p[i,j]].varValue==1]


######################

v= 4.0 #km/h

# Graph initialization
G = nx.Graph()

file_density = open("densite.txt", 'r', encoding="UTF8")
mu = dict()

id_to_pos = {}
pos_to_id = {}

id = 0

for l in file_density:
    data = l.split(",")
    for i in range(0,len(data),3):
        vertex = (float(data[i][1:]),float(data[i+1]))
        number = int(data[i+2][6:])

        id_to_pos[id] = vertex
        pos_to_id[vertex] = id
        mu[id] = number

        id+=1

file_density.close()


# Ajout des noeuds
G.add_nodes_from([i for i in range(id)])


file_graph = open("graph.txt", 'r', encoding="UTF8")
file_graph.readline()

for l in file_graph:
    data = l.split(" ")
    weight = float(data[1])
    extremity1 = pos_to_id[(float(data[2].split(",")[0][1:]),float(data[2].split(",")[1][0:]))]
    extremity2 = pos_to_id[(float(data[3].split(",")[0][1:]),float(data[3].split(",")[1][0:]))]

    G.add_edge(extremity1, extremity2, weight=v*weight)

file_graph.close()


from random import randint

safe = [randint(0,len(G)) for i in range(10)]

print(solver(100,G,safe,60000.0))

补充说明：Qmax用于将条件约束转为线性约束，例如当$p_{i,j}=1$时，约束$q_i = t_{i,j} + q_j$会被转换为线性约束$q_i \geq t_{i,j} + q_j - Q_{max} (1-p_{i,j})$（对应代码中约束2）。

问题排查方向 #

1. Qmax取值不合理 #

• 取值过大：远超q变量实际可能的最大值，会导致约束松弛度过高，引发求解器数值稳定性问题，难以收敛到有效解。
• 取值过小：若Qmax不足以覆盖$q_i - q_j - G[i][j]["weight"]$的最大可能负值，当$p_{i,j}=0$时，约束$q_i - q_j - Qmax*p_{i,j} >= G[i][j]["weight"] - Qmax$会无法满足，直接导致模型不可行。

建议计算合理下界：$Qmax >= \max_{(i,j) \in E} (G[i][j]["weight"] + q_{max})$，其中$q_{max}$是任意不安全节点到安全节点的最长路径权重。

2. 约束2的线性化不完整 #

用户期望$p_{i,j}=1$时$q_i = t_{i,j} + q_j$，但当前仅实现了下界约束，缺少上界约束：

# 需添加的上界约束
maet += q_var[q[i]] - q_var[q[j]] + Qmax*p_var[p[i,j]] <= G[i][j]["weight"] + Qmax

仅保留下界会导致q变量取值无上限，模型可能无界或求解器无法锁定最优解。

3. safe节点生成错误 #

代码中randint(0, len(G))会生成len(G)这个不存在的节点（图节点范围是0~len(G)-1），导致unsafe列表计算错误，约束5引用无效变量，引发模型异常。修正为：

safe = [randint(0, len(G)-1) for i in range(10)]

4. 求解器配置问题 #

PuLP默认求解器可能对大规模整数规划支持不足，可添加日志查看求解状态：

maet.solve(PULP_CBC_CMD(msg=True))

通过日志判断模型是否不可行、无界或求解超时，也可尝试更换商业求解器（如Gurobi、CPLEX）提升求解能力。

5. 约束1逻辑偏差 #

约束1退化为马尔可夫链期望递推时（无$p_{i,j}=1$），表达式$q_i - \sum M_{ij}q_j >= \sum M_{ij}t_{ij}$可能过紧，导致模型不可行。需验证该约束的数学推导是否符合问题需求。

调试建议 #

• 构造小规模测试用例（手动创建小图、指定已知可行解），验证模型基础逻辑。
• 输出LP文件检查约束与变量：

  maet.writeLP("maet_model.lp")
  

  直观排查变量引用、系数、约束方向等错误。
• 逐步注释约束，定位导致模型不可行的具体约束。

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