如何在Pyomo+GLPK的目标函数中引入绝对值求解节点链路成本优化

问题背景 #

需要求解节点间链路建设的最优成本，目标为最小化总成本。当前存在以下问题：

• 初始模型可求解变量，但成本最优值不符合预期，原因是需对每项成本取绝对值
• 直接使用abs()会使问题变为非线性，无法通过GLPK求解
• 尝试引入新变量p处理绝对值约束时，出现参数d的索引错误，需修正代码实现绝对值的正确线性化处理

成本项计算逻辑：cost+=(model.a[i,k]-model.a[j,k])*model.c*model.d[i,j]，其中：

• 当model.a[j,k]=1时，该值为负
• 当model.a[i,k]=1时，该值为正

初始实现代码 #

from pyomo.environ import *
 
# 创建Concrete模型
model = ConcreteModel()

# 集合定义
model.i = Set(initialize=[1,2,3,4,5,6,7,8,9,10,11,12,13], doc='节点集合')
model.k = Set(initialize=['Orange','SFR', 'Bouygues'], doc='运营商集合')

# 参数定义
model.c = Param(initialize=25, doc='将现有链路转换为骨干链路的成本（欧元/公里）')

links={
        (1, 2) : 1.8,
        (1, 7) : 1,
        (1, 13) : 5.4,
        (2, 8) : 2.3,
        (2, 3) : 1.7,
        (2, 5) : 7,
        (2, 7) : 2,
        (2, 12) : 3,
        (3, 4) : 2,
        (3, 10) : 6.5,
        (4, 5) : 1,
        (4, 6) : 2,
        (5, 8) : 5,
        (5, 10) : 1,
        (5, 11) : 1.5,
        (6, 11) : 2.1,
        (7, 12) : 2,
        (8, 9) : 2,
        (8, 13) : 0.7,
        (9, 10) : 1.1,
        (10, 11) : 1,
        (12, 13) : 2.5,
        }
model.d=Param(model.i, model.i,default=0, initialize=links, doc='节点间距离（单位：10公里）')

# 变量定义
model.a = Var(model.i, model.k, within=Binary, doc='二进制变量，节点i是否属于运营商k（0=不属于，1=属于）')

# 约束条件
def allocation_rule(model, i):
  return sum(model.a[i,k] for k in model.k) == 1
model.allocation = Constraint(model.i, rule=allocation_rule, doc='每个节点只能属于一个运营商')

def minimum_rule(model, k):
  return sum(model.a[i,k] for i in model.i) >= 2
model.minimum = Constraint(model.k, rule=minimum_rule, doc='每个运营商至少拥有2个节点')

# 目标函数
def totalcost(model):
    cost=0
    for i in model.i:
        for j in model.i:
            if model.d[i,j]!=0:
                for k in model.k:
                      cost+=(model.a[i,k]-model.a[j,k])*model.c*model.d[i,j]
    return cost
model.z = Objective(rule=totalcost, sense=minimize, doc='最小化三个子网骨干链路的建设成本')        

## 求解与输出
optimizer = SolverFactory("glpk",executable='/usr/bin/glpsol')     # 创建GLPK求解器实例
optimizer.solve(model)                # 求解模型
model.display()

尝试修改后的错误代码及问题 #

修改代码引入变量p处理绝对值，但出现错误：ERROR:pyomo.core:Rule failed for Param 'd' with index (1, 2)

错误代码：

from pyomo.environ import *

# 创建Concrete模型
model = ConcreteModel()

# 集合定义（错误：重复定义覆盖了节点集合）
model.i = Set(initialize=[1,2,3,4,5,6,7,8,9,10,11,12,13], doc='节点集合')
model.i = Set(initialize=['Orange','SFR', 'Bouygues'], doc='运营商集合')

# 参数定义
model.c = Param(initialize=25, doc='将现有链路转换为骨干链路的成本（欧元/公里）')

links={
    (1, 2) : 1.8,
    (1, 7) : 1,
    (2, 3) : 1.7,
    (2, 5) : 7,
    (2, 7) : 2,
    (2, 12) : 3,
    (3, 4) : 2,
    (3, 10) : 6.5,
    (4, 5) : 1,
    (4, 6) : 2,
    (5, 8) : 5,
    (5, 10) : 1,
    (5, 11) : 1.5,
    (6, 11) : 2.1,
    (7, 12) : 2,
    (8, 9) : 2,
    (8, 13) : 0.7,
    (9, 10) : 1.1,
    (10, 11) : 1,
    (12, 13) : 2.5,
    (1, 13) : 5.4,
    (2, 8) : 2.3,
    }  

model.d=Param(model.i, model.i,default=0, initialize=links, doc='节点间距离（单位：10公里）')

# 变量定义（错误：p的维度不匹配）
model.a = Var(model.i, model.k, within=Binary, doc='二进制变量，节点i是否属于运营商k（0=不属于，1=属于）')
model.p = Var(model.i,model.k, within=(0.0,None), doc='骨干链路建设成本p_ij')

# 约束条件（错误：约束规则未正确接收索引参数，遍历逻辑错误）
def allocation_rule(model, i):
  return sum(model.a[i,k] for k in model.k) == 1
model.allocation = Constraint(model.i, rule=allocation_rule, doc='每个节点只能属于一个运营商')

def minimum_rule(model, k):
  return sum(model.a[i,k] for i in model.i) >= 2
model.minimum = Constraint(model.k, rule=minimum_rule, doc='每个运营商至少拥有2个节点')

