How to solve continuous-variable quadratic program and do warm-stating
Awaiting user input
I am going to solve a continuous-variable Qubo problem by Gurobi. However, it gives me always the binary solutions. Does any one know how to resolve the issue such that I can get continuous-value solutions? Moreover, I do not know how to warm start the algorithm? I need to give some initial points as initial solutions to the code. Would you please help me how to do warm-starting?
Please find the code in the following
import gurobipy as gp
import numpy as np
gurobipy_model = gp.Model("gurobi")
from qiskit_optimization.translators import from_gurobipy
from qiskit_optimization.algorithms import GurobiOptimizer
Mu = np.array([[-1.71712467],
[-1.80263717],
[-2.01173903],
[ 0.39771078],
[-1.52744357],
[-0.67927752],
[-0.32842294],
[ 1.02471285],
[-0.52852204],
[ 1.10914181],
[ 0.05156664],
[-0.2383112 ]])
Sigma = np.array([[ 0. , -0.27204021, 0. , 0. , -0.60856918,
-0.83651528, 0. , 0. , 0. , 0. ,
0. , 0. ],
[-0.27204021, 0. , -0.69628892, 0. , 0. ,
0. , 0. , 0. , -0.83430803, 0. ,
0. , 0. ],
[ 0. , -0.69628892, 0. , -0.43617927, -0.87927083,
0. , 0. , 0. , 0. , 0. ,
0. , 0. ],
[ 0. , 0. , -0.43617927, 0. , 0. ,
0. , 0. , -0.14450387, 0.97839393, 0. ,
0. , 0. ],
[-0.60856918, 0. , -0.87927083, 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
-0.03960356, 0. ],
[-0.83651528, 0. , 0. , 0. , 0. ,
0. , -0.61454932, 0.77178708, 0. , 0. ,
0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
-0.61454932, 0. , 0. , 0. , 0.99996394,
0. , -0.71383757],
[ 0. , 0. , 0. , -0.14450387, 0. ,
0.77178708, 0. , 0. , 0. , 0.39742965,
0. , 0. ],
[ 0. , -0.83430803, 0. , 0.97839393, 0. ,
0. , 0. , 0. , 0. , 0. ,
-0.67260794, 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0.99996394, 0.39742965, 0. , 0. ,
0. , -0.28825178],
[ 0. , 0. , 0. , 0. , -0.03960356,
0. , 0. , 0. , -0.67260794, 0. ,
0. , 0.76377815],
[ 0. , 0. , 0. , 0. , 0. ,
0. , -0.71383757, 0. , 0. , -0.28825178,
0.76377815, 0. ]])
x = [gurobipy_model.addVar(vtype=gp.GRB.CONTINUOUS, name="x%s" % i, lb=0, ub=1) for i in range(len(Sigma))]
objective = sum([Mu[i][0] * x[i] for i in range(len(Mu))])
objective -= 2 * sum(
[Sigma[i][j] * x[i] * x[j] for i in range(len(Mu)) for j in range(len(Mu))]
)
gurobipy_model.setObjective(objective,gp. GRB.MAXIMIZE)
qp2 = from_gurobipy(gurobipy_model)
GurobiOptimizer().solve(qp2)
The above code leads to the following solutions:
<OptimizationResult: fval=13.271412780000002, x0=1.0, x1=1.0, x2=1.0, x3=0.0, x4=1.0, x5=1.0, x6=1.0, x7=0.0, x8=1.0, x9=0.0, x10=1.0, x11=0.0, status=SUCCESS>
-
Hi Shadi,
There is nothing wrong with your code. Your variables are continuous but the optimal solution just naturally has integer values. If there was a better solution with non-integer values then it would find it - but there isn't.
You can provide a valid starting extreme point to warm start, either through PStart, DStart, or through VBasis, CBasis. Note that this disables LP presolve by default. For models where presolve greatly reduces the problem size, this might hurt performance. To allow presolve, you need to set parameter LPWarmStart to 2.
- Riley
0 -
Thanks for the reply and your great help.
I expect to have at least some continuous values where I stop Gurobi where the MIPGap is still large. For example, in the following code:
import numpy as np
from numpy.linalg import inv
import quimb as qu
import quimb.tensor as qtn
import networkx as nx
import itertools
import tqdm
from scipy.sparse import csc_matrix
from itertools import product
import numba as nb
import copy
import matplotlib.pyplot as plt
from scipy import optimize
Adj_G = np.array([[0. , 0.05219673, 0. , 0.16303623, 0. ,
0.67521474, 0. , 0. , 0. , 0. ,
0. , 0. , 0.09561873, 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.55474926],
[0.05219673, 0. , 0. , 0. , 0. ,
0. , 0. , 0.49808213, 0. , 0.53599452,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0.04826614,
0. , 0. , 0. , 0.51098376, 0. ,
0. ],
[0. , 0. , 0. , 0.9012575 , 0.2084062 ,
0. , 0. , 0. , 0. , 0.62677833,
0.19434701, 0. , 0. , 0. , 0. ,
0. , 0. , 0.37677331, 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. ],
[0.16303623, 0. , 0.9012575 , 0. , 0. ,
0. , 0. , 0. , 0. , 0.88774306,
0. , 0.84029947, 0. , 0. , 0. ,
0. , 0. , 0. , 0.41976977, 0. ,
0. , 0. , 0. , 0. , 0. ,
0. ],
[0. , 0. , 0.2084062 , 0. , 0. ,
0.87952865, 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0.92751916, 0. , 0. ,
0. , 0. , 0.86352158, 0.64232798, 0. ,
0. ],
[0.67521474, 0. , 0. , 0. , 0.87952865,
0. , 0.5057895 , 0. , 0. , 0. ,
0.80035643, 0. , 0. , 0.24959981, 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. ],
[0. , 0. , 0. , 0. , 0. ,
0.5057895 , 0. , 0.45582632, 0.26956553, 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.64960479, 0. ,
0. , 0. , 0. , 0. , 0.66420765,
0. ],
[0. , 0.49808213, 0. , 0. , 0. ,
0. , 0.45582632, 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0.05319689, 0. , 0. , 0.33955944,
0. , 0. , 0.79072681, 0. , 0. ,
0. ],
[0. , 0. , 0. , 0. , 0. ,
0. , 0.26956553, 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.95428089, 0.80772909, 0.5564221 , 0. , 0. ,
0.78341913],
[0. , 0.53599452, 0.62677833, 0.88774306, 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0.45334776,
0.00274732, 0. , 0. , 0. , 0. ,
0. ],
[0. , 0. , 0.19434701, 0. , 0. ,
0.80035643, 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0.02048179, 0. , 0. , 0. ,
0. , 0. , 0.43266406, 0.56560779, 0. ,
0. ],
[0. , 0. , 0. , 0.84029947, 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0.89235174, 0. , 0. ,
0.47988292, 0. , 0. , 0. , 0. ,
0.66394688, 0. , 0. , 0. , 0. ,
0.40707133],
[0.09561873, 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0.89235174, 0. , 0.38101796, 0. ,
0. , 0. , 0.24381057, 0. , 0. ,
0.82976059, 0. , 0. , 0. , 0. ,
0. ],
[0. , 0. , 0. , 0. , 0. ,
0.24959981, 0. , 0. , 0. , 0. ,
0. , 0. , 0.38101796, 0. , 0.30463268,
0. , 0. , 0. , 0.85735497, 0. ,
0. , 0. , 0. , 0. , 0.14984741,
0. ],
[0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.30463268, 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0.66473423, 0.78331818, 0. , 0.55290625,
0.10918973],
[0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0.47988292, 0. , 0. , 0. ,
0. , 0.64922511, 0.54986399, 0. , 0. ,
0. , 0.29496819, 0. , 0. , 0.98257993,
0. ],
[0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0.05319689, 0. , 0. ,
0.02048179, 0. , 0. , 0. , 0. ,
0.64922511, 0. , 0. , 0.66010914, 0. ,
0. , 0. , 0. , 0.57230448, 0. ,
0. ],
[0. , 0. , 0.37677331, 0. , 0.92751916,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0.24381057, 0. , 0. ,
0.54986399, 0. , 0. , 0. , 0. ,
0. , 0.53402091, 0. , 0. , 0. ,
0. ],
[0. , 0. , 0. , 0.41976977, 0. ,
0. , 0.64960479, 0. , 0. , 0. ,
0. , 0. , 0. , 0.85735497, 0. ,
0. , 0.66010914, 0. , 0. , 0.66693744,
0. , 0. , 0. , 0. , 0. ,
0. ],
[0. , 0.04826614, 0. , 0. , 0. ,
0. , 0. , 0.33955944, 0. , 0.45334776,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.66693744, 0. ,
0. , 0. , 0. , 0. , 0. ,
0.79111979],
[0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.95428089, 0.00274732,
0. , 0.66394688, 0.82976059, 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0.89126883, 0. , 0. , 0. ,
0. ],
[0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.80772909, 0. ,
0. , 0. , 0. , 0. , 0.66473423,
0.29496819, 0. , 0.53402091, 0. , 0. ,
0.89126883, 0. , 0. , 0. , 0. ,
0. ],
[0. , 0. , 0. , 0. , 0.86352158,
0. , 0. , 0.79072681, 0.5564221 , 0. ,
0.43266406, 0. , 0. , 0. , 0.78331818,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. ],
[0. , 0.51098376, 0. , 0. , 0.64232798,
0. , 0. , 0. , 0. , 0. ,
0.56560779, 0. , 0. , 0. , 0. ,
0. , 0.57230448, 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0.49486887,
0. ],
[0. , 0. , 0. , 0. , 0. ,
0. , 0.66420765, 0. , 0. , 0. ,
0. , 0. , 0. , 0.14984741, 0.55290625,
0.98257993, 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.49486887, 0. ,
0. ],
[0.55474926, 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.78341913, 0. ,
0. , 0.40707133, 0. , 0. , 0.10918973,
0. , 0. , 0. , 0. , 0.79111979,
0. , 0. , 0. , 0. , 0. ,
0. ]])
#Find the adjacency matrix of the graph
G = nx.Graph() # or DiGraph, MultiGraph, MultiDiGraph, etc
G.add_weighted_edges_from(terms)
Adj_G = (nx.adjacency_matrix(G)).todense()
#Generate single-term entries of Qubo matrix
Mu = np.zeros((len(Adj_G),1))
for i in range(0, len(Adj_G)):
Mu[i] = sum(Adj_G[i,:])
Mu = np.transpose(Mu)
Sigma = Adj_G.copy()import gurobipy as gp
import numpy as np
import matplotlib.pyplot as plt# Create a Gurobi model
model = gp.Model()# Create variables
x = [model.addVar(vtype=gp.GRB.CONTINUOUS, name="x%s" % i, lb=0, ub=1) for i in range(len(Sigma))]
objective = sum([Mu[0][i] * x[i] for i in range(len(Sigma))])
objective -= 2 * sum(
[Sigma[i][j] * x[i] * x[j] for i in range(len(Sigma)) for j in range(len(Sigma))]
)
# Set objective: quadratic function (non-convex)
model.setObjective(objective, gp.GRB.MAXIMIZE)
model.setParam("NonConvex", 2)
# Optimize the model with callback to track objective values
objective_values = []def callback(model, where):
if where == gp.GRB.Callback.MIP:
current_obj = model.cbGet(gp.GRB.Callback.MIP_OBJBST)
objective_values.append(current_obj)
model.setParam('outputFlag',True)
# Set the callback function
model.optimize(callback)
#model.optimize()for v in model.getVars():
print('%s: %g' % (v.varName, v.x))This code gives me the following output:
Optimize a model with 0 rows, 26 columns and 0 nonzeros Model fingerprint: 0x6903cd3e Model has 65 quadratic objective terms Coefficient statistics: Matrix range [0e+00, 0e+00] Objective range [2e+00, 4e+00] QObjective range [2e-02, 8e+00] Bounds range [1e+00, 1e+00] RHS range [0e+00, 0e+00] Continuous model is non-convex -- solving as a MIP Found heuristic solution: objective -0.0000000 Presolve time: 0.01s Presolved: 65 rows, 91 columns, 195 nonzeros Variable types: 0 continuous, 91 integer (91 binary) Root relaxation: objective 3.191450e+01, 27 iterations, 0.00 seconds (0.00 work units) Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 31.91450 0 24 -0.00000 31.91450 - - 0s H 0 0 23.7008622 31.91450 34.7% - 0s H 0 0 24.6478244 31.91450 29.5% - 0s H 0 0 25.1587044 31.91450 26.9% - 0s H 0 0 25.4479741 28.63059 12.5% - 0s 0 0 28.57396 0 25 25.44797 28.57396 12.3% - 0s 0 0 27.48074 0 32 25.44797 27.48074 7.99% - 0s H 0 0 25.7761577 27.48074 6.61% - 0s Cutting planes: Gomory: 1 MIR: 2 Zero half: 21 RLT: 5 Explored 1 nodes (70 simplex iterations) in 0.03 seconds (0.00 work units) Thread count was 8 (of 8 available processors) Solution count 6: 25.7762 25.448 25.1587 ... -0 Optimal solution found (tolerance 1.00e-04) Best objective 2.577615770146e+01, best bound 2.577615770146e+01, gap 0.0000% User-callback calls 138, time in user-callback 0.00 sec x0: 1 x1: 0 x2: -0 x3: 0 x4: 1 x5: 0 x6: -0 x7: 1 x8: 1 x9: 1 x10: 1 x11: 0 x12: 1 x13: -0 x14: 1 x15: 1 x16: 0 x17: -0 x18: 1 x19: 0 x20: -0 x21: 0 x22: 0 x23: 0 x24: 0 x25: 0
Then, I try to stop it where the Mipgap is still large by adding the command:
model.setParam("MIPGap", 3.47), However, I still get binaries.
So, even in the early stage of optimization process, it gives me binary.
I am solving maxcut problem. Even for large number of variables,
it gives me binary in the early stage of optimization!!0 -
Hi Shadi,
As you can see in the log, the presolve routine is creating a model in which the variables are binary.
See "How does presolve work?"
If you turn presolve off:
model.params.presolve=0
then you can avoid this behavior (although I don't see a benefit in doing this).
If you were to set a solution limit of 2, for example, :
model.params.SolutionLimit=2
then you will see you get fractional solutions found - however don't expect the optimal solution to be fractional.
- Riley
0 -
Thanks a lot for your great help and the nice clarification.
Sure. You are right. The pre-solver considers the variables as binaries.
I deactivate the pre-solver in the early stage of process using
import gurobipy as gp
import numpy as np
import matplotlib.pyplot as plt# Create a Gurobi model
model = gp.Model()# Create variables
x = [model.addVar(vtype=gp.GRB.CONTINUOUS, name="x%s" % i, lb=0, ub=1) for i in range(len(Sigma))]
objective = sum([Mu[0][i] * x[i] for i in range(len(Sigma))])
objective -= 2 * sum(
[Sigma[i][j] * x[i] * x[j] for i in range(len(Sigma)) for j in range(len(Sigma))]
)
# Set objective: quadratic function (non-convex)
model.setObjective(objective, gp.GRB.MAXIMIZE)
model.setParam("NonConvex", 2)
model.params.presolve=0
model.params.SolutionLimit=2
# Optimize the model with callback to track objective values
objective_values = []def callback(model, where):
if where == gp.GRB.Callback.MIP:
current_obj = model.cbGet(gp.GRB.Callback.MIP_OBJBST)
objective_values.append(current_obj)
model.setParam('outputFlag',True)
# Set the callback function
model.optimize(callback)
#model.optimize()for v in model.getVars():
print('%s: %g' % (v.varName, v.x))that gives me the following outputlog:
Set parameter NonConvex to value 2 Set parameter Presolve to value 0 Set parameter SolutionLimit to value 3 Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (mac64[arm]) CPU model: Apple M1 Thread count: 8 physical cores, 8 logical processors, using up to 8 threads Optimize a model with 0 rows, 20 columns and 0 nonzeros Model fingerprint: 0x439bdf76 Model has 190 quadratic objective terms Coefficient statistics: Matrix range [0e+00, 0e+00] Objective range [8e+00, 1e+01] QObjective range [2e-02, 8e+00] Bounds range [1e+00, 1e+00] RHS range [0e+00, 0e+00] Continuous model is non-convex -- solving as a MIP Found heuristic solution: objective -0.0000000 Variable types: 211 continuous, 0 integer (0 binary) Root relaxation: objective 9.623557e+01, 86 iterations, 0.00 seconds (0.00 work units) Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 96.23557 0 190 -0.00000 96.23557 - - 0s H 0 0 6.2015125 96.23557 1452% - 0s H 0 0 38.7463939 96.23557 148% - 0s Explored 1 nodes (86 simplex iterations) in 0.01 seconds (0.00 work units) Thread count was 8 (of 8 available processors) Solution count 3: 38.7464 6.20151 -0 Solution limit reached Best objective 3.874639385111e+01, best bound 9.623557229450e+01, gap 148.3730% User-callback calls 48, time in user-callback 0.00 sec x0: 2.32926e-10 x1: 1.07374e-09 x2: 2.36949e-10 x3: 2.24515e-10 x4: 1 x5: 5.20858e-10 x6: 1 x7: 3.94854e-10 x8: 1 x9: 1 x10: 2.14363e-10 x11: 5.63529e-10 x12: 1 x13: 5.96567e-10 x14: 2.28581e-10 x15: 1 x16: 3.25358e-10 x17: 1 x18: 3.79257e-10 x19: 2.68775e-10However, I still get the solutions that are very close to binaries. Would you please let me know if you think the solver rounds the solutions to binaries? If yes, is there ant way to avoid that and get fractional values. Thanks a lot for your help
0 -
Hi Shadi,
Your code:
model.params.SolutionLimit=2
The log:
Set parameter SolutionLimit to value 3
What went wrong here?
However, I still get the solutions that are very close to binaries. Would you please let me know if you think the solver rounds the solutions to binaries?
In your example you have printed out the solutions from the solver, which are not rounded, but you're asking if the solver rounds them? I'm confused.
If you are wanting a non-optimal non-integral solution then there are many ways you could get one, including just forcing it by setting lb = ub = 0.5 on all the variables:
x = [model.addVar(vtype=gp.GRB.CONTINUOUS, name="x%s" % i, lb=0.5, ub=0.5) for i in range(len(Sigma))]
Perhaps we need to address what you are ultimately trying to do, because a lot of this isn't making sense to me. Why are you wanting a non-optimal non-integer solution?
- Riley
0
Please sign in to leave a comment.
Comments
5 comments