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 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.

import networkx as netx
from matplotlib import pyplot as plt

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.

T = nx.Graph ()

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

p = xp.problem ()
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.

p = xp.problem ()

N = 5
S = range (N)

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

p.addVariable (x)

# vectors can be used whole or addressed with their index

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 vector 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 + 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 add 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 vector 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.

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 to an empty string 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