在Pyomo中为混合整数规划问题设置复杂比例约束

我面临一个多目标混合整数规划问题，优化逻辑分为两阶段：

• 第一阶段：最大化分配至最优路线的wagon数量
• 第二阶段：在第一阶段确定的总量基础上，最小化运输成本

已知各出发站的wagon可用量、各到达站的容量约束，已实现基础约束与两阶段优化逻辑，但无法在Pyomo中添加Type2分配总量等于Type1分配总量的1/3的比例约束，特此求助。

由于决策变量x是整数类型，直接使用分数会引入浮点数精度问题，因此需将比例约束转化为整数等式：3 × Type2分配总量 = Type1分配总量。具体实现步骤如下：

1. 筛选Type1和Type2对应的路线索引
   从df_route中分别提取rem_type为Type1和Type2的行索引，用于计算两类wagon的总分配量。
2. 添加比例约束到模型
   基于筛选出的索引构建整数等式约束，确保Type2总量严格等于Type1总量的1/3。
3. 约束全局生效
   该约束属于全局约束，需在第一阶段最大化总量时就加入模型，确保后续第二阶段优化也满足比例要求。

import pandas as pd
import numpy as np
from pyomo.environ import *
from tqdm import tqdm  # 补充原代码缺失的导入

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)

# unload
df_unload = pd.DataFrame({'rem_type' : ['Type1', 'Type2', 'Type1', 'Type2'], 
                          'station_depart' : ['A', 'A', 'F', 'F'],
                          'volume' : [5, 6, 4, 7]})

# repair
df_vrp = pd.DataFrame({'rem_type' : ['Type1', 'Type2'], 
                       'station_arrive' : ['B', 'C'],
                       'capacity' : [15, 30]})

# create routes
routes_unload = dict()   
routes_vrp = dict()    

# create model
model = ConcreteModel("OP")

# decision vars
model.x = Var([idx for idx in df_route.index], domain=NonNegativeIntegers)

# --- 新增：筛选Type1和Type2的路线索引 ---
type1_indices = df_route[df_route['rem_type'] == 'Type1'].index.tolist()
type2_indices = df_route[df_route['rem_type'] == 'Type2'].index.tolist()

# --- 新增：添加比例约束 ---
model.type_ratio = Constraint(expr=3 * sum(model.x[idx] for idx in type2_indices) == sum(model.x[idx] for idx in type1_indices))

# obj function to maximize the number of allocated wags to routes
model.Size = Objective(expr=sum([model.x[i] for i in model.x]), sense=maximize)

# obj function to minimize cost of those wags from the first solution
model.Cost = Objective(
    expr=sum([df_route.loc[idx, "cost"] * model.x[idx] for idx in df_route.index]),
    sense=minimize,
)

model.Size.activate()
model.Cost.deactivate()

# routes with indexes for creating constraints
for idx in tqdm(df_route.index):
    vrp = (df_route.loc[idx, "rem_type"], df_route.loc[idx, "station_arrive"])
    if vrp not in routes_vrp:
        routes_vrp[vrp] = list()
    routes_vrp[vrp].append(idx)

    t_unload = (df_route.loc[idx, "rem_type"], df_route.loc[idx, "station_depart"])
    if t_unload not in routes_unload:
        routes_unload[t_unload] = list()
    routes_unload[t_unload].append(idx)

# constraints on the arrive/repair station
model.vrp = ConstraintList()
for v in df_vrp.index:
    t = (df_vrp.loc[v, "rem_type"], df_vrp.loc[v, "station_arrive"])
    if t in routes_vrp:
        model.vrp.add(
            sum([model.x[idx] for idx in routes_vrp[t]]) <= df_vrp.loc[v, "capacity"]
        )

# constraints on the depart station (how many wagons we have for allocation)
model.unload = ConstraintList()
for u in df_unload.index:
    t = (df_unload.loc[u, "rem_type"], df_unload.loc[u, "station_depart"])
    if t in routes_unload:
        model.unload.add(
            sum([model.x[idx] for idx in routes_unload[t]]) <= df_unload.loc[u, "volume"]
        )

results = SolverFactory("glpk").solve(model)

results.write()

# the result from the first objective function as the constraint to the second one
model.fix_count = Constraint(expr=sum([model.x[i] for i in model.x]) == model.Size())

model.Size.deactivate()
model.Cost.activate()

results = SolverFactory("glpk").solve(model)

results.write()

# create df with results of allocation to the routes
vals = []
for i in tqdm(model.x):
    vals.append(model.x[i].value)
df_results = df_route.dropna(
    subset=["rem_type", "station_depart", "station_arrive", "cost"],
    how="all",
)

df_results["total"] = vals

• 比例约束转化为整数等式是为了避免浮点数精度问题，确保整数规划求解器能正确处理约束。
• 若允许比例存在微小误差，可将等式约束改为不等式约束，例如3 * Type2总量 ≥ Type1总量 - ε和3 * Type2总量 ≤ Type1总量 + ε，其中ε为非负整数误差值。

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