Using a for loop for a single or specific value in variable list in constraint formulation
Answeredtime = list()
for x in range(0,int(num_t)):
time.append(x)
Resources=steel.addVars(res_list,time,vtype=GRB.INTEGER,name="Resources")
Tasks=steel.addVars(task_list,time,vtype= GRB.BINARY,name="Tasks")
My res_list is also a list
My time is also a list
Now When I am writing this,
Endresources=steel.addConstrs((gp.quicksum(Resources[res,t] for t in x)==1 for res in res_list3),name="Endresources")
where
x=[]
x.append(time[95])
But In Constraints it is showing,
Endresources[150],Endresources[151]........
Instead it should show:
Endresources[150,95],Endresources[151,95]......
How can I use "time" in Constraint instead of "x"
Endresources=steel.addConstrs((gp.quicksum(Resources[res,t] for t in x )==1 for res in res_list3),name="Endresources")
-
-
Gurobi Staff
Why would we require TWO constraints to execute this single constraint?
You need to formulate the OR statement in a way that is accepted by Gurobi. Gurobi does not directly accept a constraint of the form expression1 OR expression2 == 1. Thus, you have to reformulate it in a way that is accepted by Gurobi. This is what I did in my previous post.
Also, if we are adding a new binary variable, in the model who will decide if b=0 or b=1?
That's the beauty of mathematical optimization. The algorithm will decide how to set \(b\). You can imagine that in worst case one could just try all possible values for \(b\) to find the optimal solution.
For such type of constraint, is it compulsory to use binary variable?
Yes, for OR, AND, XOR and other "logical constraints", it is necessary to introduce auxiliary binary variables of some sort. There are standard techniques to do so but it is always better to exploit some knowledge to make the formulation as tight as possible, which is what we did here. We use only 1 auxiliary binary variable and 2 constraints instead of 1.
And yes, the posts are similar. I will delete them.
Thank you.
0 -
-
Hi Jaromil,
Huge thanks for this wonderful explanation. This seems very detailed one.
I understood the logic behind it.
constr1[1]: Tasks[145,0] + Tasks[145,1] + Tasks[145,2] + Tasks[145,3] ............
+ Tasks[145,92] + Tasks[145,93] + Tasks[145,94] + Tasks[145,95] - b[1]
= 0
constr2[1]: Tasks[151,0] + Tasks[151,1] + Tasks[151,2] + Tasks[151,3].............
+ Tasks[151,92] + Tasks[151,93] + Tasks[151,94] + Tasks[151,95] + b[1]
= 1So now these both are my separate constraints?
i.e b [1] can be 0 or 1 and depending upon that my any one constraint will be executed.
0 -
-
-
So now these both are my separate constraints?
i.e b [1] can be 0 or 1 and depending upon that my any one constraint will be executed.
Correct.
0 -
-
Hi Jaromil,
Below is the log file Gurobi has generated, I want to understand what it means.
I am showing initial few lines & last few lines. Could you please let me know what it actually means?
This is the link of whole log file:
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 10 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 23700 rows, 31686 columns and 160840 nonzeros
Model fingerprint: 0xeaffa06e
Variable types: 96 continuous, 31590 integer (14982 binary)
Coefficient statistics:
Matrix range [2e-01, 2e+01]
Objective range [3e+01, 9e+01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+01]
Presolve removed 16519 rows and 16544 columns
Presolve time: 0.16s
Presolved: 7181 rows, 15142 columns, 59961 nonzeros
Variable types: 0 continuous, 15142 integer (8109 binary)
Root relaxation: objective 6.448684e+04, 9074 iterations, 2.58 seconds (6.46 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 64486.9769 0 5634 - 64486.9769 - - 2s
0 0 65108.1536 0 5742 - 65108.1536 - - 6s
0 0 65118.0083 0 5700 - 65118.0083 - - 6s
0 0 65118.1585 0 5704 - 65118.1585 - - 6s
0 0 65118.2783 0 5695 - 65118.2783 - - 6s
0 0 66678.6542 0 5731 - 66678.6542 - - 10s
0 0 66693.8042 0 5701 - 66693.8042 - - 10s
0 0 66693.8234 0 5695 - 66693.8234 - - 10s
0 0 66773.8467 0 5851 - 66773.8467 - - 11s
0 0 66779.7182 0 5860 - 66779.7182 - - 11s
0 0 66779.8701 0 5859 - 66779.8701 - - 11s
0 0 66931.5449 0 5761 - 66931.5449 - - 13s
0 0 66931.9096 0 5764 - 66931.9096 - - 13s
0 0 66995.7259 0 5961 - 66995.7259 - - 14s
0 0 67000.6985 0 5968 - 67000.6985 - - 14s
0 0 67000.9042 0 5981 - 67000.9042 - - 14s
0 0 67059.7897 0 5908 - 67059.7897 - - 16s
0 0 67060.4346 0 6044 - 67060.4346 - - 16s
0 0 67088.9935 0 5855 - 67088.9935 - - 17s
0 0 67089.1326 0 5865 - 67089.1326 - - 17s
0 0 67100.8836 0 6078 - 67100.8836 - - 17s
0 0 67101.4149 0 6062 - 67101.4149 - - 17s
0 0 67112.8153 0 5984 - 67112.8153 - - 18s
H 0 0 216804.50000 67112.8153 69.0% - 18s
0 0 67113.8293 0 6006 216804.500 67113.8293 69.0% - 20s
H 0 0 184673.50000 67113.8293 63.7% - 20s229292 12986 infeasible 102 79015.1632 78815.2107 0.25% 654 21950s
231068 11445 infeasible 102 79015.1632 78830.5837 0.23% 652 21991s
233265 10222 infeasible 107 79015.1632 78845.2081 0.22% 648 22030s
235308 8717 cutoff 111 79015.1632 78857.1854 0.20% 644 22062s
237424 6900 78982.8200 133 647 79015.1632 78873.1162 0.18% 640 22092s
239744 5332 cutoff 93 79015.1632 78889.0746 0.16% 636 22115s
241610 3488 cutoff 112 79015.1632 78906.7292 0.14% 632 22139s
243783 1712 cutoff 127 79015.1632 78931.0110 0.11% 628 22160s
246168 0 cutoff 98 79015.1632 78953.0283 0.08% 623 22170s
Cutting planes:
Gomory: 31
Cover: 805
Implied bound: 7
Projected implied bound: 11
Clique: 15
MIR: 79
StrongCG: 24
Flow cover: 198
Inf proof: 23
Zero half: 249
RLT: 29
Relax-and-lift: 64
Explored 248217 nodes (153438435 simplex iterations) in 22171.08 seconds (15367.39 work units)
Thread count was 12 (of 12 available processors)
Solution count 10: 79015.2 79015.2 79015.2 ... 79052.7
Optimal solution found (tolerance 1.00e-04)
Best objective 7.901516666667e+04, best bound 7.901516666667e+04, gap 0.0000%0 -
-
-
The logfile means that Gurobi was able to successfully solve your model to optimality. Indeed, since the final gap is 0%, it is the true global solution, i.e., there is guaranteed no better solution than the one reported. The optimal solution value is 7.90151666e+4.
For more details, please refer to the MIP Logging section of the documentation.
0 -
-
Hey,
Jaromił Najman , Is it possible to have an in-between infeasible solution in the "Current Node" Column?
i.e: objective of the associated relaxation value?
229292 12986 infeasible 102 79015.1632 78815.2107 0.25% 654 21950s
0 -
-
-
A current node can be infeasible because it is only a subpart of the original model. It does not mean that the whole problem is infeasible. To get a better understanding of how a B&B algorithm work, I would recommend to have a look at our Learning series, in particular the B&B introduction.
0 -
-
Many Thanks Jaromil for the explanation.
However, I find that my algorithm takes much time (5 hrs approx.) to finish.
The same model on a different PC with same specifications takes around 30 minutes.
Why does my MILP take so longer time?
I haven't used any parameter.
0 -
-
-
However, I find that my algorithm takes much time (5 hrs approx.) to finish.
The same model on a different PC with same specifications takes around 30 minutes.
The underlying hardware plays a big role when running computations. This does not apply to optimization algorithms only. It is expected that a more powerful PC takes only a fraction of time compared to a personal notebook.
Why does my MILP take so longer time?
Note that MILPs are not easy problems and it is in general not expected that a given problem is solved within a few minutes. Depending on your application, you could increase the desired MIPGap. The MIPGap describes the relative difference between the objective value \(z_P\) of the best feasible solution found so far and the best proven bound \(z_D\). A MIPGap=0.01 would mean that the algorithm would stop once a relative gap of 1% is proven, i.e., \(z_P\) could be at most improved by 1% if at all. For more complex problems and most application a MIPGap larger than the default gap of 0.01%, e.g., 1%, is acceptable.
0 -
-
Many thanks Jaromil for this explanation.
With your time & explanation, I was able to make the model feasible for the very first time ever.
However, the results are not as expected.
So, I was analyzing some constraints & there arose 1 doubt i.e obvious that how the constraint is working:
The below constraint:
Resourcebalance = {}
for res in res_list1:
for t in time:
if t == 0:
if 1 <= res and res <= 3:
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 2, name="Resourcebalane[%s,%d]"%(res,t))
elif 4<= res and res<=5:
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 1, name="Resourcebalane[%s,%d]"%(res,t))
#else:
# Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 0, name="Resourcebalane[%s,%d]"%(res,t))
else:
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))As you can see that the indent of Resourcebalance constraint is below for task...loop, i.e outside theta loop,
But, my sumexpression inside the constraint formulation contains theta, how it will take the values of theta then?
I mean to say how the constraint is executing?
0 -
-
-
But, my sumexpression inside the constraint formulation contains theta, how it will take the values of theta then?
The sumexpression object is added to Resources in the Resourcebalance constraint. Theta is only used to construct the sumexpression object with correct indices. Did you have a look at the Resourcebalance constraints in the LP file and checked whether they are correct?
I mean to say how the constraint is executing?
What exactly do you mean by "executing"? A constraint is not executing. A constraint is given and will be satisfied by any feasible point.
0 -
-
Yes, I checked my resource Balance constraints in LP file, and it seems correct to me.
But the output does not seem to be correct. Can you please check if my indent is correct in this constraint?
0 -
-
-
Yes, I checked my resource Balance constraints in LP file, and it seems correct to me.
Then there is no reason to assume that something is wrong with it.
But the output does not seem to be correct. Can you please check if my indent is correct in this constraint?
Which output do you mean?
The indent looks correct.
0 -
-
By the Output I mean that When I saw my values of Resource variables after optimization, I can see that Resource balance constraint is violated.
0 -
-
-
Can you show the corresponding solution point values and the respective constraint? How did you access the solution point values?
0 -
-
Hi Jaromil,Jaromił Najman
Apologies for the late reply as something unfortunate happened.
I will show you why the constraint is not working:
I have made one change in my time list :
time = list()
for x in range(1,int(num_t + 1)):
time.append(x)
#VARIABLES
Resources=steel.addVars(res_list1,time,vtype=GRB.INTEGER,name="Resources")
Tasks=steel.addVars(task_list,time,vtype= GRB.BINARY,name="Tasks")
EnergyResource=steel.addVars(res_list4,time,vtype=GRB.CONTINUOUS,name="EnergyResource")
#CONSTRAINT
Resourcebalance = {}
for res in res_list1:
for t in time:
if t == 0:
if 1 <= res and res <= 3:
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 2, name="Resourcebalane[%s,%d]"%(res,t))
elif 4<= res and res<=5:
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 1, name="Resourcebalane[%s,%d]"%(res,t))
#else:
# Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 0, name="Resourcebalane[%s,%d]"%(res,t))
else:
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)):
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))Resourcebalane[1,1]: Resources[1,1] - Resources[1,96] + Tasks[1,1]
+ Tasks[2,1] + Tasks[3,1] + Tasks[4,1] + Tasks[5,1] + Tasks[6,1]
+ Tasks[7,1] + Tasks[8,1] + Tasks[9,1] + Tasks[10,1] + Tasks[11,1]
+ Tasks[12,1] + Tasks[13,1] + Tasks[14,1] + Tasks[15,1] + Tasks[16,1]
+ Tasks[17,1] + Tasks[18,1] + Tasks[19,1] + Tasks[20,1] + Tasks[21,1]
+ Tasks[22,1] + Tasks[23,1] + Tasks[24,1] = 0Here instead of Resources[1,96] Should not it should be Resources[1,0] ?
Only then the constraint will execute properly.
I think the formulation is correct.But dont know why it is happening.I realize may be its because the time list starts from 1, but then I need to give the initial values at 0, Is it possible to do so?
0 -
-
-
>Here instead of Resources[1,96] Should not it should be Resources[1,0] ?
It cannot be t = 0 because 0 is not part of your time list as you already said. Note that with the new time definition there is no variable Resources[1,0], because your time list starts with 1.
What is the issue with time = 0 in your formulation? I don't see any issue with having t starting at 0. In general, it does not matter at what value t starts as long as you adjust all indices appropriately. Having t starting at 0 is often the most convenient way.
Note that you are dealing with t=0 in a special way anyway, see
if t == 0:
#...
else:
#...Note that in your current implementation with t starting at 1, the constraints
Resourcebalance[1,0]: Resources[1,0] = 2
...are not present because there is no t = 0.
0 -
-
Hi Jaromil,
What I expect is I have to give initial value i.e Resources [1,0] but the optimization should start from time step 1. As you can see below, I am using Energy resource variable in my obj assignment, is it possible to start the optimization from time step 1 or do I need to make some changes in Resource balance constraint?
res = res_cat2idx['EN'][0]
Energyusage = {}
for ti in time:
sumexpression = gp.LinExpr(0)
for task_cat, task_list1 in tasks.items():
if 'TR' in task_cat: # transportation has no energy consumption
continue
for task in task_list1:
duration = task_duration[task]
for theta in range(0, min(duration, ti) + 1):
sumexpression.add(rtn_profile[res][task][theta]* Tasks[task,ti-theta] )
Energyusage[res,ti] = steel.addConstr(EnergyResource[res,ti] ==sumexpression,name="Energyusage[%s,%d]"%(res,ti))
obj=gp.quicksum(price_energy[t]*EnergyResource[res,t] for res in reslist2 for t in time)0 -
-
-
Hi Margi,
What I expect is I have to give initial value i.e Resources [1,0] but the optimization should start from time step 1.
In this case, you should still define your time list such that it holds 0 so you define the Resource[res,0] variables.
You can then just skip the t=0 case when defining any constraint that should not use it, e.g.,
for ti in time:
if ti == 0: # skip the t = 0 case
continue
sumexpression = gp.LinExpr(0)
# ...0 -
-
Hi Jaromil,
I understood your point but if I do so, will other constraints get affected?
By that I mean,
Taskexecution=steel.addConstrs((gp.quicksum(Tasks[task,t] for t in time)==1 for task in heat2tasks),name="Taskexecution")
Task [1,0] should never take place, I mean Tasks should start from time step 1 but as time list starts from 0, will it affect this constraint even though I am optimizing from time step 1?
0 -
-
-
Of course, every constraint using time will be affected. That is why I said that you will have to skip it. In the particular case, you posted, you could use an if condition
Taskexecution=steel.addConstrs((gp.quicksum(Tasks[task,t] for t in time if t != 0)==1 for task in heat2tasks),name="Taskexecution")
Alternatively, you could have a second list called time1 which starts with 1 and use this list whenever you want to avoid using t=0. So you would have
time = [0,1,2,3,...,95]
time1 = [1,2,3,4,...,95]You can then still define all your variables over the list time but use time1 to construct constraints where you explicitly have to avoid t=0.
0 -
-
Hi Jaromil,
Thanks I will try to implement this.
In between there is one weird error i am facing,
time1 = list()
for x in range(1,int(num_t + 1)):
time1.append(x)
energy_price = [33.28, 30.00, 29.15, 28.49, 34.66, 50.01, #14 sept nord pool day ahead SE1
71.52, 77.94, 81.97, 87.90, 92.76, 94.98,
92.31, 90.03, 90.09, 87.44, 85.37, 79.97,
79.92, 77.83, 76.28, 65.06, 53.07, 34.16]
price_energy = []
for price in energy_price:
price_energy = price_energy + [int(price)] * int(60 / rtn_t0) #prices for 96 time slots
reslist2=[]
reslist2.append(res_cat2idx['EN'][0])
obj=gp.quicksum(price_energy[t]*EnergyResource[res,t] for res in reslist2 for t in time1)
Here price_energy and time1 has same length , still how does it happen?
0 -
-
-
You should try debugging such issues via removing parts of the code.
Please remove \(\texttt{price_energy[t]}\) from \(\texttt{obj}\) and re-run. If it works then, something is off with your \(\texttt{price_energy}\) list. You should then write a simple \(\texttt{for}\)-loop over time1 and add some print statement to check which index is failing.
If removing \(\texttt{price_energy[t]}\) from \(\texttt{obj}\) does not solve the issue, then you should have a look at your definition of \(\texttt{EnergyResource[res,t]}\).
0 -
-
Hi Jaromil,
I tried every possible way but somehow its causing errors.
Is it any possible way that I start my time step from 1 and give my initial values at 0 specifically, in the Resource balance constraint?
I asked my seniors and they suggested, its simple: "you start your optimization from time step 1 and give initial value at 0"
Is there any way of doing so?
0 -
-
-
Is it any possible way that I start my time step from 1 and give my initial values at 0 specifically, in the Resource balance constraint?
You are already doing this. Please look at your code:
Resourcebalance = {}
for res in res_list1:
for t in time:
if t == 0:
if 1 <= res and res <= 3:
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 2, name="Resourcebalane[%s,%d]"%(res,t))
elif 4<= res and res<=5:
Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 1, name="Resourcebalane[%s,%d]"%(res,t))
#else:
# Resourcebalance[res,t] = steel.addConstr(Resources[res,t] == 0, name="Resourcebalane[%s,%d]"%(res,t))
else:
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)):
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))For t=0, you are setting constraints of the form Resources[res,0] = some value. Those are initial values. Then, you are adding all other constraints starting from t=1 (the else case). Now, you have to adjust all your other constraints to follow a similar encoding when required.
0 -
-
HelloJaromił Najman,
I will tell you what is happening, I saw the .sol file. I am attaching one file for your reference if you can have a look at it. I am sure you will understand what difficulty I am facing only because of this constraint.
The constraint that has been formulated is correct but it needs little modifications but I dont know how to do it.
First of all, Task should start from time step 1. i.e Task[1,0] should never take place , How do I define this ?
Because this is a production problem, I have used time step 0 just to give initial values, not for the start of optimisation horizon.
If I am using the constraint in the similar way,
below file is the result of that which is wrong.
Please kindly have a look.
Should I use a different time list or what kind of modifications shall I make so that this criterion is achieved?
If it's not clear please let me know.
-kind regards,
Margi
0 -
-
-
Hi Margi,
As mentioned earlier, you can avoid specific indices via \(\texttt{if}\)-conditions
#[...]
else:
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)):
if (t-theta) != 0:
sumexpression.add(rtn_profile[res][task][theta]* Tasks[task,t-theta] )
#[...]This way, you avoid using variables Tasks[task,0] in all constraints, while still preserving its formulation.
Is this what you are looking for?
Best regards,
Jaromił0 -
-
Hi Jaromil,
First of all, huge thanks for helping me out. Please help me in this last stage as I feel low now by building such a huge model. You have been a moral support really, Many thanks.
Also, I tried doing so, but I am hereby sharing a link of .lp file, which seems not correct. Please have a look.
https://drive.google.com/file/d/1El5fNk39oDyvz0fNgbBnhvZBg6husGhv/view?usp=sharing
Also I am sharing a reference file of the same constraint which is formulated via cplex solver.Please have a look. I can give my 200 % in this formulation of constraint but please help me.
https://drive.google.com/file/d/1El5fNk39oDyvz0fNgbBnhvZBg6husGhv/view?usp=sharing
- Many Many thanks Jaromil.
0 -
-
-
Hi Margi,
Both links point to the same file.
Could please you also share the code you used to generate this Gurobi LP file (it it's OK for you)? This would make it way easier for me to see what might be going wrong.
Do I understand correctly that it's the Resourcebalance constraints which are not correctly formulated?
0 -
-
Hello Jaromił Najman,
I am sharing link again:
1. Link to the code:
https://drive.google.com/file/d/1uTfz0DbcQUz-EAY7AmaLT6zAWKi9y5kc/view?usp=sharing
2 Link to the reference file of the same constraint which is formulated via cplex solver
"Do I understand correctly that it's the Resourcebalance constraints which are not correctly formulated?"
Yes, the problem lies in 3 steps I think:
- What should be the time list start from?
- I want my all tasks starting from t=1 but as we can see the formulation of resource balance constraint, it has term Task[1,0], How to deal with it?
- objective function should definitely start from t=1
Jaromil, once you will see the reference file too I am sure you will understand what's going wrong.
Again,huge thanks for helping me, I think I reached feasibility but as I am not getting expected results of variables, its because of this constraint.
Thanks from bottom of my heart.
0 -
Post is closed for comments.
Comments
210 comments