How to solve 0/1-QP by Gurobi

I will introduce how to use gurobi to solve a QP with 0/1 variables on python env.

Gurobi is a canonical solver to solve integer programming. I postulate you have mastered basic operations for python. This blog mainly contains three part:

1) Present QP with 0/1 variables;

2) Introduce how to install Gurobi and apply Academic license;

3) Give a specific code to solve given problem.

Now, let's begin.

I Problem model

We focus on this QP problem with 0/1 variables(01-QP):

                                                                  (1)

The main troble of this problem is 0/1 variables. Strictly, this is a NP hard problem. However, there are many commercial  software that can solve it with heuristic algorithm, such as Gurobi and cplex. By the way, Gurobi is the most famous solver.

II. How to install Gurobi and apply Academic license

I recommaned you to use colab provied by Google, since Gurobi is very nicely  embedded in colab. Next we just need to apply a academic license to solve large scale problem.

Following the instruction step by step of https://support.gurobi.com/hc/en-us/articles/4409582394769-Google-Colab-Installation-and-Licensing

More specifically, apply an account in https://www.gurobi.com/account/    

Sign in and click on the following red circle

Then click on the following red circle 

 Notice that your school maybe unable to allow to get a License sometimes.

III Sovle it.

In this part, we solve 01-QP(1) by python use gurobi on colab with m=n=50.

!pip install gurobipy  # install gurobipy, if not already installed
import gurobipy as gp
# Create environment with WLS license
import pandas as pd
from gurobipy import GRB

e = gp.Env(empty=True)
e.setParam('WLSACCESSID', '****************')
e.setParam('WLSSECRET', '***************')
e.setParam('LICENSEID', ******)
e.start()
# # Create the model within the Gurobi environment
model = gp.Model(env=e)

m = gp.Model("miqp", env=e)

# m = gp.Model("qp")
import numpy as np
# Create variables
n = 50
np.random.seed(2)
A0 = 10 * np.random.uniform(-2 / 10, 2 / 3, size=(n, n))
A = np.ceil(A0)
x_or = 10 * np.random.uniform(-5, 5, size=(n, 1))
x_or = np.ceil(x_or)
id_in = (np.where(x_or > 0.5))[0]
x_or[id_in] /= (np.linalg.norm(x_or[id_in], axis = 1)).reshape((len(id_in),1))
x_or *= (x_or > 0)
b = np.ceil(np.dot(A, x_or))

# A =   10000 * np.array([[2, 0, 9], [7, 1, 10], [0, 1, 1]])
# b0 = 10000 * np.array([9, 11, 2])
Q = np.dot(A.T, A)
b1 = np.dot(b.T, A)
x = m.addMVar(n, vtype=GRB.BINARY)
print("Test{}.".format(np.dot(b.T, b)))
obj = x @ Q  @ x - 2 * b1 @ x + np.dot(b.T, b)
m.setObjective(obj, gp.GRB.MINIMIZE)
m.update()
# # Add constraint: x + 2 y + 3 z >= 4
# m.addConstr(x + 2 * y + 3 * z >= 4, "c0")

# # Add constraint: x + y >= 1
# m.addConstr(x + y >= 1, "c1")

m.optimize()

print('Obj: %g' % obj.getValue())
x_best = (x.x).reshape((n,1))

print("error:{}.".format(np.linalg.norm(x_or - x_best)))

That's the end. Have a good day, dude.