Advanced examples: callbacks and problem querying, modifying, and analysis

Visualize the branch-and-bound tree of a problem

This example shows how to visualize the BB tree of a problem after (partially) solving it. It is assumed here that all branches are binary.

We first define a message callback for running code whenever the Optimizer wants to print a message. The callback receives four arguments: the problem and callback data and, most importantly, the message to be printed and an information number. The callback prints the output message prefixed by a time stamp related to the creation of the problem. As the message could be on multiple lines, it is split into multiple substrings, one per line.

import networkx as nx
import time
from matplotlib import pyplot as plt

def message_addtime (prob, data, msg, info):
    """Message callback example: print a timestamp before the message from the optimizer"""
    if msg:
        for submsg in msg.split('\n'):
            print("{0:6.3f}: [{2:+4d}] {1}".format(time.time() - start_time, submsg, info))

We then define a recursive function that computes the cardinality of a subtree rooted at a node i. This is necessary as the visualization of the BB tree is more balanced when the subtree size is taken into account. The card_subtree array, which is filled here, is used then for computing the width of each visualized subtree.

def postorder_count(node):
    """
    Recursively count nodes to compute the cardinality of a subtree for
    each node
    """

    card = 0

    if node in left.keys(): # see if node has a left key
        postorder_count(left[node])
        card += card_subtree[left[node]]

    if node in right.keys():
        postorder_count(right[node])
        card += card_subtree[right[node]]

    card_subtree[node] = 1 + card

We also define a function that determines the position of each node depending on the cardinality of the subtree rooted at the node.

def setpos(T, node, curpos, st_width, depth):

    """
    Set position depending on cardinality of each subtree
    """

    # Special condition: we are at the root
    if node == 1:
        T.add_node(node, pos=(0.5, 1))

    alpha = .1 # use a convex combination of subtree comparison and
               # depth to assign a width to each subtree

    if node in left.keys():

        # X position in the graph should not just depend on depth,
        # otherwise we'd see a long and thin subtree and it would just
        # look like a path

        leftwidth = st_width * (alpha * .5 + (1 - alpha) * card_subtree[left[node]]
                    / card_subtree[node])
        leftpos = curpos - (st_width - leftwidth) / 2

        T.add_node(left[node], pos=(leftpos, - depth))
        T.add_edge(node, left[node])
        setpos(T, left[node], leftpos, leftwidth, depth + 1)

    if node in right.keys():

        rightwidth = st_width * (alpha * .5 + (1 - alpha) * card_subtree[right[node]]
                    / card_subtree[node])
        rightpos = curpos + (st_width - rightwidth) / 2

        T.add_node(right[node], pos=(rightpos, - depth))
        T.add_edge(node, right[node])
        setpos(T, right[node], rightpos, rightwidth, depth + 1)

This is the only operation we need to be carried out at every node: given a node number, newnode, and its parent, parent, we store the information in the left and right arrays so that at the end of the BB we have an explicit BB tree stored in these arrays.

def storeBBnode(prob, Tree, parent, newnode, branch):
    # Tree is the callback data, and it's equal to T

    if branch == 0:
        left[parent] = newnode
    else:
        right[parent] = newnode

We now set up the BB tree data and create a problem. We read it from a local file, but any user problem can be read and analyzed. We set the node callback with addcbnewnode so that we can collect information at each new node. We also save the initial time for use by message_addtime, the function that is called every time the problem prints out a message.

T = nx.Graph()

left  = {}
right = {}
card_subtree = {}
pos = {}

start_time = time.time()

p = xp.problem()
p.addcbmessage(message_addtime)

p.read('sampleprob.mps.gz')
p.addcbnewnode(storeBBnode, T, 100)
p.controls.maxnode=40000 # Limit the number of nodes inserted in the graph
p.solve()

postorder_count(1)      # assign card_subtree to each node
setpos(T, 1, 0.5, 1, 0) # determine the position of each node
                         # depending on subtree cardinalities

pos = nx.get_node_attributes(T, 'pos')

nx.draw(T, pos) # create BB tree representation
plt.show()      # display it; you can zoom indefinitely and see all subtrees

Query and modify a simple problem

This example shows how to change an optimization problem using the Xpress Python interface.

x = xp.var()
y = xp.var()

cons1    =   x + y >= 2
upperlim = 2*x + y <= 3

p = xp.problem()

p.addVariable(x,y)
p.setObjective((x-4)**2 + (y-1)**2)
p.addConstraint(cons1, upperlim)

p.write('original', 'lp')

After saving the problem to a file, we change two of its coefficients. Note that the same operations can be carried out with a single call to p.chgmcoef([cons1,1],[x,0],[3,4]).

p.chgcoef(cons1, x, 3) # coefficient of x in cons1    becomes 3
p.chgcoef(1, 0, 4)     # coefficient of y in upperlim becomes 4

p.write('changed', 'lp')

Change a problem after solution

Construct a problem using addVariable and addConstraint, then use the Xpress API routines to amend the problem with rows and quadratic terms.

import xpress as xp

