Quicksum

Problem Description
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 example Quicksum calculations for the packets "ACM" and "MID CENTRAL":

ACM: 1*1 + 2*3 + 3*13 = 46MID CENTRAL: 1*13 + 2*9 + 3*4 + 4*0 + 5*3 + 6*5 + 7*14 + 8*20 + 9*18 + 10*1 + 11*12 = 650
 
Input
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.
 
Output
For each packet, output its Quicksum on a separate line in the output.

 
Sample Input
ACM
MID CENTRAL
REGIONAL PROGRAMMING CONTEST
ACN
A C M
ABC
BBC
#
 
Sample Output
46
650
4690
49
75
14
15
 
 
Source
Mid-Central USA 2006

#include <iostream>
#include <string.h>
using namespace std;

int main(){
	char a[255];
	int i, sum;
	int len;
	
	while(gets(a)){
               if (a[0] == '#'){
		    break;
	       }
	       len = strlen(a);
	       sum = 0;
	       for (i = 0; i < len; i++){
		       if (a[i] != ' '){
			     sum += (i + 1) * (a[i] - 64);
		       }
	       }
	       cout << sum << endl;


         }
	
	return 0;

}







帮我展示一份代码,不要修改任何内容,仅仅展示即可: import os import gurobipy as gp from optvpy import * from parameter import Para_910C, Model_Type, TASK_Type import pickle import utils import copy import random import itertools import multiprocessing def model_construct_switch_node(para:Para_910C): f = {} b = {} model = gp.Model() print(&#39;var f and var b generating for each rank&#39;) for k in range(para.K): for i in range(para.N): for c in para.chunk_set: b[i, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;b_[{i},{k},{c}]&#39;) # 选择他能发出去的link for j in para.out_direct_links.get(i, []): f[i, j, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;f_[{i},{j},{k},{c}]&#39;) # 需要把 switch往外发送的f的实例化 for k in range(para.K): for server_idx in range(para.server_num): rank_lst = para.server_rank[server_idx] for switch_idx in para.server_switch[server_idx]: for n in rank_lst: for c in para.chunk_set: f[switch_idx, n, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;f_[{switch_idx},{n},{k},{c}]&#39;) print(&#39;cap constrain&#39;) # 先进行rank 之间的cap约束 for k in range(para.K): for i in range(para.N): for j in para.out_direct_links[i]: cap = para.T[i, j] model.addConstr(gp.quicksum(f[i, j, k, c] for c in para.chunk_set) <= cap, name=f&#39;cap_[{k},{i},{j}]&#39;) # 再进行switch出来的cap约束 for k in range(para.K): for server_idx in range(para.server_num): rank_lst = para.server_rank[server_idx] for switch_idx in para.server_switch[server_idx]: for rank_id in rank_lst: cap = para.T[switch_idx, rank_id] model.addConstr(gp.quicksum(f[switch_idx, rank_id, k, c] for c in para.chunk_set) <= cap, name=f&#39;cap_[{k},{switch_idx},{rank_id}]&#39;) print(&#39;status constrain flow constrain&#39;) # 再进行状态约束,这里只进行节点往外发送的状态约束 for k in range(para.K): for i in range(para.N): init_chunk_set = para.rank_chunk_dict[i] for c in para.chunk_set: # 确定前一个时刻的状态从哪里取值 tmp_var = None if k == 0 and c in init_chunk_set: # 在0时刻上,且c属于原始chunk的,给f的右端项一个1,其余的直接从b中get值出来 tmp_var = 1 else: tmp_var = b.get((i, k - 1, c), 0) # 发送约束 for target_j in para.out_direct_links[i]: model.addConstr(f[i, target_j, k, c] <= tmp_var, name=f&#39;status_[{k},{i},{target_j},{c}]&#39;) # 状态平衡 in_expr = gp.quicksum( f.get((source_j, i, k - para.delta[source_j, i], c), 0) for source_j in para.in_direct_links[i]) model.addConstr(b[i, k, c] == tmp_var + in_expr, name=f&#39;flow_[{k},{i},{c}]&#39;) # switch cons # todo 可以考虑加这样一个约束,对于c来说,当rank1能不能接收switch接收c取决于rank2-7有没有c for server_idx in range(para.server_num): switch_lst = para.server_switch[server_idx] rank_lst = para.server_rank[server_idx] for k in range(para.K): for c in para.chunk_set: for switch_idx in switch_lst: send_expr = gp.quicksum([f.get((rank_id, switch_idx, k, c), 0) for rank_id in rank_lst]) recv_expr = gp.quicksum([f.get((switch_idx, rank_id, k, c), 0) for rank_id in rank_lst]) model.addConstr(send_expr == recv_expr, name = f&#39;switch_[{switch_idx},{k},{c}]&#39;) for rank_lst_idx in range(0, len(rank_lst), 2): rank_0 = rank_lst[rank_lst_idx] rank_1 = rank_lst[rank_lst_idx + 1] other_rank_lst = [rank_id for rank_id in rank_lst if rank_id != rank_0 and rank_id!= rank_1] # 对于rank1来说,其是否能从每一个switch中接收c都取决于 other rank有没有发送c for switch_idx in switch_lst: other_rank_send_expr = gp.quicksum([f.get((other_rank, switch_idx, k, c), 0) for other_rank in other_rank_lst]) model.addConstr(f[switch_idx, rank_0, k, c] <= other_rank_send_expr, name =f&#39;control_[{k},{c},{switch_idx},{rank_0}]&#39;) model.addConstr(f[switch_idx, rank_1, k, c] <= other_rank_send_expr, name =f&#39;control_[{k},{c},{switch_idx},{rank_1}]&#39;) # objective # all2all 的 obj设计上应该考虑 将chunk数量均分,然后给每一个chunk分配好对应的位置 expr = gp.LinExpr() for rank_id, cur_chunk_lst in para.final_obj_dict.items(): for c in cur_chunk_lst: for k in range(para.K): for source_j in para.in_direct_links[rank_id]: if k + para.delta[source_j, rank_id] >= para.K: continue expr += (para.K - k + 1) * f[source_j, rank_id, k, c] model.setObjective(expr, gp.GRB.MAXIMIZE) model.optimize() # model.write(&#39;tmp.lp&#39;) # env = OPTVEnv() # model_optv = OPTVModel(env) # model_optv.Read(&#39;tmp.lp&#39;) # model_optv.ComputeIIS() f_res = {key: val.X for key, val in f.items() if val.X} return f_res, model.Runtime def model_construct_hyper_edge_no_sio_switch_linear(para:Para_910C, K = None, initial_rank_dict = None, final_obj_dict = None): # 对于rolling horizon场景下适用 # initial rank chunk && final obj dict print("linear") if not initial_rank_dict: initial_rank_dict = copy.deepcopy(para.rank_chunk_dict) if not final_obj_dict: final_obj_dict = copy.deepcopy(para.final_obj_dict) if not K: K = para.K f = {} b = {} model = gp.Model() print(&#39;var f and var b generating for each rank&#39;) for k in range(K): for i in para.all_rank_set: # para.all_rank_set={0,3,16,19} for c in para.chunk_set: # {0, 1, 2, 3, 12, 13, 14, 15, 64, 65, 66, 67, 76, 77, 78, 79} b[i, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;b_({i},{k},{c})&#39;) # 选择他能发出去的link for j in para.out_direct_links.get(i, []): # para.out_direct_links={0: [16, 3, 19], 3: [0, 16, 19], 16: [0, 3, 19], 19: [0, 16, 3]}。也就是rank需要发送到的rank f[i, j, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;f_({i},{j},{k},{c})&#39;) print(&#39;hyper edge band upper bound&#39;) # 对于每一个rank来说,其给switch发出的大小受总量的限制 for k in range(K): for i in para.all_rank_set: # hccs 发送固有上限 hyper_links = [hyper_link for hyper_link in para.out_direct_links[i] if hyper_link != para.pair_dict.get(i)] ##!!pair_dict={},是空的,或许是遗留问题?out_direct_links里本来就没有直连的 model.addConstr(gp.quicksum(f[i, j, k, c] for c in para.chunk_set for j in hyper_links) <= para.hyper_bound, name=f&#39;f_jc_hccs_({i},{k})&#39;) # hccscap约束 # sio 发送固有上限 # 可能存在没有sio链路的情况 ## 确实是,目前没有sio链路,都是hyperedge发送。 nei_rank = para.pair_dict.get(i) if nei_rank: model.addConstr(gp.quicksum(f[i, nei_rank, k, c] for c in para.chunk_set) <= para.T[i, nei_rank], name=f&#39;f_jc_sio_({i},{nei_rank},{k})&#39;) # 目前T里面都是0 print(&#39;status constrain flow constrain&#39;) # 再进行状态约束,这里只进行节点往外发送的状态约束 for k in range(K): for i in para.all_rank_set: init_chunk_set = initial_rank_dict[i] # 这个是i rank自己的chunk for c in para.chunk_set: # 总共需要发送or接受的chunk,也就是所有会使用的到的chunk # 确定前一个时刻的状态从哪里取值 tmp_var = None if k == 0 and c in init_chunk_set: # 在0时刻上,且c属于原始chunk的,给f的右端项一个1,其余的直接从b中get值出来 tmp_var = 1 # 代表初始的b,也就是可以发送 else: tmp_var = b.get((i, k - 1, c), 0) # 这里的get 0,有点危险。k-1时刻有,k时刻才能发送 # 发送约束 ### 对于一个i rank,他的j rank从para.out_direct_links[i]取就好了。 for target_j in para.out_direct_links[i]: model.addConstr(f[i, target_j, k, c] <= tmp_var, name=f&#39;status_({i},{k},{target_j},{c})&#39;) # 状态平衡 ### 我有一个疑问,这里没有过滤c,也就是说source的rank可能给i发送他自己的c?虽然是这样子不会出错。 in_expr = gp.quicksum( f.get((source_j, i, k - para.delta[source_j, i], c), 0) for source_j in para.in_direct_links[i]) model.addConstr(b[i, k, c] == tmp_var + in_expr, name=f&#39;flow_({i},{k},{c})&#39;) # 给每一个link上都构造一个flag变量y # todo 给接收端也加上一个上限 新增一个recv flag, 由f in 决定 flag = {} flag_in = {} cap_var = {} z_var = {} o_var = {} bigM = 1000 # 不需要 for k in range(K): for i in para.all_rank_set: cap_var[k,i] = model.addVar(lb = 0, ub = para.hyper_bound, vtype = gp.GRB.CONTINUOUS, name=f&#39;cap_({k},{i})&#39;) hyper_links = [hyper_link for hyper_link in para.out_direct_links[i] if hyper_link != para.pair_dict.get(i)] hyper_in_links = [hyper_link for hyper_link in para.in_direct_links[i] if hyper_link != para.pair_dict.get(i)] sio_link = para.pair_dict.get(i) # rank出口的hccs约束 flag变量表示每一个hyper link是否有发送 for hyper_link in hyper_links: flag[k, i, hyper_link] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;flag_({k},{i},{hyper_link})&#39;) # y变量。 sum_expr = gp.quicksum(f[i,hyper_link, k, c] for c in para.chunk_set) # y <= sum, sum <= y*M model.addConstr(flag[k,i, hyper_link] <= sum_expr, name=f&#39;y_l_out_({k},{i},{hyper_link})&#39;) model.addConstr(sum_expr <= 100 * flag[k, i, hyper_link], name=f&#39;y_r_out_({k},{i},{hyper_link})&#39;) # rank入口的hccs约束 flag_in表示每一个hyper link是否有接收 for hyper_in_link in hyper_in_links: flag_in[k, i, hyper_in_link] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;flag_in_({k},{i},{hyper_in_link})&#39;) sum_in_expr = gp.quicksum(f[hyper_in_link, i, k - para.delta[hyper_in_link, i], c] for c in para.chunk_set) model.addConstr(flag_in[k, i, hyper_in_link] <= sum_in_expr, name=f&#39;y_l_in_({k},{i},{hyper_in_link})&#39;) model.addConstr(sum_in_expr <= 100 * flag_in[k, i, hyper_in_link], name=f&#39;y_r_in_({k},{i},{hyper_in_link})&#39;) ### 上面其实做了两件事:1.k时刻,i给j发送的chunk数不超过;2.k时刻,到达i的chunk数不超过。 ### 新的cap计算 # cap两侧的约束 for hyper_link in hyper_links: o_var[k,i,hyper_link] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;o_({k},{i},{hyper_link})&#39;) sum_f_expr = gp.quicksum(f[i,hyper_link,k,c] for c in para.chunk_set) model.addConstr(sum_f_expr <= cap_var[k,i], name=f&#39;cap1_({i},{hyper_link},{k})&#39;) model.addConstr(cap_var[k, i] <= sum_f_expr + (1 - o_var[k,i,hyper_link])*bigM, name=f&#39;cap2_({i},{hyper_link},{k})&#39;) # o_var的约束 sum_o_expr = gp.quicksum(o_var[k,i,hyper_link] for hyper_link in hyper_links) model.addConstr(1 <= sum_o_expr, name=f&#39;o_({i},{k})&#39;) # non-linear cons: hccs 并发降速 sum_flag_expr = gp.quicksum(flag[k, i, hyper_link] for hyper_link in hyper_links) # model.addConstr(sum_flag_expr * cap_var[k,i] <= para.hyper_bound) # todo:不需要了!!!!!!!!!! # 总数上限 # 6/7 switch num 出口总数限制 model.addConstr(sum_flag_expr <= 7, name=f&#39;up_out_({i},{k})&#39;) # 入口总数限制 model.addConstr(gp.quicksum(flag_in[k,i, hyper_in_link] for hyper_in_link in hyper_in_links) <= 7, name=f&#39;up_in_({i},{k})&#39;) # 每一个link上的带宽上限取决于并发降速的带宽 todo:这是什么约束?抽象模型里没见到 for hyper_link in hyper_links: c_sum_expr = gp.quicksum(f[i, hyper_link, k, c] for c in para.chunk_set) model.addConstr(c_sum_expr <= cap_var[k,i], name=f&#39;hccsup_({i},{hyper_link},{k})&#39;) if sio_link is not None: # 针对sio的带宽进行一次限制 # sum < z * 7 7:switch num upper bound z_var[k, i] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;z_({k},{i})&#39;) model.addConstr(sum_flag_expr <= z_var[k,i] * 7, name=f&#39;z_({i},{k})&#39;) # z是一个指示变量,表示是否有hccs发送 # sum_f 小于 siocap rank_to_switch_expr = gp.quicksum(f[i, sio_link, k, c] for c in para.chunk_set) model.addConstr(rank_to_switch_expr <= para.T[i, sio_link], name=f&#39;sioup1_({i},{sio_link},{k})&#39;) # sum_f 小于 cap+1-z)M model.addConstr(rank_to_switch_expr <= cap_var[k, i] + (1-z_var[k,i])*bigM, name=f&#39;sioup2_({i},{sio_link},{k})&#39;) # sio_cap = sio_cap - z(sio_cap - cap) # rank_to_switch_expr = gp.quicksum(f[i, sio_link, k, c] for c in para.chunk_set) # model.addConstr(rank_to_switch_expr <= para.T[i, sio_link] - z_var[k,i]*(para.T[i, sio_link] - cap_var[k,i])) # objective # all2all 的 obj设计上应该考虑 将chunk数量均分,然后给每一个chunk分配好对应的位置 expr = gp.LinExpr() for rank_id, cur_chunk_lst in final_obj_dict.items(): for c in cur_chunk_lst: for k in range(K): for source_j in para.in_direct_links[rank_id]: if k + para.delta[source_j, rank_id] >= K: continue expr += (K - k + 1) * f[source_j, rank_id, k, c] model.setObjective(expr, gp.GRB.MAXIMIZE) """ # model.Params.MIPGap = 0.05 model.optimize() # model.write(&#39;tmp.lp&#39;) # env = OPTVEnv() # model_optv = OPTVModel(env) # model_optv.Read(&#39;tmp.lp&#39;) # model_optv.ComputeIIS() f_res = {key: val.X for key, val in f.items() if val.X} cap_res = {key: val.X for key, val in cap_var.items() if val.X} flag_res = {key: val.X for key, val in flag.items() if val.X} z_res = {key: val.X for key, val in z_var.items() if val.X} """ # return f_res, model.Runtime return model def model_construct_hyper_edge_no_sio_switch(para:Para_910C, K = None, initial_rank_dict = None, final_obj_dict = None): # 对于rolling horizon场景下适用 # initial rank chunk && final obj dict if not initial_rank_dict: initial_rank_dict = copy.deepcopy(para.rank_chunk_dict) if not final_obj_dict: final_obj_dict = copy.deepcopy(para.final_obj_dict) if not K: K = para.K f = {} b = {} model = gp.Model() print(&#39;var f and var b generating for each rank&#39;) for k in range(K): for i in para.all_rank_set: # para.all_rank_set={0,3,16,19} for c in para.chunk_set: # {0, 1, 2, 3, 12, 13, 14, 15, 64, 65, 66, 67, 76, 77, 78, 79} b[i, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;b_({i},{k},{c})&#39;) # 选择他能发出去的link for j in para.out_direct_links.get(i, []): # para.out_direct_links={0: [16, 3, 19], 3: [0, 16, 19], 16: [0, 3, 19], 19: [0, 16, 3]}。也就是rank需要发送到的rank f[i, j, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;f_({i},{j},{k},{c})&#39;) print(&#39;hyper edge band upper bound&#39;) # 对于每一个rank来说,其给switch发出的大小受总量的限制 for k in range(K): for i in para.all_rank_set: # hccs 发送固有上限 hyper_links = [hyper_link for hyper_link in para.out_direct_links[i] if hyper_link != para.pair_dict.get(i)] ##!!pair_dict={},是空的,或许是遗留问题?out_direct_links里本来就没有直连的 model.addConstr(gp.quicksum(f[i, j, k, c] for c in para.chunk_set for j in hyper_links) <= para.hyper_bound) # hccscap约束 # sio 发送固有上限 # 可能存在没有sio链路的情况 ## 确实是,目前没有sio链路,都是hyperedge发送。 nei_rank = para.pair_dict.get(i) if nei_rank: model.addConstr(gp.quicksum(f[i, nei_rank, k, c] for c in para.chunk_set) <= para.T[i, nei_rank]) # 目前T里面都是0 print(&#39;status constrain flow constrain&#39;) # 再进行状态约束,这里只进行节点往外发送的状态约束 for k in range(K): for i in para.all_rank_set: init_chunk_set = initial_rank_dict[i] for c in para.chunk_set: # 确定前一个时刻的状态从哪里取值 tmp_var = None if k == 0 and c in init_chunk_set: # 在0时刻上,且c属于原始chunk的,给f的右端项一个1,其余的直接从b中get值出来 tmp_var = 1 # 代表初始的b,也就是可以发送 else: tmp_var = b.get((i, k - 1, c), 0) # 这里的get 0,有点危险。k-1时刻有,k时刻才能发送 # 发送约束 for target_j in para.out_direct_links[i]: model.addConstr(f[i, target_j, k, c] <= tmp_var, name=f&#39;status_({k},{i},{target_j},{c})&#39;) # 状态平衡 ### 我有一个疑问,这里没有过滤c,也就是说source的rank可能给i发送他自己的c?虽然是这样子不会出错。 in_expr = gp.quicksum( f.get((source_j, i, k - para.delta[source_j, i], c), 0) for source_j in para.in_direct_links[i]) model.addConstr(b[i, k, c] == tmp_var + in_expr, name=f&#39;flow_({k},{i},{c})&#39;) # 给每一个link上都构造一个flag变量y # todo 给接收端也加上一个上限 新增一个recv flag, 由f in 决定 flag = {} flag_in = {} cap_var = {} z_var = {} for k in range(K): for i in para.all_rank_set: cap_var[k,i] = model.addVar(lb = 0, ub = para.hyper_bound, vtype = gp.GRB.CONTINUOUS) hyper_links = [hyper_link for hyper_link in para.out_direct_links[i] if hyper_link != para.pair_dict.get(i)] hyper_in_links = [hyper_link for hyper_link in para.in_direct_links[i] if hyper_link != para.pair_dict.get(i)] sio_link = para.pair_dict.get(i) # rank出口的hccs约束 flag变量表示每一个hyper link是否有发送 for hyper_link in hyper_links: flag[k, i, hyper_link] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;flag_({k},{i},{hyper_link})&#39;) # y变量。 sum_expr = gp.quicksum(f[i,hyper_link, k, c] for c in para.chunk_set) # y <= sum, sum <= y*M model.addConstr(flag[k,i, hyper_link] <= sum_expr) model.addConstr(sum_expr <= 100 * flag[k, i, hyper_link]) # rank入口的hccs约束 flag_in表示每一个hyper link是否有接收 for hyper_in_link in hyper_in_links: flag_in[k, i, hyper_in_link] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;flag_in_({k},{i},{hyper_in_link})&#39;) sum_in_expr = gp.quicksum(f[hyper_in_link, i, k - para.delta[hyper_in_link, i], c] for c in para.chunk_set) model.addConstr(flag_in[k, i, hyper_in_link] <= sum_in_expr) model.addConstr(sum_in_expr <= 100 * flag_in[k, i, hyper_in_link]) # non-linear cons: hccs 并发降速 sum_flag_expr = gp.quicksum(flag[k, i, hyper_link] for hyper_link in hyper_links) model.addConstr(sum_flag_expr * cap_var[k,i] <= para.hyper_bound) # 总数上限 # 6/7 switch num 出口总数限制 model.addConstr(sum_flag_expr <= 7) # 入口总数限制 model.addConstr(gp.quicksum(flag_in[k,i, hyper_in_link] for hyper_in_link in hyper_in_links) <= 7) # 每一个link上的带宽上限取决于并发降速的带宽 for hyper_link in hyper_links: c_sum_expr = gp.quicksum(f[i, hyper_link, k, c] for c in para.chunk_set) model.addConstr(c_sum_expr <= cap_var[k,i]) if sio_link is not None: # 针对sio的带宽进行一次限制 # sum < z * 7 7:switch num upper bound z_var[k, i] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;z_({k},{i})&#39;) model.addConstr(sum_flag_expr <= z_var[k,i] * 7) # z是一个指示变量,表示是否有hccs发送 # sio_cap = sio_cap - z(sio_cap - cap) rank_to_switch_expr = gp.quicksum(f[i, sio_link, k, c] for c in para.chunk_set) model.addConstr(rank_to_switch_expr <= para.T[i, sio_link] - z_var[k,i]*(para.T[i, sio_link] - cap_var[k,i])) # objective # all2all 的 obj设计上应该考虑 将chunk数量均分,然后给每一个chunk分配好对应的位置 expr = gp.LinExpr() for rank_id, cur_chunk_lst in final_obj_dict.items(): for c in cur_chunk_lst: for k in range(K): for source_j in para.in_direct_links[rank_id]: if k + para.delta[source_j, rank_id] >= K: continue expr += (K - k + 1) * f[source_j, rank_id, k, c] model.setObjective(expr, gp.GRB.MAXIMIZE) """ # model.Params.MIPGap = 0.05 model.optimize() # model.write(&#39;tmp.lp&#39;) # env = OPTVEnv() # model_optv = OPTVModel(env) # model_optv.Read(&#39;tmp.lp&#39;) # model_optv.ComputeIIS() f_res = {key: val.X for key, val in f.items() if val.X} cap_res = {key: val.X for key, val in cap_var.items() if val.X} flag_res = {key: val.X for key, val in flag.items() if val.X} z_res = {key: val.X for key, val in z_var.items() if val.X} return f_res, model.Runtime """ return model def model_construct_hyper_edge(para:Para_910C): f = {} b = {} model = gp.Model() print(&#39;var f and var b generating for each rank&#39;) for k in range(para.K): for i in range(para.N): for c in para.chunk_set: b[i, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;b_[{i},{k},{c}]&#39;) # 选择他能发出去的link for j in para.out_direct_links.get(i, []): f[i, j, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;f_[{i},{j},{k},{c}]&#39;) print(&#39;virtual switch node gene&#39;) for k in range(para.K): for c in para.chunk_set: for switch_idx in para.all_switch_set: for out_rank in para.switch_out_links.get(switch_idx, []): f[switch_idx, out_rank, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;f_[{switch_idx},{out_rank},{k},{c}]&#39;) print(&#39;rdma switch&#39;) if para.rdma_switch: # 如果此时是多机环境,则需要对其rdma switch进行建模 # 此时先实例化送rdma switch的发送 for k in range(para.K): for c in para.chunk_set: for rank_id in range(para.N): f[para.rdma_switch, rank_id, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;f_[{para.rdma_switch},{rank_id},{k},{c}]&#39;) # print(&#39;cap constrain&#39;) # 先进行rank 之间的cap约束 # 分为三个部分:rank, sio-switch, rdma-switch 的out link # for k in range(para.K): # # rank # # 考虑总约束的情况下,这里可以注释 # for i in range(para.N): # for j in para.out_direct_links[i]: # # cap = para.T[i, j] # model.addConstr(gp.quicksum(f[i, j, k, c] for c in para.chunk_set) <= cap, name=f&#39;cap_[{k},{i},{j}]&#39;) # # # sio switch # # 这里要看一下 可以进行注释 # for switch_idx in para.all_switch_set: # out_links = para.switch_out_links[switch_idx] # for j in range(out_links): # cap = para.T[switch_idx, j] # model.addConstr(gp.quicksum(f[switch_idx, j, k, c] for c in para.chunk_set) <= cap, name=f&#39;cap_[{k},{switch_idx},{j}]&#39;) # # # rdma switch # if para.rdma_switch: # out_links = list(range(para.N)) # for j in out_links: # cap = para.T[para.rdma_switch, j] # model.addConstr(gp.quicksum(f[para.rdma_switch, j, k, c] for c in para.chunk_set) <= cap, # name=f&#39;cap_[{k},{para.rdma_switch},{j}]&#39;) print(&#39;hyper edge band upper bound&#39;) # 对于每一个rank开说,其给switch发出的大小受总量的限制 for k in range(para.K): for i in range(para.N): # 此处使用 int类型表示其是否是一个 hyper link hyper_links = [hyper_link for hyper_link in para.out_direct_links[i] if type(hyper_link) == int] model.addConstr(gp.quicksum(f[i, j, k, c] for c in para.chunk_set for j in hyper_links) <= para.hyper_bound) print(&#39;status constrain flow constrain&#39;) # 再进行状态约束,这里只进行节点往外发送的状态约束 for k in range(para.K): for i in range(para.N): init_chunk_set = para.rank_chunk_dict[i] for c in para.chunk_set: # 确定前一个时刻的状态从哪里取值 tmp_var = None if k == 0 and c in init_chunk_set: # 在0时刻上,且c属于原始chunk的,给f的右端项一个1,其余的直接从b中get值出来 tmp_var = 1 else: tmp_var = b.get((i, k - 1, c), 0) # 发送约束 for target_j in para.out_direct_links[i]: model.addConstr(f[i, target_j, k, c] <= tmp_var, name=f&#39;status_[{k},{i},{target_j},{c}]&#39;) # 状态平衡 in_expr = gp.quicksum( f.get((source_j, i, k - para.delta[source_j, i], c), 0) for source_j in para.in_direct_links[i]) model.addConstr(b[i, k, c] == tmp_var + in_expr, name=f&#39;flow_[{k},{i},{c}]&#39;) print(&#39;switch balance&#39;) # 在每一个switch上保持进出一致:k,c # 主要是给sio 虚拟switch的约束 for switch_idx in para.all_switch_set: out_links = para.switch_out_links[switch_idx] in_links = para.switch_in_links[switch_idx] for k in range(para.K): for c in para.chunk_set: in_expr = gp.quicksum(f[i, switch_idx, k, c] for i in in_links) out_expr = gp.quicksum(f[switch_idx, j, k, c] for j in out_links) model.addConstr(in_expr == out_expr, name = &#39;balance&#39;) # # todo non-linear # # 给每一个link上都构造一个flag变量y # flag = {} # cap_var = {} # z_var = {} # for k in range(para.K): # for i in range(para.N): # cap_var[k,i] = model.addVar(lb = 0, ub = para.hyper_bound, vtype = gp.GRB.CONTINUOUS) # z_var[k,i] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;z_[{k},{i}]&#39;) # hyper_links = [hyper_link for hyper_link in para.out_direct_links[i] if type(hyper_link) == int] # sio_link = [hyper_link for hyper_link in para.out_direct_links[i] if type(hyper_link) != int][0] # # for hyper_link in hyper_links: # flag[k, i, hyper_link] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;flag_[{k},{i},{hyper_link}]&#39;) # # sum_expr = gp.quicksum(f[i,hyper_link, k, c] for c in para.chunk_set) # # # y <= sum, sum <= y*M # model.addConstr(flag[k,i, hyper_link] <= sum_expr) # model.addConstr(sum_expr <= 100 * flag[k, i, hyper_link]) # # sum_flag_expr = gp.quicksum(flag[k,i, hyper_link] for hyper_link in hyper_links) # # non-linear cons # model.addConstr(sum_flag_expr * cap_var[k,i] <= para.hyper_bound) # # # 总数上限 # model.addConstr(sum_flag_expr <= 7) # # # 每一个link上的带宽 # for hyper_link in hyper_links: # c_sum_expr = gp.quicksum(f[i, hyper_link, k, c] for c in para.chunk_set) # model.addConstr(c_sum_expr <= cap_var[k,i]) # # # 针对sio的带宽进行一次限制 # switch_idx = para.rank_to_switch[rank_id] # model.addConstr(sum_flag_expr <= z_var[k,i] * (para.N + 1) ) # # rank_to_switch_expr = gp.quicksum(f[rank_id, switch_idx, k, c] for c in para.chunk_set) # model.addConstr(rank_to_switch_expr <= para.T[rank_id, switch_idx] - z_var[k,i](para.T[rank_id, switch_idx] - cap_var[k,i])) # 将rnk0 - rank1的情况ban掉 for switch_idx in para.all_switch_set: rank0, rank1 = para.switch_out_links[switch_idx] for k in range(para.K): for c in para.chunk_set: model.addConstr(f[rank0, rank1, k, c] == 0) model.addConstr(f[rank1, rank0, k, c] == 0) # objective # all2all 的 obj设计上应该考虑 将chunk数量均分,然后给每一个chunk分配好对应的位置 expr = gp.LinExpr() for rank_id, cur_chunk_lst in para.final_obj_dict.items(): for c in cur_chunk_lst: for k in range(para.K): for source_j in para.in_direct_links[rank_id]: if k + para.delta[source_j, rank_id] >= para.K: continue expr += (para.K - k + 1) * f[source_j, rank_id, k, c] model.setObjective(expr, gp.GRB.MAXIMIZE) model.optimize() # model.write(&#39;tmp.lp&#39;) # env = OPTVEnv() # model_optv = OPTVModel(env) # model_optv.Read(&#39;tmp.lp&#39;) # model_optv.ComputeIIS() f_res = {key: val.X for key, val in f.items() if val.X} # cap_res = {key: val.X for key, val in cap.items() if val.X} # flag_res = {key: val.X for key, val in flag.items() if val.X} with open(&#39;ttt.pkl&#39;, &#39;wb&#39;) as f: pickle.dump((f_res, cap_res, flag_res), f) exit() return f_res, model.Runtime def model_construct(para:Para_910C): # 使用switch balance的模型 f = {} b = {} model = gp.Model() print(&#39;var f and var b generating for each rank&#39;) for k in range(para.K): for i in range(para.N): for c in para.chunk_set: b[i, k, c] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=f&#39;b_[{i},{k},{c}]&#39;) # 选择他能发出去的link for j in para.out_direct_links.get(i, []): f[i,j,k,c] = model.addVar(lb = 0, ub = 1, vtype = gp.GRB.BINARY, name = f&#39;f_[{i},{j},{k},{c}]&#39;) # 需要把 switch往外发送的f的实例化 for k in range(para.K): for server_idx in range(para.server_num): rank_lst = para.server_rank[server_idx] for switch_idx in para.server_switch[server_idx]: for n in rank_lst: for c in para.chunk_set: f[switch_idx, n, k, c] = model.addVar(lb = 0, ub = 1, vtype = gp.GRB.BINARY, name =f&#39;f_[{switch_idx},{n},{k},{c}]&#39;) print(&#39;cap constrain&#39;) # 先进行rank 之间的cap约束 for k in range(para.K): for i in range(para.N): for j in para.out_direct_links[i]: cap = para.T[i,j] model.addConstr(gp.quicksum(f[i,j,k,c] for c in para.chunk_set) <= cap, name = f&#39;cap_[{k},{i},{j}]&#39;) # 再进行switch出来的cap约束 for k in range(para.K): for server_idx in range(para.server_num): rank_lst = para.server_rank[server_idx] for switch_idx in para.server_switch[server_idx]: for rank_id in rank_lst: cap = para.T[switch_idx, rank_id] model.addConstr(gp.quicksum(f[switch_idx, rank_id, k, c] for c in para.chunk_set) <= cap, name=f&#39;cap_[{k},{switch_idx},{rank_id}]&#39;) print(&#39;status constrain flow constrain&#39;) # 再进行状态约束,这里只进行节点往外发送的状态约束 for k in range(para.K): for i in range(para.N): init_chunk_set = para.rank_chunk_dict[i] for c in para.chunk_set: # 确定前一个时刻的状态从哪里取值 tmp_var = None if k == 0 and c in init_chunk_set: # 在0时刻上,且c属于原始chunk的,给f的右端项一个1,其余的直接从b中get值出来 tmp_var = 1 else: tmp_var = b.get((i, k - 1, c), 0) # 发送约束 for target_j in para.out_direct_links[i]: model.addConstr(f[i,target_j,k,c] <= tmp_var, name = f&#39;status_[{k},{i},{target_j},{c}]&#39;) # 状态平衡 in_expr = gp.quicksum(f.get((source_j, i, k-para.delta[source_j, i], c), 0) for source_j in para.in_direct_links[i]) model.addConstr(b[i, k, c] == tmp_var + in_expr, name = f&#39;flow_[{k},{i},{c}]&#39;) # switch cons for server_idx in range(para.server_num): switch_lst = para.server_switch[server_idx] rank_lst = para.server_rank[server_idx] for switch_idx in switch_lst: for rank_id in rank_lst: other_rank_lst = [other_rank for other_rank in rank_lst if other_rank != rank_id] for k in range(para.K): for c in para.chunk_set: recv_expr = gp.quicksum([f.get((switch_idx, other_rank, k, c), 0) for other_rank in other_rank_lst]) model.addConstr(f[rank_id, switch_idx, k, c] == recv_expr, name = f&#39;switch_[{switch_idx},{rank_id},{k},{c}]&#39;) # model.addConstr(f[0, &#39;s0&#39;, 0, 0] == 1, name = &#39;tmp&#39;) # model.addConstr(f[&#39;s0&#39;,3, 0, 0] == 1, name=&#39;tmp&#39;) # todo 以上述为例,s0 给 3发出东西后,按道理 其也应该等于 除了1 之外的其他rank给s0发送一个chunk 0才能维持所有的switch约束保持正确 因此有bug # objective expr = gp.LinExpr() for rank_id in range(para.N): initial_chunk = para.rank_chunk_dict[rank_id] for k in range(para.K): for j in range(para.N): if rank_id == j: continue for source_j in para.in_direct_links[j]: if k + para.delta[source_j, j] >= para.K: continue for c in initial_chunk: expr += (para.K - k + 1) * f[source_j, j, k, c] model.setObjective(expr, gp.GRB.MAXIMIZE) model.optimize() # model.write(&#39;tmp.lp&#39;) # env = OPTVEnv() # model_optv = OPTVModel(env) # model_optv.Read(&#39;tmp.lp&#39;) # model_optv.ComputeIIS() f_res = {key:val.X for key, val in f.items() if val.X} return f_res def rolling_horizon_solve(para:Para_910C, iter_K): # todo 形式上只要将原模型中的已经到达最终位置的 计入到每一个rank的initial chunk就可以了? # todo 要不要修改一下obj,可以修改一下! final_obj_dict = copy.deepcopy(para.final_obj_dict) initial_rank_chunks = copy.deepcopy(para.rank_chunk_dict) rolling_iter = 0 f_final_res = {} final_solving_time = 0 while final_obj_dict: f_res, solving_time = model_construct_hyper_edge_no_sio_switch(para, iter_K, initial_rank_chunks, final_obj_dict) f_res_filter = utils.res_filter(f_res, initial_rank_chunks, final_obj_dict) initial_rank_chunks, final_obj_dict = utils.iter_update(f_res_filter, initial_rank_chunks, final_obj_dict) for i,j,k,c in f_res.keys(): f_final_res[i,j,k + rolling_iter * iter_K, c] = 1 final_solving_time += solving_time rolling_iter += 1 print(f_final_res) return f_final_res, final_solving_time def script(): with open(&#39;tmp.pkl&#39;, &#39;rb&#39;) as f: data = pickle.load(f) res = data[&#39;res&#39;] para = data[&#39;para&#39;] new_res = utils.switch_prepare(res) new_res = utils.res_filter(new_res, para.rank_chunk_dict, para.final_obj_dict) if para.task_mode == TASK_Type.REDUCE_SCATTER: # 针对reduce scatter进行进一步的 结果改进! # todo 这里有点问题,从0 new_res = utils.reduce_scatter_inverse(new_res, para.rank_chunk_dict, para.final_obj_dict) tmp = list(new_res.keys()) tmp.sort(key = lambda x:x) for i,j,k,c in tmp: # 算一下c是哪一个chunk上来的 source = None for rank in range(para.N): chunks = para.rank_chunk_dict.get(rank, []) if c in chunks: source = rank break print(i,j,k,c,"----", source) if source != i: print(&#39;attention&#39;) ki_dict = {} for (i, j, k, c) in tmp: ki_dict[k, i] = ki_dict.get((k, i), []) + [(i, j, k, c)] for (k, i), f_lst in ki_dict.items(): nei_rank = para.pair_dict.get(i) send_dict = {} for _, j, _, c in f_lst: send_dict[j] = send_dict.get(j, []) + [c] hyper_link_num = sum([1 for j in send_dict.keys() if j != nei_rank]) sio_send_num = len(send_dict.get(nei_rank, [])) max_sum = max([len(vals) for vals in send_dict.values()]) if hyper_link_num: split_bound = para.hyper_bound / hyper_link_num for hyper_rank, chunks in send_dict.items(): assert split_bound >= len(chunks) assert sio_send_num <= split_bound else: split_bound = para.T[0,1] assert sio_send_num <= para.T[0,1] print(f&#39;epoch:{k}, rank:{i}, hyper_send:{hyper_link_num}, split_bound:{split_bound}, sio_send_num:{sio_send_num}, max_sum:{max_sum}&#39;) obj = max([key[2] for key in new_res.keys()]) + 1 return new_res, obj def find_factors(n): return [i for i in range(1, n + 1) if n % i == 0] def process_instance(i): # server_num采样 server_num = random.randint(1, 4) if server_num != 1: P = random.choice(all_lists) elif server_num == 1: P = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] else: raise Exception("unknown server num") # task_mode = TASK_Type.ALL_GATHER task_mode = random.choice([TASK_Type.REDUCE_SCATTER, TASK_Type.ALL_GATHER]) # Input_size采样 num = random.choice([32,64,128]) Input_size = num * 1024 * 1024 # Bytes # C采样 if task_mode == TASK_Type.REDUCE_SCATTER or task_mode == TASK_Type.ALL2ALL: C = random.choice([len(P) * server_num,len(P) * server_num *2,len(P) * server_num*3]) elif task_mode == TASK_Type.ALL_GATHER: C = random.randint(1, 10) else: raise Exception("unknown task mode") # HCCSCAP采样 HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) while HCCSCAP <= 0: print("chunk太大了,一个epoch发不完") C += 1 HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) ### HCCSCAP不做上限约束 # while HCCSCAP > 10: # print("chunk太小,HCCSCAP太大") # Input_size *= 2 # HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) print("HCCSCAP:", HCCSCAP) # para与instance生成 para = Para_910C() file_name = f&#39;{topo_name}_{server_num}_{P}_{HCCSCAP}_{Input_size / 1024 / 1024 / 1024}_{C}_{model_mode.value}_{task_mode.value}&#39; if model_mode == Model_Type.HYPER_EDGE: para.init_hyper_edge_no_sio(topo_name, tau, C, K, Input_size, server_num, switch, task_mode, P) m = model_construct_hyper_edge_no_sio_switch_linear(para) else: raise Exception("wrong model_mode") ## 保存 m.write(filename=f"./instance/train/910C/instance_{i}_{file_name}.lp") # with open(f&#39;./para/train/910c_para_{i}_{file_name}.pkl&#39;, &#39;wb&#39;) as f: # pickle.dump({&#39;file_name&#39;: file_name, # &#39;para&#39;: para}, f) def process_test_instance(i): # server_num采样 server_num = random.randint(1, 4) if server_num != 1: P = random.choice(all_lists) elif server_num == 1: P = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] else: raise Exception("unknown server num") # task_mode = TASK_Type.ALL_GATHER task_mode = TASK_Type.REDUCE_SCATTER # task_mode = TASK_Type.ALL2ALL # task_mode = random.choice([TASK_Type.REDUCE_SCATTER, TASK_Type.ALL_GATHER, TASK_Type.ALL2ALL]) # task_mode = random.choice([TASK_Type.REDUCE_SCATTER, TASK_Type.ALL_GATHER, TASK_Type.ALL2ALL]) # Input_size采样 num = random.choice([32, 64, 128]) Input_size = num * 1024 * 1024 # Bytes # C采样 if task_mode == TASK_Type.REDUCE_SCATTER or task_mode == TASK_Type.ALL2ALL: C = random.choice([len(P) * server_num, len(P) * server_num * 2]) elif task_mode == TASK_Type.ALL_GATHER: C = random.randint(1, 10) else: raise Exception("unknown task mode") # HCCSCAP采样 HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) while HCCSCAP <= 0: print("chunk太大了,一个epoch发不完") C += 1 HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) ### HCCSCAP不做上限约束 # while HCCSCAP > 10: # print("chunk太小,HCCSCAP太大") # Input_size *= 2 # HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) print("HCCSCAP:", HCCSCAP) # para与instance生成 para = Para_910C() file_name = f&#39;{topo_name}_{server_num}_{P}_{HCCSCAP}_{Input_size / 1024 / 1024 / 1024}_{C}_{model_mode.value}_{task_mode.value}&#39; if model_mode == Model_Type.HYPER_EDGE: para.init_hyper_edge_no_sio(topo_name, tau, C, K, Input_size, server_num, switch, task_mode, P) m = model_construct_hyper_edge_no_sio_switch_linear(para) # m = model_construct_hyper_edge_no_sio_switch(para) else: raise Exception("wrong model_mode") ## 保存 m.write(filename=f"./instance/test_zhongqi_reduce_scatter/910C/instance_{i}_{file_name}.lp") # m.write(filename=f"./instance/test/910C/instance_{i}_{file_name}.lp") # with open(f&#39;./para/test/910c_para_{i}_{file_name}.pkl&#39;, &#39;wb&#39;) as f: # pickle.dump({&#39;file_name&#39;: file_name, # &#39;para&#39;: para}, f) if __name__ == &#39;__main__&#39;: instanceNum = 1 tau = 4.4e2 / 1e6 # Epoch duration 1e-3 for 910A K = 8 # An upperbound of the total epoch number switch = 7 # 不动 model_mode = Model_Type.HYPER_EDGE topo_name = "910C_2" # switch_node random.seed(1145) para_folder = "./para" if not os.path.exists(para_folder): os.makedirs(para_folder, exist_ok=True) if not os.path.exists(os.path.join(para_folder, "train")): os.makedirs(os.path.join(para_folder, "train"), exist_ok=True) if not os.path.exists(os.path.join(para_folder, "test")): os.makedirs(os.path.join(para_folder, "test"), exist_ok=True) instance_folder = "./instance" if not os.path.exists(instance_folder): os.makedirs(instance_folder, exist_ok=True) if not os.path.exists(os.path.join(instance_folder, "train/910C")): os.makedirs(os.path.join(instance_folder, "train/910C"), exist_ok=True) if not os.path.exists(os.path.join(instance_folder, "test/910C")): os.makedirs(os.path.join(instance_folder, "test/910C"), exist_ok=True) ### P采样 all_lists = [] for r in range(1, 17): for p in itertools.combinations(range(16), r): if len(list(p)) <= 3: all_lists.append(list(p)) ## 实例数目 # todo:目前的test还是all gather trainNum = 1 # trainNum = 200 testNum = 10 processes = [] max_processes = 32 # 定义最多同时开启的进程数量为10 if False: with multiprocessing.Pool(processes=max_processes) as pool: pool.map(process_instance, range(trainNum)) if True: with multiprocessing.Pool(processes=max_processes) as pool: pool.map(process_test_instance, range(testNum)) # for i in range(trainNum): # p = multiprocessing.Process(target=process_instance, args=(i,)) # p.start() # processes.append(p) """ for i in range(trainNum): server_num = random.randint(1, 4) P = random.choice(all_lists) # task_mode = TASK_Type.REDUCE_SCATTER task_mode = TASK_Type.ALL_GATHER # C采样 if task_mode == TASK_Type.REDUCE_SCATTER or task_mode == TASK_Type.ALL2ALL: C = len(P) * server_num # todo:待修改 # factors = find_factors(len(P)* server_num) # 找到len(P)* server_num的所有因子(todo:我感觉需要排除1,最大因子,以及一些过大的C) # C = random.choice(factors) elif task_mode == TASK_Type.ALL_GATHER: C = random.randint(1,20) # todo:待修改 else: raise Exception("unknown task mode") # HCCSCAP采样 num = random.randint(16, 512) Input_size = num * 1024 * 1024 # Bytes HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) while HCCSCAP <= 0: print("chunk太大了,一个epoch发不完") Input_size /= 2 HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) while HCCSCAP >= 10: print("chunk太小,HCCSCAP太大") Input_size *= 2 HCCSCAP = round(140 * 1024 ** 3 / (Input_size / C) * tau) print("HCCSCAP:", HCCSCAP) # para与instance生成 para = Para_910C() file_name = f&#39;{topo_name}_{server_num}_{P}_{Input_size / 1024 / 1024 / 1024}_{C}_{model_mode.value}_{task_mode.value}&#39; if model_mode == Model_Type.HYPER_EDGE: para.init_hyper_edge_no_sio(topo_name, tau, C, K, Input_size, server_num, switch, task_mode, P) m = model_construct_hyper_edge_no_sio_switch_linear(para) else: # model_mode == &#39;switch&#39; raise Exception("这里没改好,而且似乎也不需要") para.init(topo_name, tau, C, K, Input_size, server_num, 6, task_mode) f_res, solving_time = model_construct_switch_node(para) ## 保存 m.write(filename=f"./instance/train/910c_instance_{i}_{file_name}.lp") with open(f&#39;./para/train/910c_para_{i}_{file_name}.pkl&#39;, &#39;wb&#39;) as f: pickle.dump({&#39;file_name&#39;: file_name, &#39;para&#39;: para}, f) """
最新发布
11-11
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(&#39;scenarios/scenarios_S5.csv&#39;) scenarios = scenarios_df[&#39;Scenario&#39;].tolist() mu_s = {row[&#39;Scenario&#39;]: row[&#39;mu&#39;] for _, row in scenarios_df.iterrows()} prob_hist = {row[&#39;Scenario&#39;]: row[&#39;Probability&#39;] 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 = [&#39;C1&#39;, &#39;C2&#39;] # 二类车主 D = [&#39;D1&#39;, &#39;D2&#39;] # 一类车主 J = [&#39;J1&#39;, &#39;J2&#39;] # 维修中心 L = [&#39;L1&#39;, &#39;L2&#39;] # 回收价格选项 H = [&#39;H1&#39;, &#39;H2&#39;] # 4S店 F = [&#39;F1&#39;, &#39;F2&#39;] # 运营公司 K = [&#39;K1&#39;, &#39;K2&#39;] # 综合利用企业 G = [&#39;G1&#39;, &#39;G2&#39;] # 拆车厂 M = [&#39;M1&#39;, &#39;M2&#39;] # 梯次利用中心 N = [&#39;N1&#39;, &#39;N2&#39;] # 再生利用企业 A = [&#39;Co&#39;, &#39;Mn&#39;, &#39;Ni&#39;, &#39;Li&#39;, &#39;Others&#39;] # 金属材料 # 参数赋值 # 废旧动力电池单位剩余容量 L_bar = 0.7 # 回收量参数(吨) Q1 = {&#39;L1&#39;: 112, &#39;L2&#39;: 75} # 一类车主回收量 Q2 = {&#39;L1&#39;: 112, &#39;L2&#39;: 75} # 二类车主回收量 Q3 = {&#39;L1&#39;: 75, &#39;L2&#39;: 50} # 拆车厂回收量 # 经济参数(元) pr = {&#39;L1&#39;: 4000, &#39;L2&#39;: 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 = {&#39;Co&#39;: 350000, &#39;Mn&#39;: 6600, &#39;Ni&#39;: 4100, &#39;Li&#39;: 449000, &#39;Others&#39;: 25000} # 市场价格(元/吨) alpha_a = {&#39;Co&#39;: 0.0422 * 0.98, &#39;Mn&#39;: 0.0109 * 0.98, &#39;Ni&#39;: 0.3427 * 0.98, &#39;Li&#39;: 0.0419 * 0.85, &#39;Others&#39;: 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}" ) # KM之间的流量 (依赖场景的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}" ) # KM之间的流量 (依赖场景的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 = { &#39;X_j&#39;: X_j_val, &#39;X_k&#39;: X_k_val, &#39;X_m&#39;: X_m_val, &#39;X_n&#39;: X_n_val, &#39;Z_l&#39;: Z_l_val, &#39;obj&#39;: master_obj, &#39;sub_obj&#39;: master_obj + stage1_cost_val, &#39;active_scenarios&#39;: 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 = { &#39;X_j&#39;: X_j_val, &#39;X_k&#39;: X_k_val, &#39;X_m&#39;: X_m_val, &#39;X_n&#39;: X_n_val, &#39;Z_l&#39;: Z_l_val, &#39;obj&#39;: total_obj_estimate, &#39;sub_obj&#39;: sub_obj, &#39;active_scenarios&#39;: 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 = { &#39;X_j&#39;: X_j_val, &#39;X_k&#39;: X_k_val, &#39;X_m&#39;: X_m_val, &#39;X_n&#39;: X_n_val, &#39;Z_l&#39;: Z_l_val, &#39;obj&#39;: master_obj, &#39;sub_obj&#39;: master_obj + stage1_cost_val, &#39;active_scenarios&#39;: active_scenarios.copy() } print(f"找到可行解,目标值: {best_sol[&#39;obj&#39;]:,.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[&#39;active_scenarios&#39;]) master.optimize() if master.status == GRB.OPTIMAL: # 提取概率分布 p_s_val = {s: p_s[s].X for s in best_solution[&#39;active_scenarios&#39;]} # 提取关键物流变量(以第一个激活场景为例) s0 = best_solution[&#39;active_scenarios&#39;][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[&#39;X_j&#39;][j] < 1 else ("合营" if best_solution[&#39;X_j&#39;][j] == 1 else "新建") print(f"{j}: {mode}") # 价格选择 for l in L: if best_solution[&#39;Z_l&#39;][l] > 0.5: print(f"\n选择回收价格: {pr[l]}元 (选项{l})") # 概率分布 print("\n最坏情况概率分布:") for s in best_solution[&#39;active_scenarios&#39;]: 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[&#39;X_j&#39;][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[&#39;Z_l&#39;][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[&#39;sub_obj&#39;]:,.2f}元") print(f"总系统利润: {best_solution[&#39;obj&#39;]:,.2f}元") # 矩约束验证 exp_mu = sum(p_s_val.get(s, 0) * mu_s[s] for s in best_solution[&#39;active_scenarios&#39;]) exp_sigma2 = sum(p_s_val.get(s, 0) * (mu_s[s] ** 2) for s in best_solution[&#39;active_scenarios&#39;]) 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、付费专栏及课程。

余额充值