Help on MILP problem. No solution found or slow moving bound! HELP!

[mortenkofoed]

Hi ! 

I am having real problems with the performance of Gurobi. I am modelling a heat storage tank connected to a heat pump. The model is a MILP.  In my approach i am doing a rolling horizon with 1 hour decision time and 9 hour prediction. In the first 16 runs, everything runs without problems or significant cost of time. However, running the 17'th is a nightmare. In the end of the 9 hour prediction horizon a domestic hot water tapping is occurring, hence the heat pump must be turn on in most of the timestep in the 17'th run. I have tried warm start, MIPFocus, Presolve, Cuts, and more nothing seems to work. I am running with 24 cores so computational power should not be a problem.  

See attached mps file. Please let me know if you need more! 

Read MPS format model from file test17.mps

Reading time = 0.01 seconds

DHWHP: 523 rows, 656 columns, 2334 nonzeros

Changed value of parameter Presolve to 2

   Prev: -1  Min: -1  Max: 2  Default: -1

Changed value of parameter PreSparsify to 1

   Prev: -1  Min: -1  Max: 1  Default: -1

Changed value of parameter MIPFocus to 2

   Prev: 0  Min: 0  Max: 3  Default: 0

Read MIP start from file test17_Solable.mst

Optimize a model with 523 rows, 656 columns and 2334 nonzeros

Model has 123 general constraints

Variable types: 615 continuous, 41 integer (41 binary)

Coefficient statistics:

  Matrix range     [1e-04, 3e+03]

  Objective range  [1e+00, 3e+01]

  Bounds range     [1e+00, 1e+02]

  RHS range        [2e-03, 1e+02]

MIP start did not produce a new incumbent solution

Presolve added 113 rows and 17 columns

Presolve time: 0.07s

Presolved: 636 rows, 673 columns, 2581 nonzeros

Variable types: 492 continuous, 181 integer (181 binary)

Presolve removed 105 rows and 105 columns

Presolved: 531 rows, 568 columns, 2369 nonzeros

Extra one simplex iteration after uncrush

Thank you so much for help!!

Best regards 

Morten

  

[Tobias Achterberg]

The model is pretty tough due to the MAX constraints. It is already very hard to find a
feasible solution.

I found a feasible solution by the follwing procedure (which you can do from the Python
interactive shell):

1. Read in the model:
m = read('test17.mps')
2. Convert it into the feas-relax model:
m.feasRelaxS(0, False, False, True)
3. Set the SOS big-M parameters to a smaller value:
m.setParam('PreSOS1BigM', 10)
m.setParam('PreSOS2BigM', 10)
4. Solve the feas-relax model:
m.optimize()
This should yield a solution with zero objective (i.e., without relaxation).
5. Write out this solution:
m.write('test17.sol')

Now, you can use this solution as a MIP start to your original problem. From the command
line, this would be:

gurobi_cl InputFile=test17.sol test17.mps

Maybe, you want to avoid again the SOS translation by setting the SOS big-M parameters:

gurobi_cl PreSOS1BigM=10 PreSOS2BigM=10 InputFile=test17.sol test17.mps


But even with this feasible solution at hand, Gurobi is only able to improve the solution
a little bit, but it is still very hard to increase the global dual bound.


Regards,

Tobias

[Tobias Achterberg]

It seems that your model behaves a little bit better if you use indicator constraints
instead of MAX constraints. Don't know why.

For example, your constraint

Delta_temp_layer_at_index_1_in_between_Layer_2_and_1 = MAX (
Delta_temp_at_index_1_in_between_Layer_2_and_1 , 0 )

can also be modeled as:

Ind_at_index_1_in_between_Layer_2_and_1 = 1 ->
Delta_temp_layer_at_index_1_in_between_Layer_2_and_1 -
Delta_temp_at_index_1_in_between_Layer_2_and_1 = 0

Ind_at_index_1_in_between_Layer_2_and_1 = 0 ->
Delta_temp_layer_at_index_1_in_between_Layer_2_and_1 = 0


But this doesn't make the model easy to solve. I just see the lower bound increasing a bit
faster.


Tobias

[mortenkofoed]

Thank you so much for your time. I guessed it probably could be the reason. I have tried a big M constraint. However, it seems that the MAX constraint worked better overall. I will try to change it to an indicator constraint, however i am not 100% sure how to do this? Could this be the way:

#First variables

 # Create variables for the temperature difference between layers. y
 Delta_temp_layer = {}
 for index in range(0,stop):
      for n in range(1,layers-2):
          Delta_temp_layer [index, n] = model.addVar( vtype=GRB.CONTINUOUS,
                                    name='Delta_temp_layer_at_index_{0}_in_between_Layer_{1}_and_{2}'.format(index, n+1,n))
   
 # Create variables for the temperature difference between layers. y
 Delta_temp = {}
 for index in range(0,stop):
     for n in range(1,layers-2):
        Delta_temp[index, n] = model.addVar(lb=-50, vtype=GRB.CONTINUOUS,
                             name='Delta_temp_at_index_{0}_in_between_Layer_{1}_and_{2}'.format(index, n+1,n))

# Create variable for indicator constraint
New_Binary_variable = model.addVars(range(0,stop), vtype=GRB.BINARY, name="New_Binary_variable")   

#Second Constraint