p = xp.problem()
N = 5
S = range(N)

x = [xp.var(vartype=xp.binary) for i in S]

p.addVariable(x)

# Vectors of variables can be used whole or addressed with an index or
# index range

c0 = xp.Sum(x) <= 10
cc = [x[i]/1.1 <= x[i+1]*2 for i in range(N-1)]

p.addConstraint(c0, cc)

p.setObjective(3 - x[0])

mysol = [0, 0, 1, 1, 1, 1.4]

# add a variable with its coefficients

p.addcols([4], [0,3], [c0,4,2], [-3, 2.4, 1.4], [0], [2], ['Y'], ['B'])
p.write("problem1", "lp")

# load a MIP solution
p.loadmipsol([0,0,1,1,1,1.4])

We now add a quadratic term x₀² - 2 x₀ x₃ + x₁³ to the second constraint. Note that the -2 coefficient for an off-diagonal element must be passed divided by two.

p.addqmatrix(cc[0], [x[0],x[3],x[3]], [x[0],x[0],x[3]], [1,-1,1])

As constraint list cc was added after c0, it is the latter which has index 0 in the problem, while cc[0] has index 1. We then add the seventh and eighth constraints:

subject to:	x0 + 2 x1 + 3x2	≥	4
	4x0 + 5x1 + 6x2 + 7 x3 + 8 x4 -3 y	≤	4.4

