Question about linearization
Answered
Hello,
This might be a simple question, but I want to really understand what is going on behind the scene. I'm trying to linearize a product between a binary and a continuous variable. The binary variable also forms a clique, as such, I write the product as the continuous variable y times the sum over the clique, which is still a binary value. Also note that under any feasible conditions the value of y should be zero, regardless of optimality.
Overall, the linerization technique is according to this reference. http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-continuous-and-binary.html
This post emphasizes on the choice of big-M. So since I know y should be around zero, I try values from 0.01 to 1.0 (+ and -) and for anything below 0.05 (approx) model becomes infeasible and if I set M to be a bit higher then it works. Can someone explain what might be happening under the hood?
Regards
Buddi
-
Hi Buddi,
Could you please post the exact bilinear term with variable bounds and the linearization you implemented? Please note that the link you posted does not work.
Best regards,
Jaromił0 -
Hello Jaromil,
Here is the particular linearized constraint. self is an instance of a Gurobi model. _x is binary, _z is continuous in the range [0, 1], e_high and e_low are the big-M's I play with.
self.addConstrs(self._z[i, j] <= e_high * quicksum(self._x[i, j, k] for k in self._types if k == 'elect')
for i in self._locations
for j in self._locations
)self.addConstrs(self._z[i, j] >= e_low * quicksum(self._x[i, j, k] for k in self._types if k == 'elect')
for i in self._locations
for j in self._locations
)self.addConstrs(self._z[i, j] <= (self._E[i, j] * self._b / self._efficiency + self._soc[i] - self._soc[j]) -
e_low * (1 - quicksum(self._x[i, j, k] for k in self._types if k == 'elect'))
for i in self._locations
for j in self._locations
)self.addConstrs(self._z[i, j] >= self._E[i, j] * self._b / self._efficiency + self._soc[i] - self._soc[j] -
e_high * (1 - quicksum(self._x[i, j, k] for k in self._types if k == 'elect'))
for i in self._locations
for j in self._locations
)0 -
I also attach a screenshot of the post.
0 -
So you are trying to linearize the term
(self._E[i, j] * self._b / self._efficiency + self._soc[i] - self._soc[j]) *
quicksum(self._x[i, j, k] for k in self._types if k == 'elect')Is this correct?
If yes, then it would be way easier to handle with the linearization if you introduce auxiliary variables for the terms
aux_cont == (self._E[i, j] * self._b / self._efficiency + self._soc[i] - self._soc[j])
aux_int == quicksum(self._x[i, j, k] for k in self._types if k == 'elect')With the above auxiliary variables, you can see that \(\texttt{aux_int}\) is an integer variable with range \(\{0,\dots,\sum_{k} 1\}\) and \(\texttt{aux_cont}\) is an auxiliary variable with bounds \([e_{low},e_{high}]\). So You actually want to linearize the term
aux_cont * aux_int
which is a multiplication of a continuous and an integer variable and thus the "binary times continuous" linearization does not work.
In general, instead of linearizing bilinear terms by hand, you can let Gurobi do the job. Depending on the rest of your model, it may be possible that you have to set the NonConvex parameter.
Best regards,
Jaromił0 -
Thank you very much. Yes, I missed that when I sum over the binaries, it turns out to be an integer value.
I post the progress above. You see that the dual bound is stuck. By 300s, the incumbent is already very good. I have a strong suspicion that the constraint that showed you before is the main culprit here, because the it is the value of self._b that is almost entirely responsible for the optimality gap at this point. There are no other (main) constraints involving this variable, except the one I write after the log.Solved with dual simplex (primal model)
0 0 1076.00000 0 195 1509.80081 1076.00000 28.7% - 88s
0 0 1076.00000 0 139 1509.80081 1076.00000 28.7% - 91s
H 0 0 1509.7998849 1076.00000 28.7% - 101s
0 0 1076.00000 0 127 1509.79988 1076.00000 28.7% - 104s
H 0 0 1509.7021733 1076.00000 28.7% - 112s
0 0 1076.00000 0 149 1509.70217 1076.00000 28.7% - 114s
0 0 1076.00000 0 135 1509.70217 1076.00000 28.7% - 124s
0 0 1076.00000 0 138 1509.70217 1076.00000 28.7% - 125s
0 0 1076.00000 0 148 1509.70217 1076.00000 28.7% - 138s
0 0 1076.00000 0 163 1509.70217 1076.00000 28.7% - 139s
0 0 1076.00000 0 148 1509.70217 1076.00000 28.7% - 152s
0 0 1076.00000 0 130 1509.70217 1076.00000 28.7% - 156s
0 0 1076.00000 0 174 1509.70217 1076.00000 28.7% - 173s
0 0 1076.00000 0 175 1509.70217 1076.00000 28.7% - 196s
H 0 0 1509.6942936 1076.00000 28.7% - 201s
0 0 1076.00000 0 190 1509.69429 1076.00000 28.7% - 203s
0 0 1076.00000 0 136 1509.69429 1076.00000 28.7% - 218s
0 0 1076.00000 0 182 1509.69429 1076.00000 28.7% - 222s
0 0 1076.00000 0 145 1509.69429 1076.00000 28.7% - 239s
0 0 1076.00000 0 187 1509.69429 1076.00000 28.7% - 245s
0 0 1076.00000 0 138 1509.69429 1076.00000 28.7% - 261s
0 0 1076.00000 0 163 1509.69429 1076.00000 28.7% - 262s
0 0 1076.00000 0 146 1509.69429 1076.00000 28.7% - 278sVariable definitions:self._b = self.addVar(name="b")
self._bi = self.addVar(lb=176, ub=self._battery_size, name="bi")self.addConstr(self._b * self._bi == 1)I watched a techtalk that presented several strategies, but I have some trouble on figuring out how I can improve the formulation of this constraint such that the best bound can be improved faster. Do you have any suggestions?
0 -
I watched a techtalk that presented several strategies, but I have some trouble on figuring out how I can improve the formulation of this constraint such that the best bound can be improved faster. Do you have any suggestions?
Since both variables are continuous, you should try to provide tight lower and upper bounds for variable \(b\). Then, both variables \(b\) and \(bi\) are bounded and the performance should improve.
0 -
Thank you Jaromil. I tried your suggestion.
H 0 0 1547.9138011 1400.00000 9.56% - 279s
0 2 1400.00000 0 35 1547.91380 1400.00000 9.56% - 282s
7 16 1400.00000 3 40 1547.91380 1400.00000 9.56% 223 287s
15 26 1400.00000 4 36 1547.91380 1400.00000 9.56% 144 290s
H 35 46 1547.9054563 1400.00000 9.56% 113 304s
45 122 1400.00000 5 84 1547.90546 1400.00000 9.56% 117 314s
121 348 1400.00000 9 25 1547.90546 1400.00000 9.56% 191 337s
H 315 348 1547.9026785 1400.00000 9.56% 101 337s
H 346 348 1540.4636923 1400.00000 9.12% 97.1 337s
348 997 1400.00000 17 20 1540.46369 1400.00000 9.12% 97.7 387s
H 765 997 1540.4634987 1400.00000 9.12% 60.1 387s
H 861 997 1540.4632668 1400.00000 9.12% 59.6 387s
H 926 997 1510.1708333 1400.00000 7.30% 58.0 387sThe best bound 1400 is what comes immediately after nrel heuristic and it never improves.
0 -
How big is your model? Could you share an LP or MPS file of your model? You can write your model to a file via the write method. Note that uploading files in the Community Forum is not possible but we discuss an alternative in Posting to the Community Forum.
If you cannot share your model, could you please post the first ~20 lines of the log output (or if it's not super long, just the whole log output)?
0 -
Please find the log output below. The LP file is available at https://www.filemail.com/d/xsadvtceoialmei.
Thread count: 10 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 351369 rows, 443590 columns and 2546834 nonzeros
Model fingerprint: 0xd5cd6d37
Model has 201422 quadratic constraints
Model has 680 general constraints
Variable types: 167938 continuous, 275652 integer (275652 binary)
Coefficient statistics:
Matrix range [1e-02, 2e+03]
QMatrix range [9e-01, 6e+01]
Elapsed time for NoRel heuristic: 6s (best bound 1400)
Found phase-1 solution: relaxation 12.0107
Found phase-1 solution: relaxation 12.0089
Found phase-1 solution: relaxation 12.0064
Found phase-1 solution: relaxation 10.0594
Elapsed time for NoRel heuristic: 11s (best bound 1400)
Found phase-1 solution: relaxation 6.03636
Found phase-1 solution: relaxation 6.00048
Found phase-1 solution: relaxation 4
Elapsed time for NoRel heuristic: 18s (best bound 1400)
Found phase-1 solution: relaxation 2
Found phase-1 solution: relaxation 0
Found heuristic solution: objective 1904.0000000
Transition to phase 2
Found heuristic solution: objective 1848.0444444
Elapsed time for NoRel heuristic: 23s (best bound 1400)
Found heuristic solution: objective 1848.0416667
Found heuristic solution: objective 1840.1608854
Found heuristic solution: objective 1838.9569444
Found heuristic solution: objective 1831.4152778
Found heuristic solution: objective 1831.4152503
Found heuristic solution: objective 1786.3637787
Found heuristic solution: objective 1779.7064027
Found heuristic solution: objective 1769.4860215
Found heuristic solution: objective 1738.6023457
Found heuristic solution: objective 1732.4248063
Found heuristic solution: objective 1723.4762510
Found heuristic solution: objective 1723.0244058
Elapsed time for NoRel heuristic: 29s (best bound 1400)
Found heuristic solution: objective 1722.7594768
Found heuristic solution: objective 1721.0040737
Found heuristic solution: objective 1716.1249711
Found heuristic solution: objective 1716.1249532
Found heuristic solution: objective 1716.0259495
Found heuristic solution: objective 1716.0123746
Found heuristic solution: objective 1714.9370344
Found heuristic solution: objective 1714.8304525
Found heuristic solution: objective 1714.7277646
Elapsed time for NoRel heuristic: 38s (best bound 1400)
Elapsed time for NoRel heuristic: 46s (best bound 1400)
Found heuristic solution: objective 1713.4013889
Found heuristic solution: objective 1700.0000000
Elapsed time for NoRel heuristic: 52s (best bound 1400)
Elapsed time for NoRel heuristic: 58s (best bound 1400)
Elapsed time for NoRel heuristic: 67s (best bound 1400)
Elapsed time for NoRel heuristic: 76s (best bound 1400)
Deterministic concurrent LP optimizer: primal and dual simplex
Showing first log only...
Root simplex log...Iteration Objective Primal Inf. Dual Inf. Time
0 1.6040000e+03 2.085614e+01 2.759639e+10 77s
Concurrent spin time: 0.03sSolved with dual simplex
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
965 1.4000000e+03 0.000000e+00 0.000000e+00 77sRoot relaxation: objective 1.400000e+03, 965 iterations, 0.17 seconds (0.19 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time0 0 1400.00000 0 28 1700.00000 1400.00000 17.6% - 78s
H 0 0 1573.2027778 1400.00000 11.0% - 78s
H 0 0 1556.1815818 1400.00000 10.0% - 80s
0 0 1400.00000 0 29 1556.18158 1400.00000 10.0% - 80s
H 0 0 1555.9731799 1400.00000 10.0% - 82s
H 0 0 1555.2870287 1400.00000 10.0% - 82s
H 0 0 1546.8884327 1400.00000 9.50% - 85s
0 0 1400.00000 0 26 1546.88843 1400.00000 9.50% - 85s
H 0 0 1511.4430422 1400.00000 7.37% - 85s
0 0 1400.00000 0 25 1511.44304 1400.00000 7.37% - 85s
0 0 1400.00000 0 37 1511.44304 1400.00000 7.37% - 91s
H 0 0 1511.0000000 1400.00000 7.35% - 92s
0 0 1400.00000 0 36 1511.00000 1400.00000 7.35% - 92s
H 0 0 1510.8966611 1400.00000 7.34% - 97s
0 0 1400.00000 0 42 1510.89666 1400.00000 7.34% - 97s
0 0 1400.00000 0 55 1510.89666 1400.00000 7.34% - 97s
0 0 1400.00000 0 33 1510.89666 1400.00000 7.34% - 100s
0 2 1400.00000 0 31 1510.89666 1400.00000 7.34% - 112s
H 1 4 1510.8960772 1400.00000 7.34% 985 115s
25 34 1400.00000 5 245 1510.89608 1400.00000 7.34% 479 120s
35 44 1400.00000 6 224 1510.89608 1400.00000 7.34% 524 129s
H 37 44 1510.8916234 1400.00000 7.34% 497 129s
H 44 44 1509.7111003 1400.00000 7.27% 422 129s
73 104 1400.00000 9 290 1509.71110 1400.00000 7.27% 285 132s
277 284 1400.00000 29 115 1509.71110 1400.00000 7.27% 194 144s
H 278 284 1509.7092667 1400.00000 7.27% 193 144s
H 279 284 1509.7068256 1400.00000 7.27% 192 144s
287 482 1400.00000 30 41 1509.70683 1400.00000 7.27% 197 151s
H 366 482 1509.7041667 1400.00000 7.27% 176 151s
485 1480 1400.00000 33 131 1509.70417 1400.00000 7.27% 144 159s
1573 2067 1500.00000 61 65 1509.70417 1400.00000 7.27% 72.7 166s
H 2255 2120 1509.7036943 1400.00000 7.27% 72.0 178s
H 2256 2120 1509.6997513 1400.00000 7.27% 72.0 178s
2317 2122 1400.00000 12 34 1509.69975 1400.00000 7.27% 71.5 195s
2318 2123 1500.00000 165 37 1509.69975 1400.00000 7.27% 71.5 213s
2319 2124 1400.00000 81 36 1509.69975 1400.00000 7.27% 71.5 221s
2320 2124 1400.00000 55 38 1509.69975 1400.00000 7.27% 71.4 234sThank you very much.
0 -
The link I sent before is broken due to the fulltop at the end.
The download link again: https://www.filemail.com/d/xsadvtceoialmei
0 -
Your model is quite large and has a lot of nonconvex equalities. Thus, it is expected that the dual bound will proceed only very (very!) slowly.
Your model has the following constraints
qc1: [ - x[0,1,elect] * soc[0] + x[0,1,elect] * soc[1]
+ 29.29444444444444 x[0,1,elect] * b ] = 0Instead of using this quadratic equality constraint, you could instead use an indicator constraint
model.addConstr((x[0,1,elect]==1) >> (- soc[0] + soc[1] + 29.29444444444444 * b == 0))
The indicator constraints means that if \(\texttt{x[0,1,elect]==1}\) then, \(\texttt{- soc[0] + soc[1] + 29.29444444444444 * b == 0}\) has to hold and otherwise not. This works only, because \(\texttt{x[0,1,elect]}\) is a binary variable.
This idea can be applied to quite a lot of your constraints.
Then, you have constraints of the form
qc63101: [ x[0,1,fuel] * soc[0] - x[0,1,fuel] * soc[1] ] = 0
You can again use indicator constraints here to formulate
model.addConstr((x[0,1,fuel]==1) >> (soc[0] - soc[1] == 0))
The same idea applies to the remaining quadratic constraints of the form
qc200741: [ - x[0,340,elect] * soc[0] + x[0,340,elect] * soc[340]
- x[0,340,elect] * d_soc[340] + 60.96111111111111 x[0,340,elect] * b ]
= 00 -
Thank you. I will try that.
So in general, can we say that indicator constraints perform better than bilinear and quadratic constraints?
0 -
So in general, can we say that indicator constraints perform better than bilinear and quadratic constraints?
Usually, yes.
0
Please sign in to leave a comment.
Comments
13 comments