基于Pyomo的优化任务约束松弛以加速计算

问题背景 #

在始发站到目的地的wagon双目标调配任务中，核心目标是最大化调配数量同时最小化成本。但实际场景变量规模较大时，因要求各到达站wagon类型保持1/3比例的硬约束（原代码中CONSTRAINT 3），导致求解耗时极长甚至无法得到最优解；移除该约束后求解速度显著提升。需通过约束松弛（允许部分站点违反比例约束，同时在目标函数中加入惩罚项）或大M法加速求解。

实现方案 #

方案1：约束松弛+惩罚项（软约束） #

将原硬等式约束替换为允许偏差的软约束，引入正负偏差变量衡量违反程度，在成本目标中加入惩罚项，平衡约束遵守程度与求解速度。

核心步骤：

• 移除原有ratio_limit硬约束
• 定义非负偏差变量：pos_dev（Type1超出比例的偏差）、neg_dev（Type1不足比例的偏差）
• 构建松弛后的约束，允许实际值与目标比例存在偏差
• 在成本目标函数中加入偏差惩罚项，惩罚系数根据业务需求设定（系数越大，越倾向于遵守比例约束）

方案2：大M法（带松弛开关的约束） #

通过二进制变量控制是否允许站点违反比例约束，配合大M值限制偏差范围，同时在目标中加入违反约束的惩罚。适合需要严格控制违规站点的场景。

核心步骤：

• 引入二进制变量relax（1表示允许违规，0表示必须遵守约束）
• 用大M值限制偏差上限，relax=0时约束强制生效，relax=1时允许偏差
• 在目标中加入二进制变量的惩罚项，惩罚违规站点

修改后的完整代码示例 #

以下代码基于原模型，实现约束松弛+惩罚项方案，保留双目标求解逻辑：

import pandas as pd
from pyomo.environ import *

data = {'rem_type' : ['Type1', 'Type1', 'Type2', 'Type2', 'Type1', 'Type1', 'Type2', 'Type2'],
        'station_depart': ['A', 'A', 'A', 'A', 'F', 'F', 'F', 'F'],
        'station_arrive': ['B', 'C', 'B', 'C', 'B', 'C', 'B', 'C'],
        'cost': [100, 103, 111, 101, 105, 114, 95, 99]}
df_route = pd.DataFrame(data)

costs = df_route.set_index(['station_depart', 'station_arrive', 'rem_type']).to_dict()['cost']

# 始发站供应量数据
df_unload = pd.DataFrame({'rem_type' : ['Type1', 'Type2', 'Type1', 'Type2'], 
                          'station_depart' : ['A', 'A', 'F', 'F'],
                          'volume' : [5, 6, 4, 7]})

demands = df_unload.set_index(['station_depart', 'rem_type']).to_dict()['volume']

# 到达站容量数据
df_vrp = pd.DataFrame({'rem_type' : ['Type1', 'Type2', 'Type1', 'Type2'], 
                       'station_arrive' : ['B', 'B', 'C', 'C'],
                       'capacity' : [15, 15, 15, 15]})

capacity = df_vrp.set_index(['station_arrive', 'rem_type']).to_dict()['capacity']

# 提取集合信息
routes = {(s, t) for (s, t, _) in costs.keys()}
start_nodes = {k[0] for k in costs.keys()}
end_nodes = {k[1] for k in costs.keys()}
all_nodes = start_nodes | end_nodes

# 目标比例设定
ratios =    {('Type1', 'Type2') : (3, 1)}

# 构建模型
model = ConcreteModel("OP")

# 集合定义
model.N = Set(initialize=sorted(all_nodes))
model.S = Set(within=model.N, initialize=sorted(start_nodes))       # 始发站集合
model.T = Set(within=model.N, initialize=sorted(end_nodes))         # 到达站集合
model.R = Set(within=model.N * model.N, initialize=sorted(routes))  # 路径集合
model.V = Set(initialize=['Type1', 'Type2'])                        # wagon类型集合
model.ratio_pairs = Set(within=model.V * model.V, initialize=ratios.keys())

# 参数定义
model.cost         = Param(model.R, model.V, initialize=costs)
model.demand       = Param(model.S, model.V, initialize=demands, default=0)
model.capacity     = Param(model.T, model.V, initialize=capacity, default=0)
model.ratio_limits = Param(model.ratio_pairs, initialize=ratios, domain=Any)
model.penalty = Param(initialize=100, domain=PositiveReals)  # 偏差惩罚系数，按需调整

# 变量定义
model.send = Var(model.R, model.V, domain=NonNegativeIntegers)  # 各路径各类型wagon调配数量
model.pos_dev = Var(model.T, model.ratio_pairs, domain=NonNegativeIntegers)  # Type1超出比例的偏差
model.neg_dev = Var(model.T, model.ratio_pairs, domain=NonNegativeIntegers)  # Type1不足比例的偏差

# 目标1：最大化调配总量
model.obj1 = Objective(expr=sum_product(model.send), sense=maximize)

# 约束1：始发站供应量限制
@model.Constraint(model.S, model.V)
def demand_limit(model, s, v):
    return sum(model.send[r, v] for r in model.R if r[0] == s) <= model.demand[s, v]

# 约束2：到达站容量限制
@model.Constraint(model.T, model.V)
def capacity_limit(model, t, v):
    return sum(model.send[r, v] for r in model.R if r[1] == t) <= model.capacity[t, v]

# 松弛后的比例约束：允许存在偏差，用变量记录违反程度
@model.Constraint(model.T, model.ratio_pairs)
def relaxed_ratio_limit(model, t, v1, v2):
    type1_total = sum(model.send[r, v1] for r in model.R if r[1] == t)
    type2_total = sum(model.send[r, v2] for r in model.R if r[1] == t)
    ratio1, ratio2 = model.ratio_limits[v1, v2]
    return type1_total * ratio2 == type2_total * ratio1 + model.pos_dev[t, v1, v2] - model.neg_dev[t, v1, v2]

# 第一阶段求解：最大化调配数量
results = SolverFactory("glpk").solve(model)
results.write()

# 第二阶段求解：保证调配数量不低于第一阶段，最小化成本+偏差惩罚
model.min_send = Constraint(expr=sum_product(model.send) >= value(model.obj1))
model.obj2 = Objective(expr=sum_product(model.cost, model.send) + model.penalty * (sum(model.pos_dev) + sum(model.neg_dev)), sense=minimize)
model.obj1.deactivate()

results = SolverFactory("glpk").solve(model)
results.write()
model.send.display()
model.pos_dev.display()
model.neg_dev.display()

关键说明 #

• 惩罚系数调整：model.penalty需根据业务需求设定。若优先遵守比例约束，设置较大值；若优先求解速度和成本，设置较小值。
• 大M法适配：若需使用大M法，可添加二进制变量model.relax = Var(model.T, model.ratio_pairs, domain=Binary)，将松弛约束修改为：

  @model.Constraint(model.T, model.ratio_pairs)
  def bigm_ratio_limit(model, t, v1, v2):
      type1_total = sum(model.send[r, v1] for r in model.R if r[1] == t)
      type2_total = sum(model.send[r, v2] for r in model.R if r[1] == t)
      ratio1, ratio2 = model.ratio_limits[v1, v2]
      M = 1000  # 根据实际最大可能流量设定大M值
      return (type1_total * ratio2 - type2_total * ratio1) <= M * model.relax[t, v1, v2]
  

  同时在目标函数中加入model.penalty * sum(model.relax)作为违规惩罚。

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