import gurobipy as gp
from gurobipy import GRB
import numpy as np
import pandas as pd
import random
import time
import copy
# 记录总开始时间
total_start_time = time.time()
# 读取场景数据
scenarios_df = pd.read_csv('scenarios/scenarios_S5.csv')
scenarios = scenarios_df['Scenario'].tolist()
mu_s = {row['Scenario']: row['mu'] for _, row in scenarios_df.iterrows()}
prob_hist = {row['Scenario']: row['Probability'] for _, row in scenarios_df.iterrows()}
# 计算矩的上下限 (基于历史数据)
mu_values = np.array(list(mu_s.values()))
mu_min = mu_values.min()
mu_max = mu_values.max()
mu_avg = mu_values.mean()
sigma2_min = (mu_values ** 2).min()
sigma2_max = (mu_values ** 2).max()
sigma2_avg = (mu_values ** 2).mean()
# 设置矩约束的上下限 (适当放宽范围)
mu_lb = min(0.1, mu_avg - 0.2) # 下限
mu_ub = max(0.9, mu_avg + 0.2) # 上限
sigma2_lb = min(0.01, sigma2_avg - 0.1) # 下限
sigma2_ub = max(0.99, sigma2_avg + 0.1) # 上限
print(f"一阶矩区间: [{mu_lb:.4f}, {mu_ub:.4f}]")
print(f"二阶矩区间: [{sigma2_lb:.4f}, {sigma2_ub:.4f}]")
# 数据准备
# 设施集合
C = ['C1', 'C2'] # 二类车主
D = ['D1', 'D2'] # 一类车主
J = ['J1', 'J2'] # 维修中心
L = ['L1', 'L2'] # 回收价格选项
H = ['H1', 'H2'] # 4S店
F = ['F1', 'F2'] # 运营公司
K = ['K1', 'K2'] # 综合利用企业
G = ['G1', 'G2'] # 拆车厂
M = ['M1', 'M2'] # 梯次利用中心
N = ['N1', 'N2'] # 再生利用企业
A = ['Co', 'Mn', 'Ni', 'Li', 'Others'] # 金属材料
# 参数赋值
# 废旧动力电池单位剩余容量
L_bar = 0.7
# 回收量参数(吨)
Q1 = {'L1': 112, 'L2': 75} # 一类车主回收量
Q2 = {'L1': 112, 'L2': 75} # 二类车主回收量
Q3 = {'L1': 75, 'L2': 50} # 拆车厂回收量
# 经济参数(元)
pr = {'L1': 4000, 'L2': 6500} # 回收价格
me1, me2 = 4000, 2000 # 单位维修收益
SU = 304 # 梯次利用补贴
theta = 0.5 # 拆车厂梯次利用率
fc_j1 = {j: 200000 for j in J} # 合营建设成本
fc_j2 = {j: 500000 for j in J} # 新建建设成本
eta1, eta2 = 1, 1 / 2 # 建设模式系数
# 容量参数(吨)
cl_j = {j: 145000 for j in J} # 维修中心容量
cl_k = {k: 145000 for k in K} # 综合利用企业
cl_m = {m: 117000 for m in M} # 梯次中心
cl_n = {n: 38000 for n in N} # 再生企业
# 运营成本(元/吨)
vc_g = {g: random.randint(180, 220) for g in G}
vc_h = {h: round(random.uniform(8.3, 32.15), 2) for h in H}
vc_j = {j: random.randint(480, 520) for j in J}
vc_f = {f: round(random.uniform(8.3, 32.15), 2) for f in F}
vc_k = {k: random.randint(370, 410) for k in K}
vc_m = {m: random.randint(15, 25) for m in M}
vc_n = {n: random.randint(370, 410) for n in N}
# 运输成本(元/吨)
np.random.seed(114514)
dc_hj = {(h, j): np.random.randint(4292 / 2, 4292) for h in H for j in J}
dc_jg = {(j, g): np.random.randint(4292 / 2, 4292) for j in J for g in G}
dc_jk = {(j, k): np.random.randint(4292 / 2, 4292) for j in J for k in K}
dc_gk = {(g, k): np.random.randint(4292 / 2, 4292) for g in G for k in K}
dc_km = {(k, m): np.random.randint(4292 / 2, 4292) for k in K for m in M}
dc_kn = {(k, n): np.random.randint(4292 / 2, 4292) for k in K for n in N}
dc_dh = {(d, h): np.random.randint(4292 / 2, 4292) for d in D for h in H}
dc_cf = {(c, f): np.random.randint(4292 / 2, 4292) for c in C for f in F}
dc_jh = {(j, h): np.random.randint(4292 / 2, 4292) for j in J for h in H}
dc_jf = {(j, f): np.random.randint(4292 / 2, 4292) for j in J for f in F}
# 金属参数
v_a = {'Co': 350000, 'Mn': 6600, 'Ni': 4100, 'Li': 449000, 'Others': 25000} # 市场价格(元/吨)
alpha_a = {'Co': 0.0422 * 0.98, 'Mn': 0.0109 * 0.98, 'Ni': 0.3427 * 0.98,
'Li': 0.0419 * 0.85, 'Others': 0.08 * 0.9} # 材料使用率
r1, r2 = 0.45, 0.45 # 电池回收比率
# ====================== C&CG 算法实现 ======================
def create_master_problem(active_scenarios):
"""创建主问题模型"""
master = gp.Model("Master_Problem")
# 第一阶段决策变量
X_j = master.addVars(J, vtype=GRB.INTEGER, name="X_j", lb=0, ub=2) # 维修中心模式
X_k = master.addVars(K, vtype=GRB.BINARY, name="X_k") # 综合利用企业
X_m = master.addVars(M, vtype=GRB.BINARY, name="X_m") # 梯次中心
X_n = master.addVars(N, vtype=GRB.BINARY, name="X_n") # 再生企业
Z_l = master.addVars(L, vtype=GRB.BINARY, name="Z_l") # 价格选择
# 辅助指示变量
ind_j1 = master.addVars(J, vtype=GRB.BINARY, name="ind_j1")
ind_j2 = master.addVars(J, vtype=GRB.BINARY, name="ind_j2")
ind_j3 = master.addVars(J, vtype=GRB.BINARY, name="ind_j3")
# 添加指示约束
for j in J:
master.addConstr(ind_j1[j] + ind_j2[j] + ind_j3[j] == 1, name=f"ind_j_{j}")
master.addConstr(1 * ind_j1[j] + 2 * ind_j2[j] == X_j[j], name=f"ind_j_combine_X_j_{j}")
# 第二阶段决策变量 (仅激活的场景)
Y_dh = {};
Y_cf = {};
Y_hj = {};
Y_jh = {};
Y_fj = {};
Y_jf = {};
Y_jg = {};
Y_jk = {};
Y_gk = {};
Y_km = {};
Y_kn = {};
Y_kna = {}
for s in active_scenarios:
Y_dh[s] = master.addVars(D, H, name=f"Y_dh_{s}", lb=0)
Y_cf[s] = master.addVars(C, F, name=f"Y_cf_{s}", lb=0)
Y_hj[s] = master.addVars(H, J, name=f"Y_hj_{s}", lb=0)
Y_jh[s] = master.addVars(J, H, name=f"Y_jh_{s}", lb=0)
Y_fj[s] = master.addVars(F, J, name=f"Y_fj_{s}", lb=0)
Y_jf[s] = master.addVars(J, F, name=f"Y_jf_{s}", lb=0)
Y_jg[s] = master.addVars(J, G, name=f"Y_jg_{s}", lb=0)
Y_jk[s] = master.addVars(J, K, name=f"Y_jk_{s}", lb=0)
Y_gk[s] = master.addVars(G, K, name=f"Y_gk_{s}", lb=0)
Y_km[s] = master.addVars(K, M, name=f"Y_km_{s}", lb=0)
Y_kn[s] = master.addVars(K, N, name=f"Y_kn_{s}", lb=0)
Y_kna[s] = master.addVars(K, N, A, name=f"Y_kna_{s}", lb=0)
# 概率变量 (决策变量)
p_s = master.addVars(active_scenarios, name="p_s", lb=0)
# 第一阶段目标函数组件
FC = gp.quicksum(eta1 * fc_j1[j] * ind_j1[j] + eta2 * fc_j2[j] * ind_j2[j] for j in J)
RC = gp.quicksum(pr[l] * Z_l[l] * (Q1[l] + Q2[l]) for l in L)
stage1_cost = FC + RC
# 第二阶段目标函数组件 (仅激活的场景)
stage2_obj = {}
for s in active_scenarios:
# 维修收益
ME_s = (1 - r1) * me1 * gp.quicksum(Y_hj[s][h, j] for h in H for j in J) + \
(1 - r2) * me2 * gp.quicksum(Y_fj[s][f, j] for f in F for j in J)
# 梯次利用收益 (依赖场景的mu)
AE_s = L_bar * SU * (
mu_s[s] * gp.quicksum(Z_l[l] * (Q1[l] + Q2[l]) for l in L) +
theta * gp.quicksum(Z_l[l] * Q3[l] for l in L)
)
# 再生利用收益
YE_s = gp.quicksum(v_a[a] * alpha_a[a] * Y_kna[s][k, n, a]
for k in K for n in N for a in A)
# 运营成本
EPC_s = (
gp.quicksum(vc_h[h] * gp.quicksum(Y_dh[s][d, h] for d in D) for h in H) +
gp.quicksum(vc_f[f] * gp.quicksum(Y_cf[s][c, f] for c in C) for f in F) +
gp.quicksum(vc_j[j] * (gp.quicksum(Y_hj[s][h, j] for h in H) +
gp.quicksum(Y_fj[s][f, j] for f in F)) for j in J) +
gp.quicksum(vc_g[g] * gp.quicksum(Y_jg[s][j, g] for j in J) for g in G) +
gp.quicksum(vc_k[k] * (gp.quicksum(Y_jk[s][j, k] for j in J) +
gp.quicksum(Y_gk[s][g, k] for g in G)) for k in K) +
gp.quicksum(vc_m[m] * gp.quicksum(Y_km[s][k, m] for k in K) for m in M) +
gp.quicksum(vc_n[n] * gp.quicksum(Y_kn[s][k, n] for k in K) for n in N)
)
# 运输成本
DC_s = (
gp.quicksum(dc_dh[d, h] * Y_dh[s][d, h] for d in D for h in H) +
gp.quicksum(dc_cf[c, f] * Y_cf[s][c, f] for c in C for f in F) +
gp.quicksum(dc_hj[h, j] * (Y_hj[s][h, j] + Y_jh[s][j, h]) for h in H for j in J) +
gp.quicksum(dc_jg[j, g] * Y_jg[s][j, g] for j in J for g in G) +
gp.quicksum(dc_jk[j, k] * Y_jk[s][j, k] for j in J for k in K) +
gp.quicksum(dc_gk[g, k] * Y_gk[s][g, k] for g in G for k in K) +
gp.quicksum(dc_km[k, m] * Y_km[s][k, m] for k in K for m in M) +
gp.quicksum(dc_kn[k, n] * Y_kn[s][k, n] for k in K for n in N) +
gp.quicksum(dc_jf[j, f] * (Y_jf[s][j, f] + Y_fj[s][f, j]) for j in J for f in F)
) * (1 / 50) # 运输成本缩放因子
# 第二阶段目标函数值
stage2_obj[s] = ME_s + AE_s + YE_s - EPC_s - DC_s
# 总目标函数 (最大化系统利润)
total_obj = -stage1_cost + gp.quicksum(p_s[s] * stage2_obj[s] for s in active_scenarios)
master.setObjective(total_obj, GRB.MAXIMIZE)
# 第一阶段约束
master.addConstr(gp.quicksum(Z_l[l] for l in L) == 1, "Single_Price")
master.addConstr(gp.quicksum(X_j[j] for j in J) >= 1, "Open_J")
master.addConstr(gp.quicksum(X_k[k] for k in K) >= 1, "Open_K")
master.addConstr(gp.quicksum(X_m[m] for m in M) >= 1, "Open_M")
master.addConstr(gp.quicksum(X_n[n] for n in N) >= 1, "Open_N")
# 模糊集约束 (矩约束) - 仅激活的场景
master.addConstr(gp.quicksum(p_s[s] * mu_s[s] for s in active_scenarios) >= mu_lb, "mu_lb")
master.addConstr(gp.quicksum(p_s[s] * mu_s[s] for s in active_scenarios) <= mu_ub, "mu_ub")
master.addConstr(gp.quicksum(p_s[s] * (mu_s[s] ** 2) for s in active_scenarios) >= sigma2_lb, "sigma2_lb")
master.addConstr(gp.quicksum(p_s[s] * (mu_s[s] ** 2) for s in active_scenarios) <= sigma2_ub, "sigma2_ub")
master.addConstr(gp.quicksum(p_s[s] for s in active_scenarios) == 1, "prob_sum")
# 第二阶段约束 (仅激活的场景)
for s in active_scenarios:
# 设施最大容量约束
for j in J:
cap_j = cl_j[j] * X_j[j]
master.addConstr(gp.quicksum(Y_hj[s][h, j] for h in H) +
gp.quicksum(Y_fj[s][f, j] for f in F) <= cap_j,
name=f"Cap_J_{j}_{s}")
for k in K:
cap_k = cl_k[k] * X_k[k]
master.addConstr(gp.quicksum(Y_jk[s][j, k] for j in J) +
gp.quicksum(Y_gk[s][g, k] for g in G) <= cap_k,
name=f"Cap_K_{k}_{s}")
for m in M:
cap_m = cl_m[m] * X_m[m]
master.addConstr(gp.quicksum(Y_km[s][k, m] for k in K) <= cap_m,
name=f"Cap_M_{m}_{s}")
for n in N:
cap_n = cl_n[n] * X_n[n]
master.addConstr(gp.quicksum(Y_kn[s][k, n] for k in K) <= cap_n,
name=f"Cap_N_{n}_{s}")
# 流量平衡约束
master.addConstr(
gp.quicksum(Y_dh[s][d, h] for d in D for h in H) == gp.quicksum(Q1[l] * Z_l[l] for l in L),
name=f"Insert_D_{s}"
)
master.addConstr(
gp.quicksum(Y_cf[s][c, f] for c in C for f in F) == gp.quicksum(Q2[l] * Z_l[l] for l in L),
name=f"Insert_C_{s}"
)
master.addConstr(
gp.quicksum(Y_gk[s][g, k] for g in G for k in K) == gp.quicksum(Q3[l] * Z_l[l] for l in L),
name=f"Insert_G_{s}"
)
# H流量平衡
for h in H:
master.addConstr(
gp.quicksum(Y_dh[s][d, h] for d in D) == gp.quicksum(Y_hj[s][h, j] for j in J),
name=f"Balance_H_{h}_{s}"
)
# HJ之间关系
for h in H:
for j in J:
master.addConstr(
Y_jh[s][j, h] == (1 - r1) * Y_hj[s][h, j],
name=f"Balance_DH_{h}_{j}_{s}"
)
# F流量平衡
for f in F:
master.addConstr(
gp.quicksum(Y_cf[s][c, f] for c in C) == gp.quicksum(Y_fj[s][f, j] for j in J),
name=f"Balance_F_{f}_{s}"
)
# JF之间关系
for f in F:
for j in J:
master.addConstr(
Y_jf[s][j, f] == (1 - r2) * Y_fj[s][f, j],
name=f"Balance_FJ_{f}_{j}_{s}"
)
# 拆车厂G流量平衡
for g in G:
master.addConstr(
gp.quicksum(Y_jg[s][j, g] for j in J) == gp.quicksum(Y_gk[s][g, k] for k in K),
name=f"Balance_G_{g}_{s}"
)
# 维修中心J流量平衡
for j in J:
master.addConstr(
gp.quicksum(Y_hj[s][h, j] for h in H) - gp.quicksum(Y_jh[s][j, h] for h in H) +
gp.quicksum(Y_fj[s][f, j] for f in F) - gp.quicksum(Y_jf[s][j, f] for f in F) ==
gp.quicksum(Y_jg[s][j, g] for g in G) + gp.quicksum(Y_jk[s][j, k] for k in K),
name=f"Balance_J_{j}_{s}"
)
# 综合利用企业K流量平衡
for k in K:
master.addConstr(
gp.quicksum(Y_jk[s][j, k] for j in J) + gp.quicksum(Y_gk[s][g, k] for g in G) ==
gp.quicksum(Y_km[s][k, m] for m in M) + gp.quicksum(Y_kn[s][k, n] for n in N),
name=f"Balance_K_{k}_{s}"
)
# K和M之间的流量 (依赖场景的mu_s)
for k in K:
master.addConstr(
gp.quicksum(Y_km[s][k, m] for m in M) ==
mu_s[s] * gp.quicksum(Y_jk[s][j, k] for j in J) +
theta * gp.quicksum(Y_gk[s][g, k] for g in G),
name=f"Balance_K&M_{k}_{s}"
)
# 再生企业N流量平衡
for k in K:
for n in N:
master.addConstr(
Y_kn[s][k, n] == gp.quicksum(Y_kna[s][k, n, a] for a in A),
name=f"Metal_Balance_{k}_{n}_{s}"
)
# 设置求解参数
master.Params.OutputFlag = 0
master.Params.Cuts = 3 # 最激进的切割生成
master.Params.Presolve = 2 # 加强预处理 (0-2)
master.Params.Heuristics = 0.8 # 启发式搜索比例 (0-1)
master.Params.MIPFocus = 1 # 1=可行解, 2=最优性证明, 3=边界改进
master.Params.NonConvex = 2 # 处理非线性目标
master.Params.MIPGap = 0.05 # 1%的MIP间隙
master.Params.ConcurrentMIP = 4 # 并行求解4个MIP模型
return master, X_j, X_k, X_m, X_n, Z_l, p_s
def create_subproblem(X_j_val, X_k_val, X_m_val, X_n_val, Z_l_val, candidate_scenarios):
"""创建子问题模型,用于找到最坏情况场景"""
sub = gp.Model("Subproblem")
# 第二阶段决策变量 (仅候选场景)
Y_dh = sub.addVars(candidate_scenarios, D, H, name="Y_dh", lb=0)
Y_cf = sub.addVars(candidate_scenarios, C, F, name="Y_cf", lb=0)
Y_hj = sub.addVars(candidate_scenarios, H, J, name="Y_hj", lb=0)
Y_jh = sub.addVars(candidate_scenarios, J, H, name="Y_jh", lb=0)
Y_fj = sub.addVars(candidate_scenarios, F, J, name="Y_fj", lb=0)
Y_jf = sub.addVars(candidate_scenarios, J, F, name="Y_jf", lb=0)
Y_jg = sub.addVars(candidate_scenarios, J, G, name="Y_jg", lb=0)
Y_jk = sub.addVars(candidate_scenarios, J, K, name="Y_jk", lb=0)
Y_gk = sub.addVars(candidate_scenarios, G, K, name="Y_gk", lb=0)
Y_km = sub.addVars(candidate_scenarios, K, M, name="Y_km", lb=0)
Y_kn = sub.addVars(candidate_scenarios, K, N, name="Y_kn", lb=0)
# 概率变量 (决策变量)
p_s = sub.addVars(candidate_scenarios, name="p_s", lb=0)
# 第二阶段目标函数组件
# 计算平均金属价值(常数)
v_avg = sum(v_a[a] * alpha_a[a] for a in A) / len(A)
stage2_obj = {}
for s in candidate_scenarios:
# 维修收益
ME_s = (1 - r1) * me1 * gp.quicksum(Y_hj[s, h, j] for h in H for j in J) + \
(1 - r2) * me2 * gp.quicksum(Y_fj[s, f, j] for f in F for j in J)
# 梯次利用收益 (依赖场景的mu)
AE_s = L_bar * SU * (
mu_s[s] * gp.quicksum(Z_l_val[l] * (Q1[l] + Q2[l]) for l in L) +
theta * gp.quicksum(Z_l_val[l] * Q3[l] for l in L)
)
# 再生利用收益
YE_s = v_avg * gp.quicksum(Y_kn[s, k, n] for k in K for n in N)
# 运营成本
EPC_s = (
gp.quicksum(vc_h[h] * gp.quicksum(Y_dh[s, d, h] for d in D) for h in H) +
gp.quicksum(vc_f[f] * gp.quicksum(Y_cf[s, c, f] for c in C) for f in F) +
gp.quicksum(vc_j[j] * (gp.quicksum(Y_hj[s, h, j] for h in H) +
gp.quicksum(Y_fj[s, f, j] for f in F)) for j in J) +
gp.quicksum(vc_g[g] * gp.quicksum(Y_jg[s, j, g] for j in J) for g in G) +
gp.quicksum(vc_k[k] * (gp.quicksum(Y_jk[s, j, k] for j in J) +
gp.quicksum(Y_gk[s, g, k] for g in G)) for k in K) +
gp.quicksum(vc_m[m] * gp.quicksum(Y_km[s, k, m] for k in K) for m in M) +
gp.quicksum(vc_n[n] * gp.quicksum(Y_kn[s, k, n] for k in K) for n in N)
)
# 运输成本
DC_s = (
gp.quicksum(dc_dh[d, h] * Y_dh[s, d, h] for d in D for h in H) +
gp.quicksum(dc_cf[c, f] * Y_cf[s, c, f] for c in C for f in F) +
gp.quicksum(dc_hj[h, j] * (Y_hj[s, h, j] + Y_jh[s, j, h]) for h in H for j in J) +
gp.quicksum(dc_jg[j, g] * Y_jg[s, j, g] for j in J for g in G) +
gp.quicksum(dc_jk[j, k] * Y_jk[s, j, k] for j in J for k in K) +
gp.quicksum(dc_gk[g, k] * Y_gk[s, g, k] for g in G for k in K) +
gp.quicksum(dc_km[k, m] * Y_km[s, k, m] for k in K for m in M) +
gp.quicksum(dc_kn[k, n] * Y_kn[s, k, n] for k in K for n in N) +
gp.quicksum(dc_jf[j, f] * (Y_jf[s, j, f] + Y_fj[s, f, j]) for j in J for f in F)
) * (1 / 100) # 运输成本缩放因子
# 第二阶段目标函数值
stage2_obj[s] = ME_s + AE_s + YE_s - EPC_s - DC_s
# 总目标函数 (最大化期望利润)
sub.setObjective(gp.quicksum(p_s[s] * stage2_obj[s] for s in candidate_scenarios), GRB.MAXIMIZE)
# 模糊集约束 (矩约束)
sub.addConstr(gp.quicksum(p_s[s] * mu_s[s] for s in candidate_scenarios) >= mu_lb, "mu_lb")
sub.addConstr(gp.quicksum(p_s[s] * mu_s[s] for s in candidate_scenarios) <= mu_ub, "mu_ub")
sub.addConstr(gp.quicksum(p_s[s] * (mu_s[s] ** 2) for s in candidate_scenarios) >= sigma2_lb, "sigma2_lb")
sub.addConstr(gp.quicksum(p_s[s] * (mu_s[s] ** 2) for s in candidate_scenarios) <= sigma2_ub, "sigma2_ub")
sub.addConstr(gp.quicksum(p_s[s] for s in candidate_scenarios) == 1, "prob_sum")
# 第二阶段约束 (每个候选场景)
for s in candidate_scenarios:
# 设施最大容量约束
for j in J:
cap_j = cl_j[j] * X_j_val[j]
sub.addConstr(gp.quicksum(Y_hj[s, h, j] for h in H) +
gp.quicksum(Y_fj[s, f, j] for f in F) <= cap_j,
name=f"Cap_J_{j}_{s}")
for k in K:
cap_k = cl_k[k] * X_k_val[k]
sub.addConstr(gp.quicksum(Y_jk[s, j, k] for j in J) +
gp.quicksum(Y_gk[s, g, k] for g in G) <= cap_k,
name=f"Cap_K_{k}_{s}")
for m in M:
cap_m = cl_m[m] * X_m_val[m]
sub.addConstr(gp.quicksum(Y_km[s, k, m] for k in K) <= cap_m,
name=f"Cap_M_{m}_{s}")
for n in N:
cap_n = cl_n[n] * X_n_val[n]
sub.addConstr(gp.quicksum(Y_kn[s, k, n] for k in K) <= cap_n,
name=f"Cap_N_{n}_{s}")
# 流量平衡约束
sub.addConstr(
gp.quicksum(Y_dh[s, d, h] for d in D for h in H) == gp.quicksum(Q1[l] * Z_l_val[l] for l in L),
name=f"Insert_D_{s}"
)
sub.addConstr(
gp.quicksum(Y_cf[s, c, f] for c in C for f in F) == gp.quicksum(Q2[l] * Z_l_val[l] for l in L),
name=f"Insert_C_{s}"
)
sub.addConstr(
gp.quicksum(Y_gk[s, g, k] for g in G for k in K) == gp.quicksum(Q3[l] * Z_l_val[l] for l in L),
name=f"Insert_G_{s}"
)
# H流量平衡
for h in H:
sub.addConstr(
gp.quicksum(Y_dh[s, d, h] for d in D) == gp.quicksum(Y_hj[s, h, j] for j in J),
name=f"Balance_H_{h}_{s}"
)
# HJ之间关系
for h in H:
for j in J:
sub.addConstr(
Y_jh[s, j, h] == (1 - r1) * Y_hj[s, h, j],
name=f"Balance_DH_{h}_{j}_{s}"
)
# F流量平衡
for f in F:
sub.addConstr(
gp.quicksum(Y_cf[s, c, f] for c in C) == gp.quicksum(Y_fj[s, f, j] for j in J),
name=f"Balance_F_{f}_{s}"
)
# JF之间关系
for f in F:
for j in J:
sub.addConstr(
Y_jf[s, j, f] == (1 - r2) * Y_fj[s, f, j],
name=f"Balance_FJ_{f}_{j}_{s}"
)
# 拆车厂G流量平衡
for g in G:
sub.addConstr(
gp.quicksum(Y_jg[s, j, g] for j in J) == gp.quicksum(Y_gk[s, g, k] for k in K),
name=f"Balance_G_{g}_{s}"
)
# 维修中心J流量平衡
for j in J:
sub.addConstr(
gp.quicksum(Y_hj[s, h, j] for h in H) - gp.quicksum(Y_jh[s, j, h] for h in H) +
gp.quicksum(Y_fj[s, f, j] for f in F) - gp.quicksum(Y_jf[s, j, f] for f in F) ==
gp.quicksum(Y_jg[s, j, g] for g in G) + gp.quicksum(Y_jk[s, j, k] for k in K),
name=f"Balance_J_{j}_{s}"
)
# 综合利用企业K流量平衡
for k in K:
sub.addConstr(
gp.quicksum(Y_jk[s, j, k] for j in J) + gp.quicksum(Y_gk[s, g, k] for g in G) ==
gp.quicksum(Y_km[s, k, m] for m in M) + gp.quicksum(Y_kn[s, k, n] for n in N),
name=f"Balance_K_{k}_{s}"
)
# K和M之间的流量 (依赖场景的mu_s)
for k in K:
sub.addConstr(
gp.quicksum(Y_km[s, k, m] for m in M) ==
mu_s[s] * gp.quicksum(Y_jk[s, j, k] for j in J) +
theta * gp.quicksum(Y_gk[s, g, k] for g in G),
name=f"Balance_K&M_{k}_{s}"
)
# 添加再生企业容量约束
for n in N:
cap_n = cl_n[n] * X_n_val[n]
sub.addConstr(gp.quicksum(Y_kn[s, k, n] for k in K) <= cap_n,
name=f"Cap_N_{n}_{s}")
# 设置求解参数
sub.Params.OutputFlag = 0
sub.Params.NonConvex = 2 # 允许求解非凸二次规划问题
sub.Params.Method = 2 # 内点法(更适合QP问题)
sub.Params.NumericFocus = 3
sub.Params.ScaleFlag = 2
sub.Params.BarHomogeneous = 1
sub.Params.BarQCPConvTol = 1e-5
return sub, p_s
def ccg_algorithm():
"""C&CG算法主函数"""
# 初始化 - 添加一个初始场景确保可行性
active_scenarios = scenarios[:min(3, len(scenarios))] # 初始添加第一个场景
all_scenarios = scenarios
best_sol = None
LB = -GRB.INFINITY
UB = GRB.INFINITY
tol = 1e3
max_iter = 10
iter_count = 0
failure_count = 0
best_UB = -GRB.INFINITY # 添加最佳上界跟踪
print("\n========== 开始C&CG算法 ==========")
print(f"最大迭代次数: {max_iter}, 收敛容差: {tol}")
# 初始主问题(包含一个场景)
master, X_j, X_k, X_m, X_n, Z_l, p_s = create_master_problem(active_scenarios)
# 添加场景去重辅助函数
def get_available_scenarios():
return [s for s in all_scenarios if s not in active_scenarios]
# 添加设施容量可行性约束
master.addConstr(
gp.quicksum(cl_j[j] * X_j[j] for j in J) >=
max(Q1[l] + Q2[l] for l in L),
"Min_Capacity_J"
)
master.addConstr(
gp.quicksum(cl_k[k] * X_k[k] for k in K) >=
max(r1 * Q1[l] + r2 * Q2[l] + Q3[l] for l in L),
"Min_Capacity_K"
)
while iter_count < max_iter and (UB - LB) > tol:
iter_count += 1
print(f"\n--- 迭代 {iter_count} ---")
print(f"激活场景数: {len(active_scenarios)}/{len(all_scenarios)}")
print(f"当前激活场景: {active_scenarios}")
# 步骤1: 求解主问题
print("求解主问题...")
master_start = time.time()
master.optimize()
master_time = time.time() - master_start
if master.status == GRB.OPTIMAL:
# 获取第一阶段决策值
X_j_val = {j: round(X_j[j].X) for j in J}
X_k_val = {k: round(X_k[k].X) for k in K}
X_m_val = {m: round(X_m[m].X) for m in M}
X_n_val = {n: round(X_n[n].X) for n in N}
Z_l_val = {l: Z_l[l].X for l in L}
# 计算第一阶段成本
FC_val = 0
for j in J:
mode = X_j_val[j]
if mode == 1:
FC_val += eta1 * fc_j1[j]
elif mode == 2:
FC_val += eta2 * fc_j2[j]
RC_val = 0
for l in L:
if Z_l_val[l] > 0.5:
RC_val += pr[l] * (Q1[l] + Q2[l])
stage1_cost_val = FC_val + RC_val
# 主问题目标值 (下界)
master_obj = master.ObjVal
LB = max(LB, master_obj)
print(f"主问题求解时间: {master_time:.2f}秒")
print(f"第一阶段决策: X_j={X_j_val}, X_k={X_k_val}, X_m={X_m_val}, X_n={X_n_val}")
print(f"价格选择: Z_l={Z_l_val}")
print(f"第一阶段成本: {stage1_cost_val:,.2f}元")
print(f"主问题目标值: {master_obj:,.2f}元 (当前下界: {LB:,.2f}元)")
# 获取可用场景列表
candidate_scenarios = get_available_scenarios()
# 步骤2: 求解子问题(找到最坏情况场景)
if not candidate_scenarios:
print("所有场景已激活,算法收敛")
UB = LB # 上下界收敛
# 保存当前解作为最优解
best_sol = {
'X_j': X_j_val,
'X_k': X_k_val,
'X_m': X_m_val,
'X_n': X_n_val,
'Z_l': Z_l_val,
'obj': master_obj,
'sub_obj': master_obj + stage1_cost_val,
'active_scenarios': active_scenarios.copy()
}
break
print(f"求解子问题(候选场景数: {len(candidate_scenarios)})...")
sub_start = time.time()
sub, p_s_sub = create_subproblem(X_j_val, X_k_val, X_m_val, X_n_val, Z_l_val, candidate_scenarios)
sub.optimize()
sub_time = time.time() - sub_start
if sub.status == GRB.OPTIMAL:
# 找到最坏情况场景
worst_scenario = None
max_prob = 0
p_s_vals = {}
for s in candidate_scenarios:
p_val = p_s_sub[s].X
p_s_vals[s] = p_val
if p_val > max_prob:
max_prob = p_val
worst_scenario = s
if worst_scenario:
print(f"添加最坏情况场景: {worst_scenario} (概率: {max_prob:.4f})")
active_scenarios.append(worst_scenario)
failure_count = 0 # 重置失败计数器
# 计算子问题目标值 (上界候选)
sub_obj = sub.ObjVal
total_obj_estimate = -stage1_cost_val + sub_obj
# 更新上界
if total_obj_estimate > best_UB:
best_UB = total_obj_estimate # 更新最佳上界
best_sol = {
'X_j': X_j_val,
'X_k': X_k_val,
'X_m': X_m_val,
'X_n': X_n_val,
'Z_l': Z_l_val,
'obj': total_obj_estimate,
'sub_obj': sub_obj,
'active_scenarios': copy.copy(active_scenarios)
}
UB = best_UB # 上界使用最佳上界值
print(f"子问题求解时间: {sub_time:.2f}秒")
print(f"第二阶段利润: {sub_obj:,.2f}元")
print(f"总目标估计值: {total_obj_estimate:,.2f}元 (当前上界: {UB:,.2f}元)")
# 计算最优间隙
gap = UB - LB
gap_percent = 100 * gap / abs(LB + 1e-6)
print(f"最优间隙: {gap:,.2f}元 ({gap_percent:.2f}%)")
# 步骤3: 更新主问题,添加新场景
print("更新主问题,添加新场景...")
master, X_j, X_k, X_m, X_n, Z_l, p_s = create_master_problem(active_scenarios)
# 添加设施容量可行性约束
master.addConstr(
gp.quicksum(cl_j[j] * X_j[j] for j in J) >=
max(Q1[l] + Q2[l] for l in L),
"Min_Capacity_J"
)
master.addConstr(
gp.quicksum(cl_k[k] * X_k[k] for k in K) >=
max(r1 * Q1[l] + r2 * Q2[l] + Q3[l] for l in L),
"Min_Capacity_K"
)
else:
print("未找到最坏情况场景,算法收敛")
UB = LB # 上下界收敛
break
else:
print(f"子问题求解失败,状态: {sub.status}")
failure_count += 1
available_scenarios = get_available_scenarios()
if available_scenarios:
# 选择第一个可用场景(避免随机重复)
new_scenario = available_scenarios[0]
print(f"添加可用场景: {new_scenario}")
active_scenarios.append(new_scenario)
master, X_j, X_k, X_m, X_n, Z_l, p_s = create_master_problem(active_scenarios)
# 添加设施容量可行性约束
master.addConstr(
gp.quicksum(cl_j[j] * X_j[j] for j in J) >=
max(Q1[l] + Q2[l] for l in L),
"Min_Capacity_J"
)
master.addConstr(
gp.quicksum(cl_k[k] * X_k[k] for k in K) >=
max(r1 * Q1[l] + r2 * Q2[l] + Q3[l] for l in L),
"Min_Capacity_K"
)
else:
print("无更多场景可添加,算法终止")
break
else:
print("主问题求解失败")
failure_count += 1
# 获取可用场景
available_scenarios = get_available_scenarios()
if available_scenarios:
# 选择第一个可用场景(避免随机重复)
new_scenario = available_scenarios[0]
print(f"添加可用场景: {new_scenario}")
active_scenarios.append(new_scenario)
master, X_j, X_k, X_m, X_n, Z_l, p_s = create_master_problem(active_scenarios)
# 添加设施容量可行性约束
master.addConstr(
gp.quicksum(cl_j[j] * X_j[j] for j in J) >=
max(Q1[l] + Q2[l] for l in L),
"Min_Capacity_J"
)
master.addConstr(
gp.quicksum(cl_k[k] * X_k[k] for k in K) >=
max(r1 * Q1[l] + r2 * Q2[l] + Q3[l] for l in L),
"Min_Capacity_K"
)
else:
print("无更多场景可添加,算法终止")
break
# 检查收敛
gap = UB - LB
if gap <= tol:
print(f"\n*** 收敛于迭代 {iter_count} ***")
print(f"下界: {LB:,.2f}元, 上界: {UB:,.2f}元, 最优间隙: {gap:,.2f}元")
break
# 最终结果
if best_sol:
print("\n*** 找到最优解 ***")
elif iter_count >= max_iter:
print("\n*** 达到最大迭代次数 ***")
if master.status == GRB.OPTIMAL:
# 保存当前解作为可行解
best_sol = {
'X_j': X_j_val,
'X_k': X_k_val,
'X_m': X_m_val,
'X_n': X_n_val,
'Z_l': Z_l_val,
'obj': master_obj,
'sub_obj': master_obj + stage1_cost_val,
'active_scenarios': active_scenarios.copy()
}
print(f"找到可行解,目标值: {best_sol['obj']:,.2f}元")
else:
print("\n*** 未找到可行解 ***")
return best_sol
# 执行C&CG算法
best_solution = ccg_algorithm()
# 结果分析(使用最佳解)
if best_solution:
# 创建完整主问题模型获取详细结果
print("构建完整主问题...")
master, X_j, X_k, X_m, X_n, Z_l, p_s = create_master_problem(best_solution['active_scenarios'])
master.optimize()
if master.status == GRB.OPTIMAL:
# 提取概率分布
p_s_val = {s: p_s[s].X for s in best_solution['active_scenarios']}
# 提取关键物流变量(以第一个激活场景为例)
s0 = best_solution['active_scenarios'][0]
Y_hj_val = {}
for h in H:
for j in J:
var = master.getVarByName(f"Y_hj_{s0}[{h},{j}]")
if var:
Y_hj_val[(h, j)] = var.X
Y_jg_val = {}
for j in J:
for g in G:
var = master.getVarByName(f"Y_jg_{s0}[{j},{g}]")
if var:
Y_jg_val[(j, g)] = var.X
Y_km_val = {}
for k in K:
for m in M:
var = master.getVarByName(f"Y_km_{s0}[{k},{m}]")
if var:
Y_km_val[(k, m)] = var.X
# 打印结果
print("\n==== 最优建设方案 ====")
print("\n维修中心建设模式:")
for j in J:
mode = "未建设" if best_solution['X_j'][j] < 1 else ("合营" if best_solution['X_j'][j] == 1 else "新建")
print(f"{j}: {mode}")
# 价格选择
for l in L:
if best_solution['Z_l'][l] > 0.5:
print(f"\n选择回收价格: {pr[l]}元 (选项{l})")
# 概率分布
print("\n最坏情况概率分布:")
for s in best_solution['active_scenarios']:
if p_s_val.get(s, 0) > 0.01: # 只显示概率大于1%的场景
print(f"场景{s}: mu={mu_s[s]:.4f}, p={p_s_val[s]:.4f}")
# 关键物流量
print(f"\n=== 物流量详情 (场景{s0}) ===")
print("\n[4S店 -> 维修中心] (Y_hj):")
for (h, j), val in Y_hj_val.items():
if val > 1e-6:
print(f" 4S店 {h} -> 维修中心 {j}: {val:.2f}")
print("\n[维修中心 -> 拆车厂] (Y_jg):")
for (j, g), val in Y_jg_val.items():
if val > 1e-6:
print(f" 维修中心 {j} -> 拆车厂 {g}: {val:.2f}")
print("\n[综合利用 -> 梯次利用] (Y_km):")
for (k, m), val in Y_km_val.items():
if val > 1e-6:
print(f" 综合利用 {k} -> 梯次利用 {m}: {val:.2f}")
# 目标函数值
print("\n=== 目标函数分解 ===")
FC_val = 0
for j in J:
mode = best_solution['X_j'][j]
if mode == 1:
FC_val += eta1 * fc_j1[j]
elif mode == 2:
FC_val += eta2 * fc_j2[j]
RC_val = 0
for l in L:
if best_solution['Z_l'][l] > 0.5:
RC_val += pr[l] * (Q1[l] + Q2[l])
stage1_cost_val = FC_val + RC_val
print(f"第一阶段成本 (FC + RC): {stage1_cost_val:,.2f}元")
print(f"第二阶段期望利润: {best_solution['sub_obj']:,.2f}元")
print(f"总系统利润: {best_solution['obj']:,.2f}元")
# 矩约束验证
exp_mu = sum(p_s_val.get(s, 0) * mu_s[s] for s in best_solution['active_scenarios'])
exp_sigma2 = sum(p_s_val.get(s, 0) * (mu_s[s] ** 2) for s in best_solution['active_scenarios'])
print(f"\n验证矩约束:")
print(f"一阶矩期望值: {exp_mu:.4f} (区间: [{mu_lb:.4f}, {mu_ub:.4f}])")
print(f"二阶矩期望值: {exp_sigma2:.4f} (区间: [{sigma2_lb:.4f}, {sigma2_ub:.4f}])")
else:
print(f"完整主问题求解失败,状态: {master.status}")
else:
print("未找到最优解")
# 计算总时间
total_time = time.time() - total_start_time
print(f"\n总运行时间: {total_time:.2f}秒")
最新发布