Nested Summation constraint for MILP problem
回答済み
Can anyone help me how to write this constraint? As I am new to python and started using gurobi for optimization.
-
Please have a look at the post Nested sum using quicksum().
0 -
I have tried a lot , I don't know I am unable to do it. I need to use continue and conditional statement as well in my constraint.
0 -
Could you provide a minimal working example of what exactly is not working?
0 -
Resourcebalance=steel.addConstrs((Resources[res,t]==Resources[res,time[time.index(t)-1]] +
gp.quicksum(rtn_profile[res][task][theta]* Tasks[task,t-theta]
for task in task_list1
if task not in rtn_profile[res] continue duration = len(rtn_profile[res][task]) – 1
for theta in range(0, min(duration, t) + 1))
for t in time
for res in res_list1),"Resourcebalance")This is what I am trying to do, and my logic is this.
Resources,Tasks=Variables
task_list1,res_list1,time=list
rtn_profile=dictionary
0 -
Your code snippet is not reproducible, i.e., one cannot copy paste and run it. Please provide at least a minimal data set for all variables, lists, and dictionaries you use.
In your quicksum \(\texttt{for}\)-loop, you use \(\texttt{continue}\) and then set duration to some value. This is not possible. Could you please elaborate what exactly this \(\texttt{if}\)-clause with \(\texttt{continue}\) and setting of the \(\texttt{duration}\) value shall do?
0 -
import numpy as np
import pandas as pd
import xlrd
import csv
from collections import OrderedDict
import mathimport gurobipy as gp
from gurobipy import GRBgroup2heats = {'1': [1, 2, 3, 4], '2': [5, 6, 7, 8], '3': [9, 10, 11, 12],
'4': [13, 14, 15, 16, 17], '5': [18, 19, 20], '6': [21, 22, 23, 24]}equip2mw = {'EAF': 85, 'AOD': 2, 'LF': 2, 'CC1': 7, 'CC2': 7}
energy_price = [23.0000, 23.1400, 24.6900, 24.450, 24.4600, 24.290,
24.8700, 30.1600, 33.3900, 34.5200, 35.4900, 32.5400,
32.3100, 29.4300, 27.2000, 26.2900, 25.9000, 27.1900,
34.9300, 45.2000, 35.4700, 33.4800, 32.1000, 27.9300]equip2process_time = {
'EAF1':{'1':80.0,'2':80.0,'3':80.0,'4':80.0,'5':85.0,'6':85.0,'7':85.0,'8':85.0,'9':90.0,'10':90.0,'11':90.0,'12':90.0,'13':85.0,'14':85.0,'15':85.0,'16':85.0,'17':80.0,'18':80.0,'19':80.0,'20':80.0,'21':80.0,'22':80.0,'23':80.0,'24':80.0},
'EAF2':{'1':80.0,'2':80.0,'3':80.0,'4':80.0,'5':85.0,'6':85.0,'7':85.0,'8':85.0,'9':90.0,'10':90.0,'11':90.0,'12':90.0,'13':85.0,'14':85.0,'15':85.0,'16':85.0,'17':80.0,'18':80.0,'19':80.0,'20':80.0,'21':80.0,'22':80.0,'23':80.0,'24':80.0},
'AOD1':{'1':75.0,'2':75.0,'3':75.0,'4':75.0,'5':80.0,'6':80.0,'7':80.0,'8':80.0,'9':95.0,'10':95.0,'11':95.0,'12':95.0,'13':85.0,'14':85.0,'15':85.0,'16':85.0,'17':85.0,'18':95.0,'19':95.0,'20':95.0,'21':80.0,'22':80.0,'23':80.0,'24':80.0},
'AOD2':{'1':75.0,'2':75.0,'3':75.0,'4':75.0,'5':80.0,'6':80.0,'7':80.0,'8':80.0,'9':95.0,'10':95.0,'11':95.0,'12':95.0,'13':85.0,'14':85.0,'15':85.0,'16':85.0,'17':85.0,'18':95.0,'19':95.0,'20':95.0,'21':80.0,'22':80.0,'23':80.0,'24':80.0},
'LF1':{'1':35.0,'2':35.0,'3':35.0,'4':35.0,'5':45.0,'6':45.0,'7':20.0,'8':20.0,'9':45.0,'10':45.0,'11':45.0,'12':45.0,'13':25.0,'14':25.0,'15':25.0,'16':25.0,'17':25.0,'18':45.0,'19':45.0,'20':30.0,'21':30.0,'22':30.0,'23':30.0,'24':30.0},
'LF2':{'1':35.0,'2':35.0,'3':35.0,'4':35.0,'5':45.0,'6':45.0,'7':20.0,'8':20.0,'9':45.0,'10':45.0,'11':45.0,'12':45.0,'13':25.0,'14':25.0,'15':25.0,'16':25.0,'17':25.0,'18':45.0,'19':45.0,'20':30.0,'21':30.0,'22':30.0,'23':30.0,'24':30.0},
'CC1':{'1':50.0,'2':50.0,'3':50.0,'4':50.0,'5':60.0,'6':60.0,'7':55.0,'8':55.0,'9':60.0,'10':60.0,'11':60.0,'12':60.0,'13':70.0,'14':70.0,'15':75.0,'16':75.0,'17':75.0,'18':60.0,'19':70.0,'20':70.0,'21':50.0,'22':50.0,'23':50.0,'24':50.0},
'CC2':{'1':50.0,'2':50.0,'3':50.0,'4':50.0,'5':60.0,'6':60.0,'7':55.0,'8':55.0,'9':60.0,'10':60.0,'11':60.0,'12':60.0,'13':70.0,'14':70.0,'15':75.0,'16':75.0,'17':75.0,'18':60.0,'19':70.0,'20':70.0,'21':50.0,'22':50.0,'23':60.0,'24':60.0}
}
stage2units = {'1': {'EAF': 2}, '2': {'AOD': 2}, '3': {'LF': 2}, '4': {'CC1': 1, 'CC2': 1}}trans_time = {'TR_S1': 10, 'TR_S2': 4, 'TR_S3': 10}
trans_time_max = {'TR_S1': 240, 'TR_S2': 240, 'TR_S3': 120}
setup_time = {'CC1': 70, 'CC2': 50}
equip2process_time['EAF'] = equip2process_time.pop('EAF1')
equip2process_time.pop('EAF2')
equip2process_time['AOD'] = equip2process_time.pop('AOD1')
equip2process_time.pop('AOD2')
equip2process_time['LF'] = equip2process_time.pop('LF1')
equip2process_time.pop('LF2')
model = 'rtn1'
rtn_t0 = 15price_energy = []
for price in energy_price:
price_energy = price_energy + [int(price)] * int(60 / rtn_t0) #prices for 96 time slots
num_t = 24 * 60 / rtn_t0 #total no of time slots i.e 96
num_stage = len((stage2units).keys()) #total no of stages i.e 4
num_groups = len((group2heats).keys()) #total no of groups i.e 6
num_heats = 0
for group, heats in group2heats.items():
num_heats += len(heats) #total no of heats i.e 24equip2num = {} #calculating each eqipment has how many units
for unit2num in stage2units.values():
for unit, num in unit2num.items():
equip2num[unit] = numres_cat2idx = OrderedDict()
r_idx = 1# equipment resources #creating a ordered dict where at 1st position there is resource EAF, 2nd position there is resource AOD , 3rd position there is resource LF, 4th position there is resource CC1, 5th position there is resource CC2
for stage in range(1, num_stage + 1):
for unit in (stage2units)[str(stage)].keys():
res_cat2idx[unit] = [r_idx]
r_idx = r_idx + 1
# intermediate products and final products, A - after, B - before
for stage in range(1, num_stage+1):
res_cat2idx['H_A_S%s' % stage] = range(r_idx, r_idx + num_heats) #'H_A_S1'=EAd1 (after stage1) 'H_B_S4'=EAs4(before stage4) #'H_A_S4'=intermediate product after stage 4 i.e final products heats
r_idx = r_idx + num_heats
if int(stage) == 1:
continue
res_cat2idx['H_B_S%s' % stage] = range(r_idx, r_idx + num_heats)
r_idx =r_idx + num_heats
# energy resource
res_cat2idx['EN'] = [r_idx]
res_num = r_idx[res_cat2idx, res_num]
tasks = OrderedDict()
i_idx = 1for stage in range(1,num_stage + 1): #calculates in order how many processing and transfer tasks are there. task_num is 156 for only 3 stages.(24 heats * 6 types of tasks for 3 stages(3 processsng and 3 transfer))
for unit in (stage2units)[str(stage)].keys(): #CC1 and CC2 tasks are not counted in task_num
if int(stage) == 4:
tasks[unit] = range(i_idx, i_idx + num_groups)
i_idx = i_idx + num_groups
else:
tasks[unit] = range(i_idx, i_idx + num_heats)
i_idx = i_idx + num_heats
tasks['TR_S%s' % stage] = range(i_idx, i_idx + num_heats)
i_idx =i_idx + num_heats
task_num = i_idx - 1
[tasks, task_num]task_duration = dict()
task_cleanup_duration = dict()for heat in range(1, num_heats + 1): #for heat 1, EAF,AOD,LF,TR_S1,TR_S2,TR_S3 , similar for heat 2...so on...
for task_type in ['EAF', 'AOD', 'LF']:
task_duration[tasks[task_type][heat - 1]]=int(math.ceil(float((equip2process_time)[task_type][str(heat)])/rtn_t0))
for task_type in ['TR_S1', 'TR_S2', 'TR_S3']:
task_duration[tasks[task_type][heat - 1]]=int(math.ceil(float((trans_time)[task_type])/rtn_t0))
for group in range(1, num_groups + 1):
for task_type in ['CC1', 'CC2']:
cast_time = [(equip2process_time)[task_type][str(heat)] for heat in group2heats[str(group)]]
total_time = float(sum(cast_time))
task_duration[tasks[task_type][group - 1]] = int(math.ceil(total_time/rtn_t0))
total_time = float(sum(cast_time) + setup_time[task_type])
task_cleanup_duration[tasks[task_type][group - 1]] = int(math.ceil(total_time/rtn_t0))
task_duration,task_cleanup_durationrtn_profile = dict()
for res_category, res_idxes in res_cat2idx.items(): #resources=res_cat2idx,created 174 blank dictionaries for 174 res_num
for res_idx in res_idxes:
rtn_profile[res_idx] = dict()for stage in range(1, num_stage + 1):
for unit in (stage2units)[str(stage)].keys(): #[RESOURCE][TASK]'[1][1-25]'(EAF) , '[2][49-73]'(AOD) ,'[3][97-121]'(LF),''[4][145-151]'(CC1) ,'[5][151-157](CC2)'
res = res_cat2idx[unit][0]
for task in tasks[unit]:
rtn_profile[res][task] = [-1] + [0] * (task_duration[task] - 1) + [1]
if stage == 4:
rtn_profile[res][task] = [-1] + [0] * (task_cleanup_duration[task] - 1) + [1]for heat_idx in range(0, num_heats):
# process task generate intermediate heat H_A_S (after stage)
for task_cat, res_cat in [('EAF', 'H_A_S1'), ('AOD', 'H_A_S2'), ('LF', 'H_A_S3')]:
task = tasks[task_cat][heat_idx]
resource = res_cat2idx[res_cat][heat_idx]
rtn_profile[resource][task] = [0] * task_duration[task] + [1]
# process task consume intermediate heat H_B_S (before stage)
for task_cat, res_cat in [('AOD', 'H_B_S2'), ('LF', 'H_B_S3')]:
task = tasks[task_cat][heat_idx]
resource = res_cat2idx[res_cat][heat_idx]
rtn_profile[resource][task] = [-1] + [0] * task_duration[task]
# transfer task transports heat
for task_cat, res_cat in [('TR_S1', ['H_A_S1', 'H_B_S2']),('TR_S2', ['H_A_S2', 'H_B_S3']),('TR_S3', ['H_A_S3', 'H_B_S4'])]:
task = tasks[task_cat][heat_idx]
resource1 = res_cat2idx[res_cat[0]][heat_idx]
rtn_profile[resource1][task] = [-1] + [0] * task_duration[task]
resource2 =res_cat2idx[res_cat[1]][heat_idx]
rtn_profile[resource2][task] = [0] * task_duration[task] + [1]heat_consume_time_in_group = dict()
heat_generate_time_in_group = dict()for unit in (stage2units)['4']:
heat_consume_time_in_group[unit] = dict()
heat_generate_time_in_group[unit] = dict()
for group, heats in group2heats.items():
task = tasks[unit][int(group)-1]
duration = 0
for heat in heats:
consume_time = int(math.floor(duration/rtn_t0))
resource = res_cat2idx['H_B_S4'][heat - 1]
rtn_profile[resource][task] = [0] * (task_duration[task] + 1)
rtn_profile[resource][task][consume_time] = -1
duration += (equip2process_time)[unit][str(heat)]
generate_time = int(math.ceil(duration/rtn_t0))
resource = res_cat2idx['H_A_S4'][heat - 1]
rtn_profile[resource][task] = [0] * (task_duration[task] + 1)
rtn_profile[resource][task][generate_time] = 1
heat_consume_time_in_group[unit][heat] = consume_time
heat_generate_time_in_group[unit][heat] = generate_timetime = list()
for x in range(1,int(num_t)+1):
time.append(x)res_list = [res for cat, cat_res in res_cat2idx.items() for res in cat_res]
task_list=[task for cat, cat_task in tasks.items() for task in cat_task]steel=gp.Model('steel RTN')
Resources=steel.addVars(res_list,time,name="Resources")
Tasks=steel.addVars(task_list,time,name="Tasks")res_list1 = [res for cat, cat_res in res_cat2idx.items() if cat != 'EN' for res in cat_res]
Resourcebalance=steel.addConstrs((Resources[res,t]==Resources[res,time[time.index(t)-1]] +
gp.quicksum(rtn_profile[res][task][theta]* Tasks[task,t-theta]
for task in task_list if task not in rtn_profile[res] continue duration = len(rtn_profile[res][task]) - 1
for theta in range(0, min(duration, t) + 1))
for t in time
for res in res_list1),"Resourcebalance")Note: duration = len(rtn_profile[res][task]) - 1 will calculate the length, but for those task that are not in rtn_profile[res] , that should be skipped otherwise the loop will break , so I have kept continue in it.
I can attach below the snippet written in cplex:(I am trying to apply same logic in Gurobi)
"
res_list = [res for cat, cat_res in resources.items() if cat != 'EN' for res in cat_res]
task_list = [task for cat, cat_task in tasks.items() for task in cat_task]
row_num = num_t * len(res_list)x_vals = []
x_rows = []
x_cols = []
y_vals = []
y_rows = []
y_cols = []
b_rhs = [0] * int(row_num)row_ = 0
for res in res_list:
for t_ in range(0, int(steel_rtn.num_t)):
y_rows.append(row_)
y_cols.append(self.pos_y_res_t(res, t_))
y_vals.append(-1)
if t_ >= 1:
y_rows.append(row_)
y_cols.append(self.pos_y_res_t(res, t_ - 1))
y_vals.append(1)
else:
b_rhs[row_] += -cal_initial_resource(res)
for task in task_list:
if task not in rtn_profile[res]:
continue
duration = len(rtn_profile[res][task]) - 1
for theta in range(0, min(duration, t_) + 1):
x_rows.append(row_)
x_cols.append(self.pos_x_task_t(task, t_ - theta))
x_vals.append(rtn_profile[res][task][theta])
row_ += 1
assert(row_num-1 == max(x_rows))
return [coo_matrix((x_vals, (x_rows, x_cols)), shape=(int(row_num), int(self.num_x))),
coo_matrix((y_vals, (y_rows, y_cols)), shape=(int(row_num), int(self.num_y))),
b_rhs]"
0 -
Please note that the idea of a minimal reproducible example is that it is quite short.
You should construct your Resourcebalance constraints and the sum expressions 1 by 1. Then, you can use a \(\texttt{continue}\) statement and assign different values to the duration variable.
Resourcebalance = {}
for res in res_list1:
for t in time:
sumexpression = gp.LinExpr(0)
for task in task_list:
if task not in rtn_profile[res]:
continue
duration = len(rtn_profile[res][task]) - 1
for theta in range(0, min(duration, t) + 1):
sumexpression.add(rtn_profile[res][task][theta]* Tasks[task,t-theta] )
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] ==
Resources[res,time[time.index(t)-1]] + sumexpression,
name="Resourcebalane[%s,%d]"%(res,t))0 -
Many thanks for the reply.
Although what I understood from above is , if there is an expression to be evaluated inside quicksum, it should be treated individually and then it shall be added.
Also, you created an empty dictionary of Resourcebalance, generally for constraint, we don't do so. Am I right?
0 -
Although what I understood from above is , if there is an expression to be evaluated inside quicksum, it should be treated individually and then it shall be added
Yes, correct. It is recommended to construct more complex terms, which involve many \(\texttt{ifs}\) and other checks, "by hand" instead of pushing everything inside a quicksum. In your particular case, you have to avoid quicksum due to the \(\texttt{duration}\) assignment.
Also, you created an empty dictionary of Resourcebalance, generally for constraint, we don't do so. Am I right?
The dictionary is initially empty and is filled with constraints in the \(\texttt{for}\)-loop via
Resourcebalance[res,t] = ...
So by accessing \(\texttt{Resourcebalance[res,t]} \) you get the respective Resourcebalance constraint.
The same would have been done automatically by addConstrs but we had to avoid it due to the quicksum complication.
0 -
I appreciate your detailed explanation. It made me clear few points that I was struggling to figure it out. I had gone though MILP tutorial on Gurobi website but these little points were missed out & it created lots of struggle in building a model.
I am sure there will be many other minute things to be looked upon. But thank you for the same.
"The same would have been done automatically by addConstrs but we had to avoid it due to the quicksum complication"
- Means, We cant use addConstrs when using quicksum, instead we should use addConstr ?
0 -
- Means, We cant use addConstrs when using quicksum, instead we should use addConstr
No, you can use quicksum in an addConstrs call, see the diet.py example.
The reason why you have to use addConstr and construct the sum by hand instead of using quicksum is that you have an assignment within your construction. In the \(\texttt{task}\) \(\texttt{for}\)-loop, you assign a different value to \(\texttt{duration}\) and use this value in the next \(\texttt{for}\)-loop.
0 -
Okay, I understood the point, I tried formulating using addConstrs for the below equation & it works:
Many thanks!
Taskexecution1=steel.addConstrs((gp.quicksum(Tasks[task,t] for t in time for task in group2tasks) ==1 for group in range(1, num_groups + 1)),name="Taskexecution1")
0 -
Surprising Output of Gurobi...!!!
Here my Constraint looks like this:
& After "print(steel.getAttr("ConstrName"))", it looks completely different. How?
0 -
The output is not surprising at all. You first print the \(\texttt{Resourcebalance}\) dictionary and then print all constraint names in your model. Multiple constraints have the same name in your case. This is because your indentation is wrong. The line
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] ==
Resources[res,time[time.index(t)-1]] + sumexpression,
name="Resourcebalane[%s,%d]"%(res,t))is executed in the \(\texttt{for theta}\) loop. This happens often when copy pasting lines from one editor to another and one has to be careful with indentation in Python
Resourcebalance = {}
for res in res_list1:
for t in time:
sumexpression = gp.LinExpr(0)
for task in task_list:
if task not in rtn_profile[res]:
continue
duration = len(rtn_profile[res][task]) - 1
for theta in range(0, min(duration, t) + 1):
sumexpression.add(rtn_profile[res][task][theta]* Tasks[task,t-theta] )
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] ==
Resources[res,time[time.index(t)-1]] + sumexpression,
name="Resourcebalane[%s,%d]"%(res,t))0 -
Hi Jaromil,
I doubleckeched the indentation issue, but I think it isnt working.
Now my loop is not running for time:(It is just picking up the last value)
0 -
Did you change anything about your \(\texttt{time}\) list? Please print the \(\texttt{time}\) list to check what is in there.
It is very hard to help when you change your code a lot and then just post screenshots. Using your code from your other post and adjusting the indentation of \(\texttt{Resourcebalance}\) works as expected.
0 -
NO, My time list is the same:
time = list()
for x in range(0,int(num_t)):
time.append(x)res_list1 = [res for cat, cat_res in res_cat2idx.items() if cat != 'EN' for res in cat_res]
Resourcebalance = {}
for res in res_list1:
for t in time:
#if t==0:
#continue
sumexpression = gp.LinExpr(0)
for task in task_list:
if task not in rtn_profile[res]:
continue
duration = len(rtn_profile[res][task]) - 1
for theta in range(0, min(duration, t) + 1):
sumexpression.add(rtn_profile[res][task][theta]* Tasks[task,t-theta] )
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == Resources[res,time[time.index(t)-1]] + sumexpression, name="Resourcebalane[%s,%d]"%(res,t))0 -
In your screenshot, the line
Resourcebalance[res,t] = ...
is incorrectly indented. In your last comment, it is indented correctly. Please double check your indentation when you copy paste from one editor to another.
0 -
Jaromil, my steel.lp file is creating but it is empty.
Can you suggest why it is happening?
0 -
Below is my steel.lp file after returning warnings:
0 -
Below is my steel.lp file after returning warnings:
This is not Gurobi output.
How did you generate the file? Do you have permissions to write and save files on your machine?
0 -
I just wrote steel.lp like this
What do you mean by permissions?
0 -
The output looks good. How did you try to open the file? Please use a standard text editor instead of Jupyter Notebook.
What do you mean by permissions?
Are you working on your own machine, i.e., you have all writing permissions, or on an external one where your file permissions might be limited.
0 -
I opened the file in Jupyter notebook.
I am working with my University machine, I am unsure about the permissions.
0 -
I opened the file in Jupyter notebook.
Could you please try opening the file in a standard text editor.
0
サインインしてコメントを残してください。
コメント
25件のコメント