使用epsilon-constraint方法在Gurobi中求解多目标问题获取Pareto前沿

前置说明 #

你的原始模型存在两个可优化点：

1. 已定义允许的电源-负荷配对pi，但未在模型中约束，会出现非法调度结果，代码中已修正为仅在pi配对范围内定义发电量变量
2. 定义的二进制变量z未在目标和约束中使用，暂时注释，后续需要可自行启用

ε-约束法实现逻辑 #

• 首先单独优化每个目标，得到每个目标的取值上下界，确定ε阈值的可调区间
• 选择成本作为主优化目标，将CO₂排放、清洁能源发电量两个目标转化为带ε约束的限制条件
• 按固定步长调整ε阈值迭代求解，过滤得到所有非支配的Pareto最优解

完整可运行代码 #

from gurobipy import Model, GRB, quicksum
import numpy as np

# 原始参数定义
power = ["Lignite", "Oil", "Gas", "RES"]
load = ["base", "middle", "peak"]
pi = [("Lignite","base"),("Lignite","middle"),("Oil","middle"),("Oil","peak"),("Gas","base"), 
       ("Gas","middle"),("Gas","peak"),("RES","base"),("RES","peak")]  
es = ["Lignite","RES"]
k = [ "cost", "CO2emission", "endogenous"]

cost = {"Lignite":30,"Oil":75,"Gas":60,"RES":90}
capacity = {"Lignite":31000,"Oil":15000,"Gas":22000,"RES":10000}
CO2emission = { "Lignite":1.44,"Oil":0.72,"Gas":0.45,"RES":0}
demand = {"base": 38400.0, "middle": 19200.0, "peak": 6400.0}

# 第一步：计算每个目标的上下界，生成payoff表
def get_bound(obj_idx):
    mdl = Model('bound_calc')
    mdl.setParam('OutputFlag', 0) # 关闭日志输出
    # 仅在允许的pi配对中定义发电量变量
    x = mdl.addVars(pi, vtype=GRB.CONTINUOUS, name="x")
    # 约束条件
    mdl.addConstrs(quicksum(x[p,l] for (p,l) in pi if p == p_i) <= capacity[p_i] for p_i in power)
    mdl.addConstrs(quicksum(x[p,l] for (p,l) in pi if l == l_i) >= demand[l_i] for l_i in load)
    # 目标定义
    if obj_idx == 0: # 最大化清洁能源发电量（对应原模型的负目标最小化）
        obj = quicksum(-1 * x[e,l] for (e,l) in pi if e in es)
    elif obj_idx == 1: # 最小化CO2排放
        obj = quicksum(x[p,l]*CO2emission[p] for (p,l) in pi)
    else: # 最小化成本
        obj = quicksum(x[p,l]*cost[p] for (p,l) in pi)
    mdl.setObjective(obj, GRB.MINIMIZE)
    mdl.optimize()
    # 返回目标的最优值和对应的变量取值
    return mdl.objVal, {v.varName: v.x for v in mdl.getVars()}

# 获取三个目标的上下界
obj_min = [0]*3
obj_max = [0]*3
for i in range(3):
    obj_min[i], _ = get_bound(i)
    # 求目标的最大值：改为最大化该目标即可
    mdl = Model('max_calc')
    mdl.setParam('OutputFlag', 0)
    x = mdl.addVars(pi, vtype=GRB.CONTINUOUS, name="x")
    mdl.addConstrs(quicksum(x[p,l] for (p,l) in pi if p == p_i) <= capacity[p_i] for p_i in power)
    mdl.addConstrs(quicksum(x[p,l] for (p,l) in pi if l == l_i) >= demand[l_i] for l_i in load)
    if i == 0:
        obj = quicksum(-1 * x[e,l] for (e,l) in pi if e in es)
    elif i ==1:
        obj = quicksum(x[p,l]*CO2emission[p] for (p,l) in pi)
    else:
        obj = quicksum(x[p,l]*cost[p] for (p,l) in pi)
    mdl.setObjective(obj, GRB.MAXIMIZE)
    mdl.optimize()
    obj_max[i] = mdl.objVal

# 第二步：设置ε步长，迭代求解
step_num = 10 # 每个非主目标的分段数，可调整，越大结果越密耗时越长
# CO2的ε取值序列
co2_eps_list = np.linspace(obj_min[1], obj_max[1], step_num)
# 清洁能源目标的ε取值序列
es_eps_list = np.linspace(obj_min[0], obj_max[0], step_num)

all_solutions = []
for co2_eps in co2_eps_list:
    for es_eps in es_eps_list:
        mdl = Model('eps_constraint')
        mdl.setParam('OutputFlag', 0)
        x = mdl.addVars(pi, vtype=GRB.CONTINUOUS, name="x")
        # 基础约束
        mdl.addConstrs(quicksum(x[p,l] for (p,l) in pi if p == p_i) <= capacity[p_i] for p_i in power)
        mdl.addConstrs(quicksum(x[p,l] for (p,l) in pi if l == l_i) >= demand[l_i] for l_i in load)
        # ε约束：非主目标不超过阈值
        mdl.addConstr(quicksum(x[p,l]*CO2emission[p] for (p,l) in pi) <= co2_eps, name="co2_eps_con")
        mdl.addConstr(quicksum(-1 * x[e,l] for (e,l) in pi if e in es) <= es_eps, name="es_eps_con")
        # 主目标：最小化成本
        mdl.setObjective(quicksum(x[p,l]*cost[p] for (p,l) in pi), GRB.MINIMIZE)
        mdl.optimize()
        if mdl.status == GRB.OPTIMAL:
            # 保存三个目标的取值
            cost_val = mdl.objVal
            co2_val = quicksum(x[p,l].x*CO2emission[p] for (p,l) in pi).getValue()
            es_val = -1 * quicksum(x[e,l].x for (e,l) in pi if e in es).getValue()
            all_solutions.append({
                "cost": cost_val,
                "co2": co2_val,
                "es_output": es_val,
                "x": {k: v.x for k,v in x.items()}
            })

# 第三步：过滤非支配解，得到Pareto前沿
def is_non_dominated(sol, sol_list):
    for s in sol_list:
        # 存在另一个解所有目标都不差于当前解，且至少一个目标更优，则当前解是被支配的
        if s["cost"] <= sol["cost"] and s["co2"] <= sol["co2"] and s["es_output"] >= sol["es_output"]:
            if s["cost"] < sol["cost"] or s["co2"] < sol["co2"] or s["es_output"] > sol["es_output"]:
                return False
    return True

pareto_solutions = [sol for sol in all_solutions if is_non_dominated(sol, all_solutions)]

# 输出Pareto前沿的目标值
print("Pareto前沿点（成本/CO2排放/清洁能源发电量）：")
for idx, sol in enumerate(pareto_solutions):
    print(f"点{idx+1}: 成本={sol['cost']:.2f}元, CO2排放={sol['co2']:.2f}吨, 清洁能源发电量={sol['es_output']:.2f}GWh")

结果调整说明 #

• 可修改step_num参数调整Pareto点的密度，数值越大采样越密，求解耗时越长
• 可根据需求更换主优化目标，只需修改ε约束对应的目标即可
• 得到的pareto_solutions列表中存储了每个Pareto点对应的调度方案，可直接提取使用

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