Bash 使用技巧大补贴

Bash 使用技巧大补贴

Bash 是我们经常与之打交道的 Shell 程序,本文针对其使用技巧进行了搜罗。相信在你看过这些内容之后,定会在 Bash 的世界里游刃有余。

  • 从历史中执行命令 有时候,我们需要在 Bash 中重复执行先前的命令。你当然可以使用上方向键来查看之前曾经运行过的命令。但这里有一种更好的方式:你可以按 Ctrl + r 组合键进入历史搜索模式,一旦找到需要重复执行的命令,按回车键即可。
  • 重复命令参数 先来看一个例子: mkdir /path/to/exampledir cd !$ 本例中,第一行命令将创建一个目录,而第二行的命令则转到刚创建的目录。这里,“!$”的作用就是重复前一个命令的参数。事实上,不仅是命令的参数可以重复,命令的选项同样可以。另外,Esc + . 快捷键可以切换这些命令参数或选项。
  • 用于编辑的快捷键
    • Ctrl + a:将光标定位到命令的开头
    • Ctrl + e:与上一个快捷键相反,将光标定位到命令的结尾
    • Ctrl + u:剪切光标之前的内容
    • Ctrl + k:与上一个快捷键相反,剪切光标之后的内容
    • Ctrl + y:粘贴以上两个快捷键所剪切的内容
    • Ctrl + t:交换光标之前两个字符的顺序
    • Ctrl + w:删除光标左边的参数(选项)或内容
    • Ctrl + l:清屏
  • 处理作业 首先,使用 Ctrl + z 快捷键可以让正在执行的命令挂起。如果要让该进程在后台执行,那么可以执行 bg 命令。而 fg 命令则可以让该进程重新回到前台来。使用 jobs 命令能够查看到哪些进程在后台执行。 你也可以在 fg 或 bg 命令中使用作业 id,如: fg %3 又如: bg %7
  • 使用置换
    • 命令置换 先看例子: du -h -a -c $(find . -name *.conf 2>&-) 注意 $() 中的部分,这将告诉 Bash 运行 find 命令,然后把返回的结果作为 du 的参数。
    • 进程置换 仍然先看例子: diff <(ps axo comm) <(ssh user@host ps axo comm) 该命令将比较本地系统和远程系统中正在运行的进程。请注意 <() 中的部分。
    • xargs 看例: find . -name *.conf -print0 | xargs -0 grep -l -Z mem_limit | xargs -0 -i cp {} {}.bak 该命令将备份当前目录中的所有 .conf 文件。
  • 使用管道 下面是一个简单的使用管道的例子: ps aux | grep init 这里,“|”操作符将 ps aux 的输出重定向给 grep init。 下面还有两个稍微复杂点的例子: ps aux | tee filename | grep init 及: ps aux | tee -a filename | grep init
  • 将标准输出保存为文件 你可以将命令的标准输出内容保存到一个文件中,举例如下: ps aux > filename 注意其中的“>”符号。 你也可以将这些输出内容追加到一个已存在的文件中: ps aux >> filename 你还可以分割一个较长的行: command1 | command2 | ... | commandN > tempfile1 cat tempfile1 | command1 | command2 | ... | commandN > tempfile2
  • 标准流:重定向与组合 重定向流的例子: ps aux 2>&1 | grep init 这里的数字代表:
    • 0:stdin
    • 1:stdout
    • 2:sterr
    上面的命令中,“grep init”不仅搜索“ps aux”的标准输出,而且搜索 sterr 输出。
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}秒")
最新发布
09-25
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值