for i in range(0,stop): 
        for lay in range(1,layers-2):
         

           #Indicator constraint

           model.addConstr(Delta_temp[i,lay] , GRB.EQUAL   , (T[i,lay+1]-T[i,lay]) , name='Delta_temp_{2}_{1}_Constr_1_at_inx{0}'.format(i,lay,lay+1))
           model.addGenConstrIndicator(New_Binary_variable, True, Delta_temp_layer[i,lay] - Delta_temp[i,lay], GRB.EQUAL, 0)

           model.addGenConstrIndicator(New_Binary_variable, False, Delta_temp_layer[i,lay], GRB.EQUAL, 0)

            #MAX Constraint
            #model.addConstr(Delta_temp[i,lay] , GRB.EQUAL   , (T[i,lay+1]-T[i,lay]),name='Delta_temp_{2}_{1}_Constr_1_at_inx{0}'.format(i,lay,lay+1))
            #model.addGenConstrMax(Delta_temp_layer[i,lay], [Delta_temp[i,lay], 0], name="maxconstr_{2}_{1}_Constr_1_at_inx{0}".format(i,lay,lay+1))
       
            ##Big M constraint
            ##model.addConstr(  T[i,lay+1]-T[i,lay], GRB.LESS_EQUAL   , phi[i,lay]*M,             name='Delta_temp_layer_{2}_{1}_Constr_1_at_inx{0}'.format(i,lay,lay+1))
            ##model.addConstr(-(T[i,lay+1]-T[i,lay]),GRB.LESS_EQUAL   , (1-phi[i,lay])*M,         name='Delta_temp_layer_{2}_{1}_Constr_2_at_inx{0}'.format(i,lay,lay+1))
            ##model.addConstr(Delta_temp_layer[i,lay],   GRB.GREATER_EQUAL, (T[i,lay+1]-T[i,lay]),  name='Delta_temp_layer_{2}_{1}_Constr_3_at_inx{0}'.format(i,lay,lay+1))
            ##model.addConstr(Delta_temp_layer[i,lay],   GRB.LESS_EQUAL   , (T[i,lay+1]-T[i,lay])+(1-phi[i,lay])*M,name='Delta_temp_layer_{2}_{1}_Constr_4_at_inx{0}'.format(i,lay,lay+1))
            ##model.addConstr(Delta_temp_layer[i,lay],   GRB.LESS_EQUAL   , phi[i,lay]*M,             name='Delta_temp_layer_{2}_{1}_Constr_5_at_inx{0}'.format(i,lay,lay+1))
       

Best regards Morten

Btw i admire your work in "Constraint Integer Programming" 

[Tobias Achterberg]

Yes, this looks correct. To be on the safe side, just export your model to an lp file and
look at the file in a text editor.

Regards,

Tobias

[mortenkofoed]

When running the script I am receiving following error: 

  File "<ipython-input-96-77af6f8c1532>", line 1, in <module>

    runfile('/Users/mortenkofoed/Dropbox/Thesis/Run2_v15/Dense_Bottom_Multi_layer_v11_MuVar.py', wdir='/Users/mortenkofoed/Dropbox/Thesis/Run2_v15')

  File "/anaconda3/lib/python3.6/site-packages/spyder_kernels/customize/spydercustomize.py", line 668, in runfile

    execfile(filename, namespace)

  File "/anaconda3/lib/python3.6/site-packages/spyder_kernels/customize/spydercustomize.py", line 108, in execfile

    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/Users/mortenkofoed/Dropbox/Thesis/Run2_v15/Dense_Bottom_Multi_layer_v11_MuVar.py", line 400, in <module>

    model.addGenConstrIndicator(New_Binary_variable, True, (Delta_temp_layer[i,lay] - Delta_temp[i,lay]), GRB.EQUAL, 0)

  File "model.pxi", line 3501, in gurobipy.Model.addGenConstrIndicator

AttributeError: 'tupledict' object has no attribute '__colno__' 

Any idea what this mean ? 

Beat regards Morten 

[Tobias Achterberg]

I guess you need to create the "New_Binary_variable" for the cross product of
range(0,stop) and range(1,layers-2), and then pass New_Binary_Variable[i,lay] to the
addGenConstrIndicator() method.

Tobias

[mortenkofoed]

You are off cause right :D. 

However the problem is still difficult to solve.. Do you know if there is any other way to reformulate the conditional constraint so it is not tough? Could you point me in a direction ??. 

[Tobias Achterberg]

It could just be that your model is hard to solve.

But maybe you are lucky and can relax your formulation. For example, if instead of

y = max(x,0)

you use

y >= x
y >= 0

then the model becomes much easier. Of course, this relaxes the problem in such
a way that you could get solutions with y > x and y > 0. Would such solutions be
meaningful in your application? Or would you be able to turn such solutions into
a meaningful solution by some postprocessing step? Maybe you can derive a
heuristic from this approach: first solve the relaxation, then add constraints

y = x

or

y = 0

for all variables in which the solution to your relaxation satisfies such a
condition, turn the remaining conditions back into

y = max(x,0)

and solve this hopefully smaller and simpler model. Of course, as I said, this
would just be a heuristic and would not necessarily provide an optimal solution.
It may even turn out that the smaller model is infeasible.


Regards,

Tobias

[mortenkofoed]

Unfortunately, relaxing the constraint in such way will not give a meningfull solution. As I am modeling a storage tank with the multilayer approach,  y >= x  , y >= 0, will induce temperatures in the storage tank that are zeros in some of the layers.. So this solution is no good. 

I guess i can increase some of the timesteps even further and thereby making the model smaller.. However also getting a less optimal solution.. Any way thank you so much for your help!!.