def absolute_rule1(model):
  return model.p >=(model.a[i,k]-model.a[j,k])*model.c*model.d[i,j]
model.absolute1 = Constraint(model.i, rule=absolute_rule1, doc='取正成本约束')

def absolute_rule2(model):
  for i in model.i:
    for j in model.i:
      if model.d[i,j]!=0:
        for k in model.k:
          return model.p >=-(model.a[i,k]-model.a[j,k])*model.c*model.d[i,j]
model.absolute2 = Constraint(model.i, rule=absolute_rule2, doc='取正成本约束')

# 目标函数（错误：累加逻辑错误，提前return）
def totalcost(model):
        cost=0
        for i in model.i:
          for j in model.i:
            if model.d[i,j]!=0:
                for k in model.k:
                    cost+=model.p
                    return cost
model.z = Objective(rule=totalcost, sense=minimize, doc='最小化三个子网骨干链路的建设成本')        

修正后的代码及说明 #

错误原因分析 #

1. 集合定义错误：重复赋值model.i，将节点集合覆盖为运营商集合，导致model.d的节点索引无法匹配
2. 变量维度错误：p的维度应为(i,j,k)，对应每个链路和运营商的绝对值成本
3. 约束规则错误：未正确遍历所有(i,j,k)组合，约束规则未接收索引参数
4. 目标函数错误：累加逻辑错误，提前返回导致未计算所有项

修正后完整代码 #

from pyomo.environ import *

# 创建Concrete模型
model = ConcreteModel()

# 集合定义：恢复正确的节点和运营商集合
model.i = Set(initialize=[1,2,3,4,5,6,7,8,9,10,11,12,13], doc='节点集合')
model.k = Set(initialize=['Orange','SFR', 'Bouygues'], doc='运营商集合')

# 参数定义
model.c = Param(initialize=25, doc='将现有链路转换为骨干链路的成本（欧元/公里）')

links={
        (1, 2) : 1.8,
        (1, 7) : 1,
        (1, 13) : 5.4,
        (2, 8) : 2.3,
        (2, 3) : 1.7,
        (2, 5) : 7,
        (2, 7) : 2,
        (2, 12) : 3,
        (3, 4) : 2,
        (3, 10) : 6.5,
        (4, 5) : 1,
        (4, 6) : 2,
        (5, 8) : 5,
        (5, 10) : 1,
        (5, 11) : 1.5,
        (6, 11) : 2.1,
        (7, 12) : 2,
        (8, 9) : 2,
        (8, 13) : 0.7,
        (9, 10) : 1.1,
        (10, 11) : 1,
        (12, 13) : 2.5,
        }
model.d=Param(model.i, model.i, default=0, initialize=links, doc='节点间距离（单位：10公里）')

# 变量定义：修正p的维度为(i,j,k)，非负实数
model.a = Var(model.i, model.k, within=Binary, doc='二进制变量，节点i是否属于运营商k（0=不属于，1=属于）')
model.p = Var(model.i, model.i, model.k, within=NonNegativeReals, doc='链路(i,j)对应运营商k的绝对值成本')

# 约束条件
def allocation_rule(model, i):
    return sum(model.a[i,k] for k in model.k) == 1
model.allocation = Constraint(model.i, rule=allocation_rule, doc='每个节点只能属于一个运营商')

def minimum_rule(model, k):
    return sum(model.a[i,k] for i in model.i) >= 2
model.minimum = Constraint(model.k, rule=minimum_rule, doc='每个运营商至少拥有2个节点')

# 绝对值线性化约束：对每个存在的链路(i,j)和运营商k
def absolute_upper_rule(model, i, j, k):
    if model.d[i,j] == 0:
        return Constraint.Skip  # 跳过无链路的节点对
    return model.p[i,j,k] >= (model.a[i,k] - model.a[j,k]) * model.c * model.d[i,j]

def absolute_lower_rule(model, i, j, k):
    if model.d[i,j] == 0:
        return Constraint.Skip  # 跳过无链路的节点对
    return model.p[i,j,k] >= -(model.a[i,k] - model.a[j,k]) * model.c * model.d[i,j]

model.absolute_upper = Constraint(model.i, model.i, model.k, rule=absolute_upper_rule, doc='绝对值上界约束')
model.absolute_lower = Constraint(model.i, model.i, model.k, rule=absolute_lower_rule, doc='绝对值下界约束')

# 目标函数：累加所有链路和运营商的绝对值成本
def totalcost(model):
    total = 0
    for i in model.i:
        for j in model.i:
            if model.d[i,j] != 0:
                for k in model.k:
                    total += model.p[i,j,k]
    return total

model.z = Objective(rule=totalcost, sense=minimize, doc='最小化三个子网骨干链路的建设成本')

# 求解与输出
optimizer = SolverFactory("glpk", executable='/usr/bin/glpsol')
optimizer.solve(model)
model.display()

修正要点说明 #

1. 集合修复：保留节点集合model.i和运营商集合model.k的正确定义，避免覆盖
2. 变量维度调整：model.p定义为(i,j,k)，对应每个链路和运营商的绝对值成本，确保维度匹配
3. 绝对值线性化：通过两个线性约束p[i,j,k] >= X和p[i,j,k] >= -X实现绝对值|X|的线性表达，其中X=(a[i,k]-a[j,k])*c*d[i,j]
4. 约束优化：跳过无链路的节点对，减少不必要的约束
5. 目标函数修正：正确累加所有p[i,j,k]变量，确保计算总成本

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