Quicksum

/*
*程序的版权和版本声明部分:
*Copyright(c)2014,烟台大学计算机学院学生
*All rights reserved.
*文件名称:
*作者:田成琳
*完成日期:2014 年 6 月 12 日
*版本号:v1.0
*对任务及求解方法的描述部分:
*问题描述:A checksum is an algorithm that scans a packet of data
           and returns a single number.The idea is that if the packet is changed,
           the checksum will also change, so checksums are often used for detecting
           transmission errors, validating document contents,
           and in many other situations where it is necessary to detect
           undesirable changes in data. For this problem, you will implement a
           checksum algorithm called Quicksum. A Quicksum packet allows only
           uppercase letters and spaces. It always begins and ends with
           an uppercase letter. Otherwise, spaces and letters can occur
           in any combination, including consecutive spaces.
           A Quicksum is the sum of the products of each character's
           position in the packet times the character's value.
           A space has a value of zero,while letters have a value equal to
           their position in the alphabet.So, A=1, B=2, etc.,
           through Z=26. Here are exampleQuicksum calculations for the packets
           "ACM" and "MID CENTRAL":
           ACM: 1*1 + 2*3 + 3*13 = 46
           MID CENTRAL: 1*13+2*9+3*4+4*0+5*3+6*5+7*14+8*20+9*18+10*1+11*12=650

*程序输入:The input consists of one or more packets followed by a line
          containing only # that signals the end of the input.
          Each packet is on a line by itself, does not begin or end with a space,
          and contains from 1 to 255 characters.

*程序输出:For each packet, output its Quicksum on a separate line in the output.
*问题分析:
*算法设计:
*/

题目分析:题中意思说输入案例有多个,每个输入案例中只包含大写字母和空格,然后根据每个大写字母的位置做相应计算。A=1,B=2,C=3...,Z=26。比如ACM:A在第一位,C第二位,M第三位,那么结果就为1*1+2*3+3*13=46。但应引起注意的是案例中有的有空格字符,空格字符看成0对待。那么实现代码如下:

#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
int main()
{
    char word[100000];
    int sum,length,i;
    while(cin.getline(word,sizeof(word))&&(string)word!="#")//运用类型转换
    {
        sum=0;
        length=strlen(word);
        for(i=0; i<length; i++)
        {
            while(word[i]==' ')//忽略空格
                i++;
            sum+=(i+1)*(int)(word[i]-64);//类型转换
        }
        cout<<sum<<endl;
    }
    return 0;
}

运行结果:


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
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值