Solution time drastically increased when in a loop?
Answered
I have a MILP model, in base scenario it can be solved within 10 seconds, see log:
Gurobi 10.0.1 (mac64[arm]) logging started Mon Sep 18 17:43:26 2023
Set parameter LogFile to value "base_scenario/results/mip_log.txt"
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (mac64[arm])
CPU model: Apple M1
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 1694 rows, 6293 columns and 24772 nonzeros
Model fingerprint: 0x92539075
Model has 8978 general constraints
Variable types: 1 continuous, 6292 integer (408 binary)
Coefficient statistics:
Matrix range [1e-02, 3e+05]
Objective range [9e-02, 2e+05]
Bounds range [1e+00, 2e+07]
RHS range [4e+03, 2e+04]
GenCon rhs range [1e+00, 1e+00]
GenCon coe range [1e+00, 1e+00]
Presolve removed 44 rows and 373 columns
Presolve time: 0.07s
Presolved: 1650 rows, 5920 columns, 23713 nonzeros
Presolved model has 4 SOS constraint(s)
Variable types: 0 continuous, 5920 integer (34 binary)
Root relaxation: objective 7.937310e+07, 6640 iterations, 0.08 seconds (0.11 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 7.9373e+07 0 624 - 7.9373e+07 - - 0s
0 0 8.0481e+07 0 650 - 8.0481e+07 - - 0s
0 0 8.0482e+07 0 691 - 8.0482e+07 - - 0s
H 0 0 8.451319e+07 8.0482e+07 4.77% - 0s
0 0 8.0486e+07 0 690 8.4513e+07 8.0486e+07 4.77% - 0s
0 0 8.0504e+07 0 807 8.4513e+07 8.0504e+07 4.74% - 0s
H 0 0 8.450392e+07 8.0504e+07 4.73% - 0s
0 0 8.0555e+07 0 1006 8.4504e+07 8.0555e+07 4.67% - 0s
0 0 8.0555e+07 0 1009 8.4504e+07 8.0555e+07 4.67% - 0s
0 0 8.0555e+07 0 1009 8.4504e+07 8.0555e+07 4.67% - 0s
0 0 8.0558e+07 0 1009 8.4504e+07 8.0558e+07 4.67% - 0s
H 0 0 8.439609e+07 8.0581e+07 4.52% - 0s
H 0 0 8.435109e+07 8.0581e+07 4.47% - 0s
0 0 8.0595e+07 0 813 8.4351e+07 8.0595e+07 4.45% - 0s
0 0 8.0607e+07 0 807 8.4351e+07 8.0607e+07 4.44% - 0s
H 0 0 8.385856e+07 8.0607e+07 3.88% - 0s
0 2 8.0607e+07 0 799 8.3859e+07 8.0607e+07 3.88% - 0s
H 33 40 8.333525e+07 8.0929e+07 2.89% 654 1s
H 72 64 8.333525e+07 8.0998e+07 2.80% 441 1s
H 988 313 8.333525e+07 8.2107e+07 1.47% 234 3s
H 1528 388 8.333525e+07 8.2287e+07 1.26% 225 3s
2244 334 8.3238e+07 14 597 8.3335e+07 8.2573e+07 0.91% 245 5s
H 2826 128 8.333525e+07 8.3288e+07 0.06% 265 6s
Cutting planes:
Gomory: 6
Implied bound: 1422
MIR: 98
Inf proof: 13
Relax-and-lift: 92
Explored 3039 nodes (763146 simplex iterations) in 6.14 seconds (10.37 work units)
Thread count was 8 (of 8 available processors)
Solution count 10: 8.33353e+07 8.33353e+07 8.33353e+07 ... 8.45132e+07
Optimal solution found (tolerance 1.00e-04)
Best objective 8.333525193000e+07, best bound 8.333524871943e+07, gap 0.0000%
However, I want to perform a sensitivity analysis on the model by changing one of the parameter. I implemented it as follows:
percentage_variations = [-0.5, -0.25, 0.0, 0.25, 0.50]
...
for change in percentage_variations:
new_p = 1 + change
....
for f in capacity.filter(like='f').columns:
value = capacity.loc[t, f]
CAP_f[(t, f)] = value * new_p
....
m.optimize()
But when it comes to the loop where change = 0.0, the model takes too long too solve, exceeding the time limit set for 20 seconds, see log below:
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (mac64[arm])
CPU model: Apple M1
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 8470 rows, 31465 columns and 123860 nonzeros
Model fingerprint: 0xaba2dfd1
Model has 44890 general constraints
Variable types: 5 continuous, 31460 integer (2040 binary)
Coefficient statistics:
Matrix range [1e-02, 3e+05]
Objective range [9e-02, 2e+05]
Bounds range [1e+00, 2e+07]
RHS range [4e+03, 2e+04]
GenCon rhs range [1e+00, 1e+00]
GenCon coe range [1e+00, 1e+00]
Presolve removed 6820 rows and 25545 columns
Presolve time: 0.74s
Presolved: 1650 rows, 5920 columns, 23713 nonzeros
Presolved model has 4 SOS constraint(s)
Variable types: 0 continuous, 5920 integer (34 binary)
Root relaxation: objective 7.885152e+07, 6279 iterations, 0.07 seconds (0.09 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 7.8852e+07 0 568 - 7.8852e+07 - - 0s
0 0 7.8961e+07 0 674 - 7.8961e+07 - - 0s
0 0 7.9955e+07 0 546 - 7.9955e+07 - - 0s
0 0 7.9955e+07 0 554 - 7.9955e+07 - - 0s
H 0 0 8.508952e+07 7.9966e+07 6.02% - 1s
0 0 7.9974e+07 0 927 8.5090e+07 7.9974e+07 6.01% - 1s
0 0 8.0000e+07 0 1033 8.5090e+07 8.0000e+07 5.98% - 1s
0 0 8.0002e+07 0 1048 8.5090e+07 8.0002e+07 5.98% - 1s
0 0 8.0002e+07 0 1044 8.5090e+07 8.0002e+07 5.98% - 1s
0 0 8.0037e+07 0 1018 8.5090e+07 8.0037e+07 5.94% - 1s
0 0 8.0044e+07 0 1018 8.5090e+07 8.0044e+07 5.93% - 1s
H 0 0 8.491365e+07 8.0044e+07 5.73% - 1s
0 2 8.0044e+07 0 1005 8.4914e+07 8.0044e+07 5.73% - 1s
H 35 40 8.488360e+07 8.0281e+07 5.42% 634 2s
H 36 40 8.475934e+07 8.0281e+07 5.28% 626 2s
H 71 70 8.467602e+07 8.0281e+07 5.19% 486 2s
H 75 70 8.466732e+07 8.0281e+07 5.18% 472 2s
H 78 77 8.463433e+07 8.0281e+07 5.14% 463 2s
H 79 77 8.455571e+07 8.0281e+07 5.06% 461 2s
H 295 237 8.355880e+07 8.0350e+07 3.84% 238 2s
H 300 237 8.355720e+07 8.0350e+07 3.84% 235 2s
H 323 237 8.336281e+07 8.0350e+07 3.61% 229 2s
H 612 388 8.336089e+07 8.0428e+07 3.52% 217 3s
H 615 388 8.335312e+07 8.0428e+07 3.51% 220 3s
H 619 388 8.332382e+07 8.0428e+07 3.48% 219 3s
H 1296 487 8.324190e+07 8.0969e+07 2.73% 211 4s
1347 544 8.1493e+07 12 520 8.3242e+07 8.1001e+07 2.69% 213 5s
H 2242 956 8.324189e+07 8.1187e+07 2.47% 179 6s
H 2318 967 8.324188e+07 8.1187e+07 2.47% 183 7s
H 2321 921 8.302042e+07 8.1187e+07 2.21% 183 7s
H 2349 882 8.296757e+07 8.1187e+07 2.15% 182 7s
H 2350 840 8.295444e+07 8.1187e+07 2.13% 182 7s
H 2352 801 8.295090e+07 8.1187e+07 2.13% 182 7s
H 3370 689 8.294013e+07 8.1187e+07 2.11% 148 8s
H 5009 747 8.293935e+07 8.1187e+07 2.11% 134 9s
H 5010 747 8.293909e+07 8.1187e+07 2.11% 134 9s
H 5018 747 8.293818e+07 8.1187e+07 2.11% 134 9s
5072 776 8.2898e+07 28 391 8.2938e+07 8.1187e+07 2.11% 134 10s
H 8139 1380 8.293621e+07 8.1510e+07 1.72% 115 11s
H 9921 1606 8.283271e+07 8.1514e+07 1.59% 105 12s
15476 3132 8.2828e+07 1087 65 8.2833e+07 8.1922e+07 1.10% 95.6 15s
H25193 437 8.269771e+07 8.1922e+07 0.94% 79.5 18s
H26042 184 8.269768e+07 8.2632e+07 0.08% 79.7 18s
30282 3091 8.2688e+07 660 25 8.2698e+07 8.2683e+07 0.02% 70.4 20s
Explored 30309 nodes (2140812 simplex iterations) in 20.01 seconds (29.34 work units)
Thread count was 8 (of 8 available processors)
Solution count 10: 8.26977e+07 8.26977e+07 8.28327e+07 ... 8.29544e+07
Time limit reached
Best objective 8.269768433600e+07, best bound 8.268332877545e+07, gap 0.0174%
I also tried m.reset(0) as an alternative, but it did not change this problem.
But if I change the order of percentage_variations so that 0.0 becomes the first item in the list, the model gets solved in 7 seconds.
I would be thankful if anyone can explain what the reason for this is.
-
Do I understand correctly that both logs should belong to the same model? This is not the case. The second log shows a different objective
Best objective 8.269768433600e+07, best bound 8.268332877545e+07, gap 0.0174%
vs
Best objective 8.333525193000e+07, best bound 8.333524871943e+07, gap 0.0000%
Also, while the statistics for the presolved models look equivalent, the initial model in the second log contains more rows, columns, and variables. If I understand correctly the models should only differ in some parameter but not in the initial size.
Do you know the Multiple Scenarios feature of Gurobi that is designed for doing some sensitive analysis? The advantage of this feature is that all variations of one model are solved in one run. See also the Gurobi Python example multiscenario.py
0 -
Thank you for your reply.
The second model is built based on the first model, with the only change being adding the percentage_variations and the for-loop. The input parameters are the same. I thought for the iteration where percentage_variations is 0, there objective would be the same, since there is no change in input parameter values and constraints etc.
Regarding multiscenario.py, if I have the following parameters and constraints:
a = [0.36, 0.27, 0.44, 0.46, 0.42]
b = [0.19, 0.28, 0.23, 0.27, 0.2]
c = [0.23, 0.18, 0.29, 0.14, 0.27]
x = model.addVar(vtype=gp.GRB.CONTINUOUS, name="x")
y = model.addVar(vtype=gp.GRB.CONTINUOUS, name="y")
z = model.addVar(vtype=gp.GRB.CONTINUOUS, name="z")
for t in len(a):
model.addConstr(x + y <= 5 * a[t] name="constraint1")
model.addConstr(2 * x - y <= 8 * b[t], name="constraint2")
model.addConstr(3 * x + 2 * y >= 10 * c[t], name="constraint2")
...In one scenario, the range of values of a should be changed to random.uniform(0.5,0.8), for b it should be random.uniform(0.2,0.4) and c random.uniform(0.1,0.2). Could this also be implemented using multiscenario?
0 -
As I understand it you want to change the RHS of constraints1 in a first scenario, the RHS of constraints2 in a second scenario, and the RHS of constraints3 in a third scenario.
This could be done as followsconstr1 = model.addConstrs((x + y <= 5 * a[t] for t in range(len(a))), name="constraint1")
constr2 = model.addConstrs((2 * x - y <= 8 * b[t] for t in range(len(b))), name="constraint2")
constr3 = model.addConstrs((3 * x + 2 * y >= 10 * c[t] for t in range(len(c))), name="constraint3")
model.NumScenarios = 4
# Scenario 0: Base model, hence, nothing to do except giving the scenario a name
model.Params.ScenarioNumber = 0
model.ScenNName = 'Base model'
# Scenario 1: change RHS of constr1 to be random numbers between 0.5 and 0.8
model.Params.ScenarioNumber = 1
model.ScenNName = 'change constr1'
for t in range(len(a)):
constr1[t].ScenNRhs = np.random.uniform(0.5,0.8)
# Scenario 2: change RHS of constr2
model.Params.ScenarioNumber = 2
model.ScenNName = 'change constr2'
for t in range(len(b)):
constr2[t].ScenNRhs = np.random.uniform(0.2,0.4)
# Scenario 3: change RHS of constr3
model.Params.ScenarioNumber = 3
model.ScenNName = 'change constr3'
for t in range(len(c)):
constr3[t].ScenNRhs = np.random.uniform(0.1,0.2)For general debugging, I would recommend writing your model into an LP-file using Model.write() before optimizing to check if this is indeed the model you intend to solve.0 -
Thank you so much for your help.
Can you explain briefly what ScenNRhs does? I am not certain how the ScenNRhs behaves. The documentation says it changes the RHS of the linear constraint. But if that is the case, wouldn't constr1[t].ScenNRhs = np.random.uniform(0.5,0.8) make it so that for scenario 1, constraint 1 becomes constr1 = model.addConstrs((x + y <= np.random.uniform(0.5,0.8)) ?
If that is the case, how should a constraint with multiple variables on both sides be modified for different scenarios? For exampleconstr1 = model.addConstrs((x + y == z * a[t] for t in range(len(a))), name="constraint_1")
0 -
To better understand how the attribute ScenNRHS can be used, you need to understand what the right-hand side (RHS) of a constraint is. Consider the constraint \(x + y - 5z = 2\). It can equivalently be written e.g. as \(x + y - 2 = 5z\). But both represent the same constraint and the RHS is 2. To identify the RHS, you need to transform the constraint s.t. all variables are on the left and all constants are on the right (which is the case for the first expression).
(You could also check the attribute RHS.)
If you have constraint \(x + y - 5z = 2\) and want to consider the alternative constraint \(x + y - 5z = 3\) in a multiple-scenario setting, you can do this by changing the attribute ScenNRHS from 2 to 3 of this constraint.Let's consider your type of example: Constraint \(x + y = 2z\) can be written as \(x + y - 2z = 0\). The RHS is 0. Now you do not want to change the RHS but the coefficient of variable \(z\). This is not as easy as changing the RHS and needs a trick. Assume in a second scenario you want to consider the constraint \(x + y = 3z\). This can be done as follows:
# reformulate == constraint as two <= and >= and constraints
constrBaseLT = model.addConstr(x+y-2*z <= 0)
constrBaseGT = model.addConstr(x+y-2*z >= 0)
# add the constraints for the different scenario from the beginning, but de-activate them
# de-activate here means that the constraints are always satisfied
constrScenLT = model.addConstr(x+y-3*z <= gp.GRB.INFINITY)
constrScenGT = model.addConstr(x+y-3*z >= -gp.GRB.INFINITY)
model.NumScenarios = 2
# Scenario 0: Base model, hence, nothing to do except giving the scenario a name
model.Params.ScenarioNumber = 0
model.ScenNName = 'Base model'
# Scenario 1: activate the alternative constraints and de-activate the original ones
constrBaseLT.ScenNRHS = gp.GRB.INFINITY
constrBaseGT.ScenNRHS = -gp.GRB.INFINITY
constrScenLT.ScenNRHS = 0
constrScenGT.ScenNRHS = 0Please also have a look at Tips and Tricks for multiple scenarios.
0
Please sign in to leave a comment.
Comments
5 comments