Skip to main content

Gurobi seems not able to solve an SOCP problem to optimal

Answered

Comments

8 comments

  • Simranjit Kaur
    Gurobi Staff Gurobi Staff

    Hi Fei Wang,

    The solution

    x = array([[0.5, 1. ],
          [0.5, 1. ]])
    y = array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])

    is not feasible for the LP model you provided unless the \(u\) variables can take negative values.

    For instance, if \(x[0,0] = 0.5\), then the constraint \(x[0,0]+u[0,0] = 0\), implies \(u[0,0]\) must be \(-0.5\). However, since the \(u\) variables have no negative lower bounds specified in the LP model, they are considered non-negative by default.

    Can you please verify if you are solving the same model and obtaining different optimal solutions?

    Regards,

    Simran

     

    1
  • Fei Wang
    Conversationalist
    First Question

    Hi Simran,

    Thank you very much! It solved my problem.  I didn't know the default lower bound is zero, 

    1
  • Fei Wang
    Conversationalist
    First Question

    I have another question:

    I found that adding the constraints to the model is the most expensive part of the code
    So it is this command  model.addConstr()   that is costing 90% of the time

    Since I am just adding some linear constraints and SOCP constraints, why adding them to the model cost so much time?

    0
  • Simranjit Kaur
    Gurobi Staff Gurobi Staff

    Hi Fei Wang,

    Can you please share an example code showing how you build your model?

    The article How do I improve the time to build the model may be helpful for you. It contains several tips to speed up model building time.

    Regards,

    Simran

    0
  • Fei Wang
    Conversationalist
    First Question

    Hi Simran,

    Here is part of the code that’s very expensive. Most of the time is spent on adding the constraints especially the first and second for loops

    model = gp.Model()
    model.Params.LogToConsole = 0
    x = model.addMVar((len(self.steiner_points), self.r), lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS, name="x")
    d = model.addVars(len(E), lb=0, vtype=GRB.CONTINUOUS, name="d")
    de = model.addVars(len(E), lb=0, vtype=GRB.CONTINUOUS, name="de")
    c = model.addMVar(len(self.ET), lb=0, vtype=GRB.CONTINUOUS, name="c")
    y = model.addMVar(len(self.ET), lb=0, ub=1, vtype=GRB.CONTINUOUS, name="y")
    z = model.addVars(len(self.ET), lb=0, vtype=GRB.CONTINUOUS, name="z")
    # ze = model.addVars(len(self.ET), lb=0, vtype=GRB.CONTINUOUS, name="z")
    w = model.addMVar((len(E), self.r), lb=-GRB.INFINITY,vtype=GRB.CONTINUOUS, name="w")
    u = model.addMVar((len(self.ET), self.r), lb=-GRB.INFINITY,vtype=GRB.CONTINUOUS, name="u")

    obj = gp.quicksum(d[i] for i in range(len(E))) + gp.quicksum(c[i] for i in range(len(self.ET)))
    model.setObjective(obj, GRB.MINIMIZE)

    start_time = time.time()
    for (i, j) in E:
    if i in self.terminal_indices_label_added_dict:
    label_terminal_i = self.steiner_tree.get_label(i)
    for k in range(self.r):
    model.addConstr(
    w[E_index_dict[(i, j)], k] == self.terminal_label_coord_dict[label_terminal_i][k]
    - x[Steiner_points_index_dict[j], k])
    model.addConstr(de[E_index_dict[(i, j)]] == gp.norm(w[E_index_dict[(i, j)], :], 2))
    model.addConstr(d[E_index_dict[(i, j)]] >= de[E_index_dict[(i, j)]])
    else:
    for k in range(self.r):
    model.addConstr(w[E_index_dict[(i, j)], k] == x[Steiner_points_index_dict[i], k] - x[
    Steiner_points_index_dict[j], k])

    model.addConstr(de[E_index_dict[(i, j)]] == gp.norm(w[E_index_dict[(i, j)], :], 2))
    model.addConstr(d[E_index_dict[(i, j)]] >= de[E_index_dict[(i, j)]])

    for (i, j) in self.ET:
    label_terminal_i = self.terminals_index_label_not_added_dict[i]
    M_i = self.M[label_terminal_i]
    # model.addConstr(z[ET_index_dict[(i, j)]] == c[ET_index_dict[(i, j)]] + M_i * (1 - y[ET_index_dict[(i, j)]]))
    for k in range(self.r):
    model.addConstr(u[ET_index_dict[(i, j)], k] == self.terminal_label_coord_dict[label_terminal_i][k]
    - x[Steiner_points_index_dict[j], k])
    model.addConstr(z[ET_index_dict[(i, j)]] ==
    gp.norm(u[ET_index_dict[(i, j)], :], 2))
    model.addConstr(c[ET_index_dict[(i, j)]] >=
    z[ET_index_dict[(i, j)]] - M_i * (1 - y[ET_index_dict[(i, j)]]))

    for i in self.terminal_indices_not_added:
    model.addConstr(gp.quicksum(y[ET_index_dict[(i, j)]] for j in self.steiner_points_degree_less_3) == 1)

    for j in self.steiner_points_degree_less_3:
    degree_j = len(self.steiner_tree.nodes[j]['neighbors'])
    model.addConstr(gp.quicksum(y[ET_index_dict[(i, j)]] for i in self.terminal_indices_not_added)
    == 3 - degree_j)

    time_solve += time.time() - start_time
    model.optimize()
    0
  • Simranjit Kaur
    Gurobi Staff Gurobi Staff

    Hi Fei Wang,

    The majority of the time spent in adding constraints to your model is perhaps consumed in the for loops like the one below:

    for k in range(self.r):
       model.addConstr( w[E_index_dict[(i, j)], k] == self.terminal_label_coord_dict[label_terminal_i][k]- x[Steiner_points_index_dict[j], k])

    Instead of looping over k, you can exploit the MVar nature of the variables and replace the above with the equivalent following constraint:

    model.addConstr( w[E_index_dict[(i, j)]] == terminal_label_coord_mat[label_terminal_i] - x[Steiner_points_index_dict[j]])
    w[E_index_dict[(i, j)] and x[Steiner_points_index_dict[j] are k dimensional arrays of variables, and this will add one constraint for each value of k.
     
    Please note that for the above constraint to work, you would need to convert the dictionary terminal_label_coord_dict to a 2D matrix terminal_label_coord_mat appropriately.
     
    Could you please make the equivalent changes for all for loops and test if this helps?
     
    Regards,
    Simran
     
     
    0
  • Fei Wang
    Conversationalist
    First Question

    Thanks Simran,

    It works. However in my code self.r is only 2. 

    So this changes won't save a lot of time, in fact, the majority of the time is spent in adding those SOCP constraints.

    Constraints as the following are taking a lot of time

     model.addConstr(de[E_index_dict[(i, j)]] == gp.norm(w[E_index_dict[(i, j)], :], 2))
    1
  • Simranjit Kaur
    Gurobi Staff Gurobi Staff

    Hi Fei Wang,

    Indeed, my previous suggestion will help save time for larger values of self.r.

    With respect to the SOCP constraints, you can try adding them using the addGenConstrNorm method. It is slightly faster than using gp.norm in my test below:

    import gurobipy as gp
    import time as time
    n = 100000
    time_solve = 0

    model = gp.Model()
    de = model.addMVar(n, name = "de")
    w = model.addMVar((n, 2), name = "w")

    # test gp.norm
    start_time = time.time()
    for i in range(n):
      s = model.addConstr(de[i] == gp.norm(w[i,:], 2))
    time_solve = time.time() - start_time
    print("time_solve_1", time_solve)  

    # test GenConstrNorm
    start_time = time.time()
    for i in range(n):
       m = model.addGenConstrNorm(de[i], w[i], 2)
    time_solve = time.time() - start_time
    print("time_solve_2", time_solve)

    Regards,

    Simran

    0

Please sign in to leave a comment.