Skip to main content

Computing problem arising in loop

Answered
Permanently deleted user
  • Edited
First Comment
First Question 
m = gp.Model('Mixed Integer Program')
m.setParam('Outputflag',0)
m.setParam('Timelimit' , 600)
X = m.addMVar(p.shape, vtype="B", ub=~ pd.isnull(p))
S = m.addMVar(p.shape, vtype="B", ub=~ pd.isnull(p))
p.fillna(1, inplace = True)
samples = []
fields = []
SS = np.arange(0, 1650, 25)
m.addConstrs(X[y,x]  >= S[y-h,x-w] for h in range(size[1]) for w in range(size[0]) for y in y_range for x in x_range)
m.addConstrs(X[y,x]  <= gp.quicksum(gp.quicksum(S[y-h,x-w] for h in range(size[1])) for w in range(size[0])) for y in y_range for x in x_range)
m.addConstr(gp.quicksum(S[y,x] for y in H for x in W) <= 0)
m.update()
objective = gp.quicksum(gp.quicksum((1-X[y,x])*(1-p[x][y]) for y in H) for x in W)
m.setObjective(objective, gp.GRB.MINIMIZE)

for area in SS:
    
    m.setParam("Timelimit", 180)
    m.remove(m.getConstrs()[-1])            
    m.addConstr(gp.quicksum(S[y,x] for y in H for x in W) <= area)
    m.update()
    m.optimize()
    obj = m.getObjective()
    d = obj.getValue()
    sample = X.X
    field = S.X           
    print(Fields: ', np.sum(field))
    print('Sample size: ', np.sum(sample))
    print('Solution: ', round(d,2))
    print('Final MIP gap value: ', m.MIPGap)
    print('--------------------------------------------------------')
    
    samples.append(sample)
    fields.append(field)


Here:
p is a dataframe consisting of nan-entries and probabilities.

H = range(p.shape[0])
W= range(p.shape[1])
x_range = range(size[0]-1, p.shape[1])
y_range = range(size[1]-1, p.shape[0])

For my current problem I have:
p has size (207,130), and has 21156 numeric entries between 0 and 1 and 5754 nan entries

size = [4,4]

Adding all constraints takes almost 6 minutes which is why we set a 10-minute timelimit for the initialization before the loop. Then in the loop we alter the last constraint by removing and adding it for the new value of area and use a timelimit of 3 minutes. The solver finds the optimum solution (MIPgap = 0) up to area = 325, and finds really good solutions (MIPgap < 0.035) within 3 minutes up to area = 550. For area = 550 finds an improved solution with a MIPgap of 0.020, but then all areas up to 1425 do not find an improvement anymore and use the solution of 550. From 1450 it finds the optimal solution again, for which the objective equals 0.  I can't figure out why this is the case, since if an S[y,x] becomes 1, up to 16 (4*4) X[y,x]'s become 1 which can only lead to a decrease in the objective value. There are still a lot of S[y,x] = 0 that force X[y,x]'s whit p[y,x] != 1 to become 1, so there is still much room for an improved solution. 

On top of that, when we change the <= sign in the last constraint to an = sign, we see that from 575 to 1425 it finds solutions with a worse objective value than for area = 550. I can't figure out why the problem becomes so much more complicated from this point for the solver, while the problem itself is fairly simple. Also, I removed the timelimit and let it run till it finds the optimal solution, which it did quite fast for area values up to 500, but become so much slower afterwards. For 575 it couldn't find the exact solution in 12 hours after which I interrupted the kernel. 

What is causing these complications and what can be done to improve the solver?

0

Comments

5 comments

Date Votes
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi Jelle,

    I won't aim to answer everything in this first response because I think it would be good to address some things first.

     

    X = m.addMVar(p.shape, vtype="B", ub=~ pd.isnull(p))  # Matrix API
S = m.addMVar(p.shape, vtype="B", ub=~ pd.isnull(p))  # Matrix API
    ...
    # term based modelling
    m.addConstrs(X[y,x]  >= S[y-h,x-w] for h in range(size[1]) for w in range(size[0]) for y in y_range for x in x_range)

    Don't mix Matrix API and term based modelling unless there are exceptional circumstances, it will result in slow model building time.  Use addVars method instead of addMVar in this case to avoid this.  Note you cannot pass a 2D pandas dataframe the ub argument of addVars, but the following should work:

    ub=(~pd.isnull(p)).to_numpy().flatten()

     

    Adding all constraints takes almost 6 minutes which is why we set a 10-minute timelimit for the initialization before the loop

    The timelimit does not work the way you are thinking it does.  It limits the time the solver may spend in the optimization routine, i.e. after you call m.optimize(), and has no relation to the time spent building the model.  (Similarly the Runtime attribute will only pertain to the time spent in optimization, not model build time).  The only time limit that takes effect is this one: 

    m.setParam("Timelimit", 180)

    and you only need to set it once (so you can set it outside of the loop).

     

    m.addConstr(gp.quicksum(S[y,x] for y in H for x in W) <= 0)
...

for area in SS:
    ...
    m.remove(m.getConstrs()[-1])            
    m.addConstr(gp.quicksum(S[y,x] for y in H for x in W) <= area)

    This is overkill for what you are wanting to achieve.  Constraint objects have a RHS attribute which can be adjusted:

    con = m.addConstr(gp.quicksum(S[y,x] for y in H for x in W) <= 0)
    ...

    for area in SS:
        ...
        con.RHS = area

    When you do this you no longer have to use m.getConstrs() and as a result can remove all instances of update() in your code (which slows it down).  You will also benefit from warmstarting from the previous solve.

    I thoroughly recommend you read the Python section of our manual, especially Attributes and Parameters, it will be a very good investment of your time if you plan to use Gurobi in any serious way.

    Once the above is addressed, let's return to this and see if the performance issue still persists?

    - Riley

    1
  • Permanently deleted user
    First Comment
    First Question

    Riley Clement

    Hi Riley,

    Thank you for your comments. By making the adjustment to the script, the program is indeed much faster in adding the constraints to the model. However, the performance issues are still there. For a timelimit of 300 seconds, we see that for area = 875 onwards the objective value sees a sudden increase again together with an illogical solution. For example, setting 25 random values of S that are 0 in the solution of the previous iteration equal to 1 would results in a lower objective value, but this is not achieved by the solver. This is probably primarily caused by the timelimit, but even when I set it at 900 seconds, the optimizer will struggle from area = 1150 onwards, finding the same illogical type of solution as before. Are there any aspects that I can change to the settings of the optimizer or my code that would allow the optimizer to find a (near) optimal solution faster for my specific problem? It finds the optimal values for smaller values of area very fast, but struggles when these become larger. You can find the current code below. On top of that, a plot of the objective value over the values of area is added.Thank you for your help!

    m = gp.Model('Mixed Integer Program')
m.setParam('Outputflag',0)
m.setParam('Timelimit' , 300) 
    X = m.addVars(p.shape[0], p.shape[1], vtype="B", ub = (~pd.isnull(p)).to_numpy().flatten())
    S = m.addVars(prob.shape[0], prob.shape[1], vtype="B", ub = (~pd.isnull(p)).to_numpy().flatten())
    p.fillna(1, inplace = True)
samples = []
fields = []
    objectives = []
SS = np.arange(0, 1650, 25)
    
m.addConstrs(X[y,x]  >= S[y-h,x-w] for h in range(size[1]) for w in range(size[0]) for y in y_range for x in x_range)
m.addConstrs(X[y,x]  <= gp.quicksum(gp.quicksum(S[y-h,x-w] for h in range(size[1])) for w in range(size[0])) for y in y_range for x in x_range)
cor = m.addConstr(gp.quicksum(S[y,x] for y in H for x in W) == 0)

objective = gp.quicksum(gp.quicksum((1-X[y,x])*(1-p[x][y]) for y in H) for x in W)
m.setObjective(objective, gp.GRB.MINIMIZE)

for area in SS:          
        cor.RHS = area
        
    m.optimize()
    obj = m.getObjective()
    d = obj.getValue()
    sample = np.array([[X[i, j].X for j in range(prob.shape[1])] for i in range(prob.shape[0])])     
        field = np.array([[S[i, j].X for j in range(prob.shape[1])] for i in range(prob.shape[0])])               
        
        print(Fields: ', np.sum(field))
    print('Sample size: ', np.sum(sample))
    print('Solution: ', round(d,2))
    print('Final MIP gap value: ', m.MIPGap)
    print('--------------------------------------------------------')
    
    samples.append(sample)
    fields.append(field)
        objectives.append(d)
    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi Jelle,

    There's a few tweaks I'd recommend for your formulation.  Noting that the objective function wants to increase the value of X variables, this family of constraints will always be satisfied by good solutions:

    X[y,x]  >= S[y-h,x-w]

    If you remove them then linear relaxations can be solved faster.

    I also note that the "cor" constraint was switched to an equality in your latest iteration of the code:

    gp.quicksum(S[y,x] for y in H for x in W) == 0

    I would change this back to "<=" as it allows the solution from one solve to be used in the next, which will help with pruning the B&B tree, and be especially useful for the harder solves towards the end of the loop.  It will also be impossible to get the huge jumps in the objective function that you observed in the graph.

    As you increase the RHS for the Cor constraint the problem becomes harder.  This seems intuitive, for small RHS values the 4x4 areas of the grid activated by a particular S variable would not likely overlap, and finding the best "placement of these 4x4 areas" seems relatively trivial.  For higher RHS values of the Cor constraint these areas overlap and the problem gets harder.

    Note that finding feasible solutions to your problem is easy, and the dual bound looks relatively strong.  The difficulty for the later iterations is in finding great solutions.  I have noticed that our No Relaxation ("NoRel") heuristic works well for finding good solutions for these harder problems.  By default it does not run but you can use NoRelHeurTime and NoRelHeurWork to change this.  I've also observed that Barrier works best for the root LP so setting Method=2 would be a good idea.

    - Riley

     

    0
  • Permanently deleted user
    First Comment
    First Question

    Riley Clement

    Hi Riley,

    Again, much thanks for all your help! These recommendations help a lot, I already implemented the Method=2, changed the equal sign and removed the set of constraints you mentioned and will have to read more into the NoRel parameters to determine which I should use and what parameter value to set. I saw that earlier you asked for the input data to run the model. Unfortunately, I can not share the data with you, but I can give you an idea of it with the following two pictures. In the first picture, the blue area displays the locations in the dataframe of non-nan entries, and the grey are contains all the nan-entries. The non-nan entries form a circular shape together. The second picture displays a histogram of the values of p. Most of the lower values are located towards the boundaries of the circle and the higher values more towards the center, which will indeed lead to overlapping areas for higher values of the RHS of Cor. From your response, I can already tell that you have a good understanding of my problem and thank you for all your help.

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Thanks Jelle, I ended up using randomly generated data but I expect the same principles apply.  Obviously a circle is a very symmetrical object, if the data values are also symmetrical then setting Symmetry=2 might help. 

    - Riley

    0

Please sign in to leave a comment.

Didn't find what you were looking for?

New post