Facility Location Problem - Standard BIP vs QUBO

Awaiting user input

Abhishek Solomon

• March 04, 2023 00:10

Hello,

I'm fairly new to the optimization world, and working on a research project that involves comparison of solving times for the standard uncapacitated facility location problem, using the following techniques:

1) As a standard BIP:

    TotalFixedCost = gp.quicksum(fixedcost[i]*FacilitySelect[i] for i in facilities)
    TotalServiceCost = gp.quicksum(servicecost[i,j]*ServiceSelect[i,j] for i in facilities for j in customers)
    TotalCost = TotalFixedCost + TotalServiceCost
    m.setObjective(TotalCost , GRB.MINIMIZE)

2) As a quadratic unconstrained binary optimization (QUBO) model with a penalty function:

 TotalFixedCost = gp.quicksum(fixedcost[i]*FacilitySelect[i] for i in facilities)
 TotalServiceCost = gp.quicksum(servicecost[i,j]*ServiceSelect[i,j] for i in facilities for j in customers)
 Cost = TotalFixedCost + TotalServiceCost
 Penalty1 = p*(gp.quicksum((1-gp.quicksum(ServiceSelect[i,j] for i in facilities))*(1-gp.quicksum(ServiceSelect[i,j] for i in facilities)) for j in customers))
 Penalty2 = p*(gp.quicksum(ServiceSelect[i,j] - ServiceSelect[i,j]*FacilitySelect[i] for i in facilities for j in customers))
 TotalCost = Cost + Penalty1 + Penalty2
 m.setObjective(TotalCost , GRB.MINIMIZE)

3) As a quadratic unconstrained binary optimization (QUBO) model, solving the xTQx objective, (with the Q matrix given as input):

    x = m.addVars(Q.shape[0], vtype=gp.GRB.BINARY)
    objective = gp.quicksum(x[i] * x[j] * Q[i, j] for i in range(Q.shape[0]) for j in range(Q.shape[0]))
    m.setObjective(objective, gp.GRB.MINIMIZE)

Using the standard BIP formulation, Gurobi is able to solve my sample problem sets instantly with 0 gap. However, the 2 QUBO methods take very long times for the same sets, and have large gaps (especially the penalty function method - in the 1000s).

Is it normal to have large optimality gaps while using the QUBO solver? Do i need to tweak any setting to tune the solver?

0

• 

  Jaromił Najman

  Gurobi Staff

  • March 07, 2023 12:56

  Hi Abishek,

  In your BIP formulation, you don't have any constraints. There is only a linear objective function, thus the model is easy to solve.

  In your QUBO with penalty function, you introduce 2 constraints via penalties. The constraints are enforced because the penalty is bigger if they are not fulfilled.

  The first constraint is in \(\texttt{Penalty1}\) and reads

  \[\begin{align*}
  \sum_{i} \text{ServiceSelect}_{i,j} = 1 \quad \forall j
  \end{align*}\]

  and the constraint hidden in \(\texttt{Penalty2}\) reads

  \[\begin{align*}
  \sum_{i,j} (\text{ServiceSelect}_{i,j}  - \text{ServiceSelect}_{i,j} *\text{FacilitySelect}_{i,j} = 0
  \end{align*}\]

  I don't know about the second QUBO formulation as no details about \(Q\) are given but I assume that something similar is happening.

  Best regards, 
  Jaromił

  0
• 

  Abhishek Solomon

  • March 07, 2023 16:21
  • Edited

  

  

  

  Thank you Jaromil. I do have constraints in the BIP model, i have just shown the objective function above. The constraints used are:

  for i in facilities:
      for j in customers:
          m.addConstr(ServiceSelect[i,j] <= FacilitySelect[i], 'Constraint1')
  
  for j in customers:
      m.addConstr(gp.quicksum((ServiceSelect[i,j]) for i in facilities) == 1 , 'Constraint2')

  For the second QUBO formulation, I am manually generating the Q matrix using logically assumed penalties based on research. Once the Q matrix is generated, i do the following:

  Q = np.array(Q)
  x = m.addVars(Q.shape[0], vtype=gp.GRB.BINARY)
  objective = gp.quicksum(x[i] * x[j] * Q[i, j] for i in range(Q.shape[0]) for j in range(Q.shape[0]))
  
  m.setObjective(objective, gp.GRB.MINIMIZE)

  What I am trying to understand is, are there any specific tuning/tweaking parameters that I need to invoke to solve QUBO models in Gurobi?

  0
• 

  Jaromił Najman

  Gurobi Staff

  • March 08, 2023 09:05
  • Edited

  In your BIP formulation, you have the constraints

  for i in facilities:
      for j in customers:
          m.addConstr(ServiceSelect[i,j] <= FacilitySelect[i], 'Constraint1')

  I don't see how these constraints translate into the penalty term

  Penalty2 = p*(gp.quicksum(ServiceSelect[i,j] - ServiceSelect[i,j]*FacilitySelect[i] for i in facilities for j in customers))

  Where does the bilinear term ServiceSelect[i,j]*FacilitySelect[i] come from?

  Could you share the logs for the 3 runs?

  What I am trying to understand is, are there any specific tuning/tweaking parameters that I need to invoke to solve QUBO models in Gurobi?

  You should definitely use the latest Gurobi version if you are solving QUBOs. You could try experimenting with the PreQLinearize parameter.

  0

Please sign in to leave a comment.