Note the new column named 'Y' is added with its index 5 (variables' indices begin at 0). The same would happen if 5 were substituted by Y.

p.addqmatrix(1, [x[0],x[3],x[3]], [x[0],x[0],x[3]], [1,-1,1])

p.addrows(qrtype=['G', 'L'],
           rhs=[4, 4.4],
           mstart=[0, 3, 9],
           mclind=[x[0],x[1],x[2], x[0],x[1],x[2],x[3],x[4], 5],
           dmatval=[1,2,3,4,5,6,7,8,-3],
           names=['newcon1', 'newcon2'])

p.solve()
p.write("amended", "lp")

slacks = []

p.calcslacks(solution=mysol, calculatedslacks=slacks)

print("slacks:", slacks)

The code below first adds five columns, then solves the problem and prints the solution, if one has been found.

p.addcols([4], [0,3], [c0,4,2], [-3,  -2,   1], [0],  [2], ['p1'], ['I'])
p.addcols([4], [0,3], [c0,4,2], [-3, 2.4, 1.4], [0], [10], ['p2'], ['C'])
p.addcols([4], [0,3], [c0,4,2], [-3,   2,   1], [0],  [1], ['p3'], ['S'])
p.addcols([4], [0,3], [c0,4,2], [-3, 2.4,   4], [0],  [2], ['p4'], ['P'])
p.addcols([4], [0,3], [c0,4,2], [-3,   2,   1], [0],  [2], ['p5'], ['R'])

p.solve()

try:
    print("new solution:", p.getSolution())
except:
    print("could not get solution, perhaps problem is infeasible")

Note that the single command below has the same effect as the four addcols calls above, and is to be preferred when adding a large number of columns for reasons of efficiency.

p.addcols([4,4,4,4,4],
           [0,3,6,9,12,15],
           [c0,4,2,c0,4,2,c0,4,2,c0,4,2,c0,4,2],
           [3, -2, 1, -3, 2.4, 1.4, 3, 2, 1, -3, 2.4, 4, 3, 2, 1],
           [0,0,0,0,0],
           [2,10,1,2,2],
           ['p1','p2','p3','p4','p5'],
           ['I','C','S','P','R'])

Comparing the coefficients of two equally sized problems

Given two problems with the same number of variables, we read their coefficient matrices into Scipy so as to compare each row for discrepancies in the coefficients. We begin by creating two Xpress problems and reading them from two files, prob1.lp and prob2.lp, though p1 and p2 might have been created with the module's modeling features.

import xpress as xp
import scipy.sparse

p1 = xp.problem()
p2 = xp.problem()

p1.read('prob1.lp')
p2.read('prob2.lp')

Next we obtain the matrix representation of the coefficient matrix for both problems. Let us suppose that, for memory reasons, we can only retrieve one million coefficients.

coef1, ind1, beg1 = [], [], []
coef2, ind2, beg2 = [], [], []

p1.getrows(beg1, ind1, coef1, 1000000, 0, p1.attributes.rows - 1)
p2.getrows(beg2, ind2, coef2, 1000000, 0, p2.attributes.rows - 1)

The function problem.getrows provides a richer output by filling up ind1 and ind2 with the Python objects (i.e. Xpress variables) corresponding to the variable indices rather than the numerical indices. We need to convert them to numerical indices using the problem.getIndex function.

ind1n = [p1.getIndex(v) for v in ind1]
ind2n = [p2.getIndex(v) for v in ind2]

The next step is to create a Compressed Sparse Row (CSR) format matrix, defined in the scipy.sparse module, using the data from problem.getrows plus the numerical indices.

Then we convert the CSR matrix to a NumPy array of arrays, so that each row is a (non-compressed) array to be compared in the loop below.

A1 = scipy.sparse.csr_matrix((coef1, ind1n, beg1))
A2 = scipy.sparse.csr_matrix((coef2, ind2n, beg2))

M1 = A1.toarray()
M2 = A2.toarray()

for i in range(min(p1.attributes.rows, p2.attributes.rows)):
    print(M1[i] != M2[i])

The result is a few vectors of size COLS with an element-wise comparison of the coefficient vector of each row, with True indicating discrepancies. A more meaningful representation can be given using other functions in NumPy.

[False False  True False False]
[False False False False False]
[False False False False  True]
[ True  True False False False]
[False False False False False]

Combining modeling and API functions

This is an example where a problem is loaded from a file, solved, then modified by adding a Global Upper Bound (GUB) constraint. Note that we do not know the structure of the problem when reading it, yet we can simply extract the list of variables and use them to add a constraint.

import xpress
p = xpress.problem()

p.read("example.lp")
p.solve()
print("solution of the original problem: ", p.getVariable(), "==>", p.getSolution())

After solving the problem, we obtain its variables through getVariable and add a constraints so that their sum cannot be more than 1.1.

x = p.getVariable()
p.addConstraint(xpress.Sum(x) <= 1.1)
p.solve()
print("New solution: ", p.getSolution())

A simple Traveling Salesman Problem (TSP) solver

A classical example of use of callbacks is the development of a simple solver for the well-known TSP problem. The aim here is not to create an efficient solver (there are far better implementations), but rather a simple solver where the user only needs to specify two callbacks: one for checking whether a given solution forms a Hamiltonian tour and one for separating a subtour elimination constraint from the current node solution.

After a successful solve (or an interrupted one with a feasible solution), the best Hamiltonian tour is displayed. Note that this section omits unnecessary details (checks of return values, exceptions, etc.) of the actual code, which can be found in the Examples/ directory.

import networkx as nx
import xpress   as xp
import re, math, sys

from matplotlib import pyplot as plt

import urllib.request as ul

filename = 'dj38.tsp'

ul.urlretrieve('http://www.math.uwaterloo.ca/tsp/world/' + filename, filename)

instance = open(filename, 'r')
coord_section = False
points = {}

G = nx.Graph()

We have downloaded an instance of the TSP and now it must be read and interpreted as it does not have a format that we know. We save in cx and cy the coordinates of all nodes in the graph, which is assumed to be complete, i.e., all nodes are connected to one another.

for line in instance.readlines():

    if re.match('NODE_COORD_SECTION.*', line):
        coord_section = True
        continue
    elif re.match('EOF.*', line):
        break

    if coord_section:
        coord = line.split(' ')
        index = int(coord[0])
        cx    = float(coord[1])
        cy    = float(coord[2])
        points[index] = (cx, cy)
        G.add_node(index, pos=(cx, cy))

The next step is to define a callback function for checking if the solution forms a Hamiltonian tour, i.e., if it connects all nodes of the graph. The callback will be passed with the method addcbpreintsol, therefore it needs to return a tuple of two values: the first value is True if the solution should be rejected, and the second is the new cutoff in case it has to be changed. This is not the case here, so None can be safely returned.

After obtaining the integer solution to be checked, the function scans the graph from node 1 to see if the solutions at one form a tour.

def check_tour(prob, G, isheuristic, cutoff):

    s = []

    prob.getlpsol(s, None, None, None)

    orignode = 1
    nextnode = 1
    card     = 0

    while nextnode != orignode or card == 0:

        FS = [j for j in V if j != nextnode
              and s[prob.getIndex(x[nextnode,j])] == 1] # forward star
        card += 1

        if len(FS) < 1:
            return (True, None) # reject solution if we can't close the loop

        nextnode = FS[0]

    # If there are n arcs in the loop, the solution is feasible

    return (card < n, None) # accept the cutoff: return second element as None

The second callback to be defined is a separator for subtour elimination constraints. It must return a nonzero value if the node is deemed infeasible by the function, zero otherwise. The function addcuts is used to insert a subtour elimination constraint.

The function works as follows: Starting from node 1, gather all connected nodes of a loop in connset. If this set contains all nodes, then the solution is valid if integer, otherwise the function adds a subtour elimination constraint in the form of a clique constraint with all arcs (i,j) for all i,j in connset.

def eliminate_subtour(prob, G):

    s = [] # initialize s as an empty list to provide it as an output parameter

    prob.getlpsol(s, None, None, None)

    orignode = 1
    nextnode = 1

    connset  = []

    while nextnode != orignode or len(connset) == 0:

        connset.append(nextnode)

        FS = [j for j in V if j != nextnode
              and s[prob.getIndex(x[nextnode, j])] == 1] # forward star

        if len(FS) < 1:
            return 0

        nextnode = FS[0]

    if len(connset) < n:

        # Add a subtour elimination using the nodes in connset (or, if
        # card(connset) > n/2, its complement)


        if len(connset) <= n/2:
            columns = [x[i,j] for i in connset for j in connset
                       if i != j]
            nArcs  = len(connset)
        else:
            columns = [x[i,j] for i in       V for j in       V
                       if not i in connset and not j in connset and i != j]
            nArcs  = n - len(connset)

        nTerms = len(columns)

        prob.addcuts([1], ['L'], [nArcs - 1], [0, nTerms], columns, [1] * nTerms)

    return 0

We now formulate the problem with the degree constraints on each node and the objective function (the cost of each arc (i,j) is assumed to be the Euclidean distance between i and j).

n = len(points)                             # number of nodes
V = range(1, n+1)                            # set of nodes
A = [(i,j) for i in V for j in V if i != j]  # set of arcs (i.e. all pairs)

x = {(i,j): xp.var(name='x_{0}_{1}'.format(i,j), vartype=xp.binary) for (i,j) in A}

conservation_in  = [xp.Sum(x[i,j] for j in V if j != i) == 1 for i in V]
conservation_out = [xp.Sum(x[j,i] for j in V if j != i) == 1 for i in V]

p = xp.problem()

p.addVariable(x)
p.addConstraint(conservation_in, conservation_out)

xind = {(i,j): p.getIndex(x[i,j]) for (i,j) in x.keys()}

# Objective function: total distance travelled
p.setObjective(xp.Sum(math.sqrt((points[i][0] - points[j][0])**2 +
                                  (points[i][1] - points[j][1])**2) *
                                   x[i,j]
                       for (i,j) in A))

p.controls.maxtime = -2000 # negative for "stop even if no solution is found"

p.addcboptnode(eliminate_subtour, G, 1)
p.addcbpreintsol(check_tour,        G, 1)

We now solve the problem, and if a solution is found it is displayed using the Python library matplotlib.

p.solve()

sol = p.getSolution()

# Read solution and store it in the graph

for (i,j) in A:
    if sol[p.getIndex(x[i,j])] > 0.5:
        G.add_edge(i,j)

# Display best tour found

pos = nx.get_node_attributes(G, 'pos')

nx.draw(G, points) # create a graph with the tour
plt.show()         # display it interactively

Another solver for TSP problems is available in example_tsp_numpy.py. The two main differences consist in the problem generation, which is now random, and in the fact that most data structures are NumPy vectors and matrices: the optimization variables, the LP solution obtained from the Branch-and-Bound, and the data used to check feasibility of the solutions.

Solving a nonconvex MIQCQP

In this example we turn the Xpress Optimizer into a solver for nonconvex MIQCQPs, i.e. problems with nonconvex quadratic objective and/or nonconvex quadratic constraints.

In order to handle nonconvex quadratic constraints, we have to reformulate the problem to a MILP so that the simplest nonlinear terms, i.e. the products of variables, are transformed into new, so-called auxiliary variables.

Product x_ix_j is assigned to a new variable w_ij so that every occurrence of that product in the problem is replaced by w_ij. Assuming l_i and u_i are the lower and upper bound on x_i, respectively, we add the linear McCormick inequalities:

• w_ij ≥ l_j x_i + l_i x_j - l_j l_i
• w_ij ≥ u_j x_i + u_i x_j - u_j u_i
• w_ij ≤ l_j x_i + u_i x_j - l_j u_i
• w_ij ≤ u_j x_i + l_i x_j - u_j l_i

The bounds on the new auxiliary variable w_ij are a function of the bounds on x_i and x_j.

Below is the code that takes care of reformulating the problem. We first have to identify all terms x_ix_j and create a dictionary linking each pair (i,j) to an auxiliary variable w_ij. The dictionary aux is used throughout the solver and contains this information. The function create_prob checks all bilinear terms and creates aux and the McCormick inequalities.

def create_prob(filename):

    [...]

    x = p.getVariable()

    aux = {}  # Dictionary containing the map (x_i,x_j) --> y_ij

    [...]

    p.addConstraint(
        [aux[i, j] >= lb[j]*x[i] + lb[i]*x[j] - lb[i] * lb[j]
         for (i, j) in aux.keys() if max(-lb[i], -lb[j]) < xp.infinity],
        [aux[i, j] >= ub[j]*x[i] + ub[i]*x[j] - ub[i] * ub[j]
         for (i, j) in aux.keys() if max(ub[i], ub[j]) < xp.infinity],
        [aux[i, j] <= ub[j]*x[i] + lb[i]*x[j] - lb[i] * ub[j]
         for (i, j) in aux.keys() if max(-lb[i], ub[j]) < xp.infinity],
        [aux[i, j] <= lb[j]*x[i] + ub[i]*x[j] - ub[i] * lb[j]
         for (i, j) in aux.keys() if max(ub[i], -lb[j]) < xp.infinity])

We also needs to tell the Optimizer that the newly created auxiliary variables and the variables that used to appear in bilinear terms should be protected against deletion by the presolver.

    securecols = list(aux.values())
    secureorig = set()

    for i, j in aux.keys():
        secureorig.add(i)
        secureorig.add(j)

    securecols += list(secureorig)

    p.loadsecurevecs(mrow=None, mcol=securecols)

The creation of a single auxiliary variable is done in addaux, where its bounds are created and, depending on whether it is the product of two variables or the square of one, it receives a different treatment.

def addaux(aux, p, i, j, lb, ub, vtype):

    # Find bounds of auxiliary first
    if i != j:
        # bilinear term
        l, u = bdprod(lb[i], ub[i], lb[j], ub[j])
    elif lb[i] >= 0:
        l, u = lb[i]**2, ub[i]**2
    elif ub[i] <= 0:
        l, u = ub[i]**2, lb[i]**2
    else:
        l, u = 0, max([lb[i]**2, ub[i]**2])

After setting the bounds on w_ij, we determine its type and create the corresponing xp.var object.

    if vtype[i] == 'B' and vtype[j] == 'B':
        t = xp.binary
    elif (vtype[i] == 'B' or vtype[i] == 'I') and \
         (vtype[j] == 'B' or vtype[j] == 'I'):
        t = xp.integer
    else:
        t = xp.continuous

    # Add auxiliaries
    aux[i, j] = xp.var(lb=l, ub=u, vartype=t,
                       name='aux_{0}_{1}'.format(
                           p.getVariable(i).name,
                           p.getVariable(j).name))

    return aux[i, j]

Quadratic constraints and the quadratic objective (if any) are converted in convQaux, where they are replaced by a linear expression containing auxiliary variables.

def convQaux(p, aux, mstart, ind, coef, row, lb, ub, vtype):

    rcols = []
    rrows = []
    rcoef = []

    for i,__ms in enumerate(mstart[:-1]):
        for j in range(mstart[i], mstart[i+1]):

            J = p.getIndex(ind[j])

            if (i, J) not in aux.keys():
                y = addaux(aux, p, i, J, lb, ub, vtype)
                p.addVariable(y)
            else:
                y = aux[i, J]

            if row < 0:  # objective
                mult = .5
            else:
                mult = 1

            if i != J:
                coe = 2 * mult * coef[j]
            else:
                coe =     mult * coef[j]

            if row < 0:
                p.chgobj([y], [coe])
            else:
                rcols.append(y)
                rrows.append(row)
                rcoef.append(coe)

    if row >= 0:

        # This is a quadratic constraint, not the objective function
        # Add linear coefficients for newly introduced variables
        p.chgmcoef(rrows, rcols, rcoef)
        # Remove quadratic matrix
        p.delqmatrix(row)

    else:

        # Objective: Remove quadratic part
        indI = []
        for i in range(len(mstart) - 1):
            indI.extend([i] * (mstart[i+1] - mstart[i]))
        # Set all quadratic elements to zero
        p.chgmqobj(indI, ind, [0] * mstart[-1])

The new problem, called a reformulation, is then solved as a MILP with a few callbacks. Given that the problem is nonconvex, we need to branch on continuous variables, those that appear in bilinear terms, and we also need to keep adding McCormick inequalities when the bounds change. This is because in branch-and-bound algorithms for nonconvex problems the linear relaxation should be exact at the extremes of the variable bound ranges.

Another callback is to decide whether to accept or not a solution that was found by the branch-and-bound: because the constraints linking w to x are missing, we must make sure that they are satisfied by a solution, and must refuse a solution that does not satisfy w_ij = x_ix_j.

def solveprob(p, aux):

    p.addcbpreintsol(cbchecksol, aux, 1)
    p.addcboptnode(cbaddcuts, aux, 3)
    p.addcbchgbranchobject(cbbranch, aux, 1)

    p.mipoptimize()

The callback functions are fundamental. The branch callback checks whether the auxiliary variables w_ij are satisfied, and if not it creates a branching object on either x_i or x_j. Due to the presolved nature of the problem at this point in the branch-and-bound, care must be applied in handling the variable indices, as they might have changed by the presolver to allow for a smaller problem.

def cbbranch(prob, aux, branch):

    lb, ub = getCBbounds(prob, len(sol))

    x = prob.getVariable()  # presolved variables

    rowmap = []
    colmap = []

    prob.getpresolvemap(rowmap, colmap)

    invcolmap = [-1 for _ in lb]

    for i, m in enumerate(colmap):
        invcolmap[m] = i

    # Check if all auxiliaries are equal to their respective bilinear
    # term. If so, we have a feasible solution

    sol = np.array(sol)

    discr = sol[Aux_ind] - sol[Aux_i] * sol[Aux_j]
    discr[Aux_i == Aux_j] = np.maximum(0, discr[Aux_i == Aux_j])
    maxdiscind = np.argmax(np.abs(discr))

    if abs(discr[maxdiscind]) < eps:
        return branch

    i,j = Aux_i[maxdiscind], Aux_j[maxdiscind]

    yind = prob.getIndex(aux[i, j])

For terms of the form w_ii = x_i², branching might still be necessary as the curve defining it is a nonconvex set.

    if i == j:

        # Test of violation is done on the original
        # space. However, the problem variables are scrambled with invcolmap

        if sol[i] > lb[i] + eps and \
           sol[i] < ub[i] - eps and \
           sol[yind] > sol[i]**2 + eps and \
           sol[yind] - lb[i]**2 <= (ub[i] + lb[i]) * (sol[i] - lb[i]) - eps:

            # Can't separate, must branch. Otherwise OA or secant
            # cut separated above should be enough

            brvarind = invcolmap[i]
            brpoint = sol[i]
            brvar = x[brvarind]
            brleft = brpoint
            brright = brpoint

            assert(brvarind >= 0)

            if brvar.vartype in [xp.integer, xp.binary]:
                brleft = math.floor(brpoint + 1e-5)
                brright = math.ceil(brpoint - 1e-5)

            b = xp.branchobj(prob, isoriginal=False)

            b.addbranches(2)

            addrowzip(prob, b, 0, 'L', brleft,  [i], [1])
            addrowzip(prob, b, 1, 'G', brright, [i], [1])

            # New variable bounds are not enough, add new McCormick
            # inequalities for y = x**2: suppose x0,y0 are the current
            # solution values for x,y, yp = x0**2 and xu,yu = xu**2 are their
            # upper bound, and similar for lower bound. Then these two
            # rows must be added, one for each branch:
            #
            # y - yp <= (yl-yp)/(xl-x0) * (x - x0)  <===>
            # (yl-yp)/(xl-x0) * x - y >= (yl-yp)/(xl-x0) * x0 - yp
            #
            # y - yp <= (yu-yp)/(xu-x0) * (x - x0)  <===>
            # (yu-yp)/(xu-x0) * x - y >= (yu-yp)/(xu-x0) * x0 - yp
            #
            # Obviously do this only for finite bounds

            ypl = brleft**2
            ypr = brright**2

            if lb[i] > -1e7 and sol[i] > lb[i] + eps:

                yl = lb[i]**2
                coeff = (yl - ypl) / (lb[i] - sol[i])

                if coeff != 0:
                    addrowzip(prob, b, 0, 'G', coeff*sol[i] - ypl,
                              [i, yind], [coeff, -1])

            if ub[i] < 1e7 and sol[i] < ub[i] - eps:

                yu = ub[i]**2
                coeff = (yu - ypr) / (ub[i] - sol[i])

                if coeff != 0:
                    addrowzip(prob, b, 1, 'G', coeff*sol[i] - ypr,
                              [i, yind], [coeff, -1])

            return b

Similarly for bilinear terms, we must choose where to branch and on which variable.

    else:

        lbi0, ubi0 = lb[i], ub[i]
        lbi1, ubi1 = lb[i], ub[i]

        lbj0, ubj0 = lb[j], ub[j]
        lbj1, ubj1 = lb[j], ub[j]

        # No cut violated, must branch
        if min(sol[i] - lb[i], ub[i] - sol[i]) / (1 + ub[i] - lb[i]) > \
           min(sol[j] - lb[j], ub[j] - sol[j]) / (1 + ub[j] - lb[j]):
            lbi1 = sol[i]
            ubi0 = sol[i]
            brvar = i
        else:
            lbj1 = sol[j]
            ubj0 = sol[j]
            brvar = j

        alpha = 0.2

        brvarind = invcolmap[brvar]
        brpoint = sol[brvar]
        brleft = brpoint
        brright = brpoint

        if x[brvarind].vartype in [xp.integer, xp.binary]:
            brleft = math.floor(brpoint + 1e-5)
            brright = math.ceil(brpoint - 1e-5)

        b = xp.branchobj(prob, isoriginal=False)

        b.addbranches(2)

        addrowzip(prob, b, 0, 'L', brleft,  [brvar], [1])
        addrowzip(prob, b, 1, 'G', brright, [brvar], [1])

        # As for the i==j case, the variable branch is
        # insufficient, so add updated McCormick inequalities.
        # There are two McCormick inequalities per changed bound:
        #
        # y >= lb[j] * x[i] + lb[i] * x[j] - lb[j] * lb[i] ---> add to branch 1
        # y >= ub[j] * x[i] + ub[i] * x[j] - ub[j] * ub[i] ---> add to branch 0
        # y <= lb[j] * x[i] + ub[i] * x[j] - lb[j] * ub[i] ---> add to branch 1 if x[brvarind] == j, 0 if x[brvarind] == i
        # y <= ub[j] * x[i] + lb[i] * x[j] - ub[j] * lb[i] ---> add to branch 1 if x[brvarind] == i, 0 if x[brvarind] == j

        addrowzip(prob, b, 0, 'G', - ubi0 * ubj0, [yind, i, j], [1, -ubj0, -ubi0])
        addrowzip(prob, b, 1, 'G', - lbi1 * lbj1, [yind, i, j], [1, -lbj1, -lbi1])

        if brvarind == i:
            addrowzip(prob, b, 0, 'L', - lbj0 * ubi0, [yind, i, j], [1, -lbj0, -ubi0])
            addrowzip(prob, b, 1, 'L', - ubj1 * lbi1, [yind, i, j], [1, -ubj1, -lbi1])
        else:
            addrowzip(prob, b, 0, 'L', - ubj0 * lbi0, [yind, i, j], [1, -ubj0, -lbi0])
            addrowzip(prob, b, 1, 'L', - lbj1 * ubi1, [yind, i, j], [1, -lbj1, -ubi1])
        return b

    # If no branching rule was found, return none
    return branch

The callback for checking a solution is straightforward: for all pairs ij, check if the corresponding identity w_ij = x_i x_j is satisfied, and if not, simply reject the solution.

def cbchecksol(prob, aux, soltype, cutoff):

    global Aux_i, Aux_j, Aux_ind

    if (prob.attributes.presolvestate & 128) == 0:
        return (1, cutoff)

    sol = []

    # Retrieve node solution
    try:
        prob.getlpsol(x=sol)
    except:
        return (1, cutoff)

    sol = np.array(sol)

    # Check if all auxiliaries are equal to their respective bilinear
    # term. If so, we have a feasible solution

    refuse = 1 if np.max(np.abs(sol[Aux_i] * sol[Aux_j] - sol[Aux_ind])) > eps else 0

    # Return with refuse != 0 if solution is rejected, 0 otherwise;
    # and same cutoff
    return (refuse, cutoff)

An important part of this nonconvex solver is a function that computes a new feasible solution. The one we attempt here is rather trivial and probably not able to find good solutions, but one could add other algorithms, which for example might just use an alternative solver, and find a feasible solution, regardless of how good.

def cbfindsol(prob, aux):

    sol = []

    try:
        prob.getlpsol(x=sol)
    except:
        return 0

    xnew = sol[:]

    # Round solution to nearest integer
    for i,t in enumerate(var_type):
        if t == 'I' or t == 'B' and \
           xnew[i] > math.floor(xnew[i] + prob.controls.miptol) + prob.controls.miptol:
            xnew[i] = math.floor(xnew[i] + .5)

    for i, j in aux.keys():
        yind = prob.getIndex(aux[i, j])
        xnew[yind] = xnew[i] * xnew[j]

    prob.addmipsol(xnew)

    return 0

The function for adding McCormick inequalities is perhaps the most important as it allows for the lower bound in the branch-and-bound to get tighter at every node. All violated inequalities are added for all pairs ij.

def cbaddmccormickcuts(prob, aux, sol):
    lb, ub = getCBbounds(prob, len(sol))

    cuts = []

    # Check if all auxiliaries are equal to their respective bilinear
    # term. If so, we have a feasible solution
    for i, j in aux.keys():

        yind = prob.getIndex(aux[i, j])

        if i == j:

            # Separate quadratic term

            if sol[yind] < sol[i]**2 - eps and \
               abs(sol[i]) < xp.infinity / 2:

                xk = sol[i]

                ox = xk
                oy = ox ** 2

                # Add Outer Approximation cut y >= xs^2 + 2xs*(x-xs)
                # <===> y - 2xs*x >= -xs^2
                cuts.append((TYPE_OA, 'G', - ox**2, [yind, i],
                             [1, -2*ox]))

            # Otherwise, check if secant can be of help: y0 - xl**2 >
            # (xu**2 - xl**2) / (xu - xl) * (x0 - xl)
            elif sol[yind] > sol[i]**2 + eps and \
                 sol[yind] - lb[i]**2 > (ub[i] + lb[i]) * (sol[i] - lb[i]) \
                 + eps and abs(lb[i] + ub[i]) < xp.infinity / 2:
                cuts.append((TYPE_SECANT, 'L',
                             lb[i]**2 - (ub[i] + lb[i]) * lb[i],
                             [yind, i], [1, - (lb[i] + ub[i])]))

        elif abs(sol[yind] - sol[i]*sol[j]) > eps:

            # Separate bilinear term, where i != j.  There might be at
            # least one cut violated

            if sol[yind] < lb[j]*sol[i] + lb[i]*sol[j] - lb[i]*lb[j] - eps:
                if lb[i] > -xp.infinity / 2 and lb[j] > -xp.infinity / 2:
                    cuts.append((TYPE_MCCORMICK, 'G', - lb[i] * lb[j],
                                 [yind, i, j], [1, -lb[j], -lb[i]]))
            elif sol[yind] < ub[j]*sol[i] + ub[i]*sol[j] - ub[i]*ub[j] - eps:
                if ub[i] < xp.infinity / 2 and ub[j] < xp.infinity / 2:
                    cuts.append((TYPE_MCCORMICK, 'G', - ub[i] * ub[j],
                                 [yind, i, j], [1, -ub[j], -ub[i]]))
            elif sol[yind] > lb[j]*sol[i] + ub[i]*sol[j] - ub[i]*lb[j] + eps:
                if ub[i] < xp.infinity / 2 and lb[j] > -xp.infinity / 2:
                    cuts.append((TYPE_MCCORMICK, 'L', - ub[i] * lb[j],
                                 [yind, i, j], [1, -lb[j], -ub[i]]))
            elif sol[yind] > ub[j]*sol[i] + lb[i]*sol[j] - lb[i]*ub[j] + eps:
                if lb[i] > -xp.infinity / 2 and ub[j] < xp.infinity / 2:
                    cuts.append((TYPE_MCCORMICK, 'L', - lb[i] * ub[j],
                                 [yind, i, j], [1, -ub[j], -lb[i]]))

    # Done creating cuts. Add them to the problem

    for (t, s, r, I, C) in cuts:  # cuts might be the empty list
        mcolsp, dvalp = [], []
        drhsp, status = prob.presolverow(s, I, C, r, prob.attributes.cols,
                                         mcolsp, dvalp)
        if status >= 0:
            prob.addcuts([t], [s], [drhsp], [0, len(mcolsp)], mcolsp, dvalp)

    return 0

Another useful component of any nonconvex solver is a procedure to tighten the variable bounds based on information that is known on other variables. For example, if new bounds are inferred on w_ij, possible tighter lower or upper bounds can be deduced on x_i and/or x_j.

def cbboundreduce(prob, aux, sol):
    cuts = []

    lb, ub = getCBbounds(prob, len(sol))

    # Check if bounds on original variables can be reduced based on
    # bounds on auxiliary ones. The other direction is already taken
    # care of by McCormick and tangent/secant cuts.

    feastol = prob.controls.feastol

    for (i,j),a in aux.items():

        auxind = prob.getIndex(a)

        lbi = lb[i]
        ubi = ub[i]
        lba = lb[auxind]
        uba = ub[auxind]

        if i == j:  # check if upper bound is tight w.r.t. bounds on
                    # x[i]

            # Forward propagation: from new independent variable
            # bounds, infer new bound for dependent variable.

            if uba > max(lbi**2, ubi**2) + feastol:
                cuts.append((TYPE_BOUNDREDUCE, 'L', max(lbi**2, ubi**2), [auxind], [1]))

            if lbi > 0 and lba < lbi**2 - feastol:
                cuts.append((TYPE_BOUNDREDUCE, 'G', lbi**2, [auxind], [1]))
            elif ubi < 0 and lba < ubi**2 - feastol:
                cuts.append((TYPE_BOUNDREDUCE, 'G', ubi**2, [auxind], [1]))

            if uba < -feastol:
                return 1  # infeasible node
            else:
                if uba < lbi**2 - feastol:
                    if lbi > 0:
                        return 1  # infeasible node
                    else:
                        cuts.append((TYPE_BOUNDREDUCE, 'G', -math.sqrt(uba), [i], [1]))
                if uba < ubi**2 - feastol:
                    if ubi < - feastol:
                        return 1
                    else:
                        cuts.append((TYPE_BOUNDREDUCE, 'L', math.sqrt(uba), [i], [1]))

            if lba > prob.controls.feastol and lbi > 0 and lbi**2 < lba - feastol:
                cuts.append((TYPE_BOUNDREDUCE, 'G', math.sqrt(lba), [i], [1]))

        else:

            tlb, tub = bdprod(lb[i], ub[i], lb[j], ub[j])

            if lba < tlb - feastol:
                cuts.append((TYPE_BOUNDREDUCE, 'G', tlb, [auxind], [1]))

            if uba > tub + feastol:
                cuts.append((TYPE_BOUNDREDUCE, 'L', tub, [auxind], [1]))

            # For simplicity let's just assume lower bounds are nonnegative

            lbj = lb[j]
            ubj = ub[j]

            if lbj >= 0 and lbi >= 0:

                if lbi*ubj < lba - feastol:
                    cuts.append((TYPE_BOUNDREDUCE, 'G', lba / ubj, [i], [1]))
                if lbj*ubi < lba - feastol:
                    cuts.append((TYPE_BOUNDREDUCE, 'G', lba / ubi, [j], [1]))

                if lbi*ubj > uba + feastol:
                    cuts.append((TYPE_BOUNDREDUCE, 'L', uba / lbi, [j], [1]))
                if lbj*ubi > uba + feastol:
                    cuts.append((TYPE_BOUNDREDUCE, 'L', uba / lbj, [i], [1]))

    # Done creating cuts. Add them to the problem

    for (t, s, r, I, C) in cuts:  # cuts might be the empty list

        mcolsp, dvalp = [], []
        drhsp, status = prob.presolverow(s, I, C, r, prob.attributes.cols,
                                         mcolsp, dvalp)
        if status >= 0:

            if len(mcolsp) == 0:
                continue
            elif len(mcolsp) == 1:
                if s == 'G':
                    btype = 'L'
                elif s == 'L':
                    btype = 'U'
                else:  # don't want to add an equality bound reduction
                    continue

                assert(dvalp[0] > 0)

                prob.chgbounds(mcolsp,[btype],[drhsp/dvalp[0]])
            else:
                prob.addcuts([t], [s], [drhsp], [0, len(mcolsp)], mcolsp, dvalp)

    return 0

There are a few other functions not shown here that are used in the example. These are functions for retrieving bounds withing a callback and other service functions. The example file provides commented code that can be used to improve the solver.