Examples of use
This chapter discusses some example Python scripts that are part of the Xpress Optimizer's Python interface. Most of them are well commented so the user can refer directly to the source for guidance.
Most of these scripts have an initial part in common, which we reproduce here but omit in all explanations below for compactness. These initial lines import the Xpress module itself and the NumPy module, which is used in some of the examples. The last line is to make the print statements, which are in Python 3 style here, work in Python 2.7 as well.
from __future__ import print_function import xpress as xp import numpy as np
Creating simple problems
Below are a few examples on how to create simple LP, MIP, MIQP, and similar problems. Note that they make use of API functions that resemble the C API functions for creating problems, and are used similarly here.
Generating a small Linear Programming problem
In this example, we create a problem and load a matrix of coefficients, a rhs, and an objective coefficient vector with the loadproblem function. We also assign names to both rows and columns (both are optional). These data correspond to the following problem with three variables and four constraints:
| minimize: | 3 x1 + 4 x2 + 5 x3 | ||
| subject to: | x1 + x3 | ≥ | -2.4 |
| 2x1 + 3x3 | ≥ | -3 | |
| 2x2 + 3x3 | = | 4 | |
| x2 + x3 | ≤ | 5 | |
| -1 ≤ x1 ≤ 3 | |||
| -1 ≤ x1 ≤ 5 | |||
| -1 ≤ x1 ≤ 8 |
p = xp.problem ()
p.loadproblem ("", # probname
['G','G','E', 'L'], # qrtypes
[-2.4, -3, 4, 5], # rhs
None, # range
[3,4,5], # obj
[0,2,4,8], # mstart
None, # mnel
[0,1,2,3,0,1,2,3], # mrwind
[1,2,2,1,1,3,3,1], # dmatval
[-1,-1,-1], # lb
[3,5,8], # ub
colnames = ['X1','X2','X3'], # column names
rownames = ['row1','row2','row3','constr_04']) # row names
p.write ("loadlp", "lp")
p.solve ()We then create another variable and add it to the problem, then modify the objective function. Note that the objective function is replaced by, not amended with, the new expression. After solving the problem, it saves it into a file called update.lp.
x = xp.var()
p.addVariable (x)
p.setObjective (x**2 + 2*x + 444)
p.solve()
p.write ("updated", "lp")A Mixed Integer Linear Programming problem
This example uses loadproblem to create a Mixed Integer Quadratically Constrained Quadratic Programming problem with two Special Ordered Sets. Note that data that is not needed is simply set as None.
The Examples directory provides similar examples for different types of problems.
p = xp.problem ()
p.loadproblem ("", # probname
['G','G','L', 'L'], # qrtypes
[-2.4, -3, 4, 5], # rhs
None, # range
[3,4,5], # obj
[0,2,4,8], # mstart
None, # mnel
[0,1,2,3,0,1,2,3], # mrwind
[1,1,1,1,1,1,1,1], # dmatval
[-1,-1,-1], # lb
[3,5,8], # ub
[0,0,0,1,1,2], # mqobj1
[0,1,2,1,2,2], # mqobj1
[2,1,1,2,1,2], # dqe
[2,3], # qcrows
[2,3], # qcnquads
[1,2,0,0,2], # qcmqcol1
[1,2,0,2,2], # qcmqcol2
[3,4,1,1,1], # qcdqval
['I','S'], # qgtype
[0,1], # mgcols
[0,2], # dlim
['1','1'], # qstype
[0,2,4], # msstart
[0,1,0,2], # mscols
[1.1,1.2,1.3,1.4]) # dref
p.solve ()Modeling examples
A simple model
This example demonstrates how variables and constraints, or arrays thereof, can be added into a problem. The script then prints the solution and all attributes/controls of the problem.
N = 4
S = range (N)
v = [xp.var (name = "y{0}".format (i)) for i in S] # set name of a variable as
m = xp.problem ()
v1 = xp.var (name = "v1", lb=0, ub=10, threshold=5, vartype = xp.continuous)
v2 = xp.var (name = "v2", lb=1, ub=7, threshold=3, vartype = xp.continuous)
vb = xp.var (name = "vb", vartype = xp.binary)
# Create a variable with name yi, where i is an index in S
v = [xp.var (name = "y{0}".format (i), lb = 0, ub = 2*N) for i in S]The call below adds both v, a vector (list) of variables, and v1 and v2, two scalar variables.
m.addVariable (vb, v, v1, v2)
c1 = v1 + v2 >= 5
m.addConstraint (c1, # Adds a list of constraints: three single constraints...
2*v1 + 3*v2 >= 5,
v[0] + v[2] >= 1,
# ... and a set of constraints indexed by all {i in
# S: i<N-1} (recall that ranges in Python are from 0
# to n-1)
(v[i+1] >= v[i] + 1 for i in S if i < N-1))
# objective overwritten at each setObjective ()
m.setObjective (xp.Sum ([i*v[i] for i in S]), sense = xp.minimize)
m.solve ()
print ("status: ", m.getProbStatus ())
print ("string: ", m.getProbStatusString ())
print ("solution:", m.getSolution ())Using IIS to investigate an infeasible problem
The problem modeled below is infeasible,
x0 = xp.var() x1 = xp.var() x2 = xp.var(vartype = xp.binary) c1 = x0 + 2 * x1 >= 1 c2 = 2 * x0 + x1 >= 1 c3 = x0 + x1 <= .5 c4 = 2 * x0 + 3 * x1 >= 0.1
The three constraints c1, c2, and c3 above are incompatible as can be easily verified. Adding all of them to a problem will make it infeasible. We use the functions to retrieve the Irreducible Infeasible Subsystems (IIS).
minf = xp.problem ("ex-infeas")
minf.addVariable (x0,x1,x2)
minf.addConstraint (c1,c2,c3,c4)
minf.solve()
minf.iisall()
print ("there are ", minf.attributes.numiis, " iis's")
miisrow = []
miiscol = []
constrainttype = []
colbndtype = []
duals = []
rdcs = []
isolationrows = []
isolationcols = []
# get data for the first IIS
minf.getiisdata (1, miisrow, miiscol, constrainttype, colbndtype,
duals, rdcs, isolationrows, isolationcols)
print ("iis data:", miisrow, miiscol, constrainttype, colbndtype,
duals, rdcs, isolationrows, isolationcols)
# Another way to check IIS isolations
print ("iis isolations:", minf.iisisolations (1))
rowsizes = []
colsizes = []
suminfeas = []
numinfeas = []
print ("iisstatus:", minf.iisstatus (rowsizes, colsizes, suminfeas, numinfeas))
print ("vectors:", rowsizes, colsizes, suminfeas, numinfeas)Modeling a problem using Python lists and vectors
We create a convex QCQP problem. We use a vector of N=5 variables and sets constraints and objective. We define all constraints and the objective function using a Python aggregate type.
N = 5
S = range (N)
v = [xp.var (name = "y{0}".format (i)) for i in S]
m = xp.problem ("problem 1")
print ("variable:", v)
m.addVariable (v)
m.addConstraint (v[i] + v[j] >= 1 for i in range (N-4) for j in range (i,i+4))
m.addConstraint (xp.Sum([v[i]**2 for i in range (N-1)]) <= N**2 * v[N-1]**2)
m.setObjective (xp.Sum ([i*v[i] for i in S]) * (xp.Sum ([i*v[i] for i in S])))
m.solve ()
print ("solution: ", m.getSolution ())A knapsack problem
Here follows an example of a knapsack problem formulated using lists of numbers. All data in the problem are lists, and so are the variables.
S = range (5) # that's the set {0,1,2,3,4}
value = [102, 512, 218, 332, 41] # or just read them from file
weight = [ 21, 98, 44, 59, 9]
x = [xp.var (vartype = xp.binary) for i in S]
profit = xp.Sum (value[i] * x[i] for i in S)
p = xp.problem ("knapsack")
p.addVariable (x)
p.addConstraint (xp.Sum (weight[i] * x[i] for i in S) <= 130)
p.setObjective (profit, sense = xp.maximize)
p.solve ()Note that the same result could have been achieved using NumPy arrays and the Xpress module's dot product as follows:
value = np.array ([102, 512, 218, 332, 41])
weight = np.array ([ 21, 98, 44, 59, 9])
x = np.array ([xp.var (vartype = xp.binary) for i in S])
profit = xp.Dot (value, x)
p = xp.problem ("knapsack")
p.addVariable (x)
p.addConstraint (xp.Dot (weight, x) <= 130)
p.setObjective (profit, sense = xp.maximize)
p.solve ()A Min-cost-flow problem using NumPy
This example solves a min-cost-flow problem using NumPy and the incidence matrix of the graph. It uses the networkx package, which might have to be installed using, for example, pip.
import networkx as netx # digraph definition V = [1,2,3,4,5] # vertices E = [[1,2], [1,4], [2,3], [3,4], [4,5], [5,1]] # arcs n = len (V) # number of nodes m = len (E) # number of arcs G = netx.DiGraph (E) # Get NumPy representation A = (netx.incidence_matrix (G, oriented = True).toarray())
We use NumPy vectors and the Xpress interface's dot product, the xpress.Dot operator. Note that although NumPy has a dot operator, especially for large models it is strongly advised to use the Xpress interface's Dot function for reasons of efficiency.
demand = np.array ([3, -5, 7, -2, -3])
cost = np.array ([23, 62, 90, 5, 6, 8])
flow = np.array ([xp.var () for i in E]) # flow variables declared on arcs
p = xp.problem ('network flow')
p.addVariable (flow)
p.addConstraint (xp.Dot (A, flow) == - demand)
p.setObjective (xp.Dot (cost, flow))
p.solve ()
for i in range (m):
print ('flow on', E[i], ':', p.getSolution (flow [i]))A nonlinear model
Let's solve a classical nonlinear problem: finding the minimum of the Rosenbrock function. For parameters a and b, minimize (a-x)2 + b (y-x2)2.
a,b = 1,100
x = xp.var (lb = -xp.infinity)
y = xp.var (lb = -xp.infinity)
p = xp.problem ()
p.addVariable (x,y)
p.setObjective ((a-x)**2 + b*(y-x**2)**2)
p.controls.xslp_solver = 0 # solve it with SLP, not Knitro
p.solve ()
print ("status: ", p.getProbStatus ())
print ("string: ", p.getProbStatusString ())
print ("solution:", p.getSolution ())Finding the maximum-area n-gon
The problem asks, given n, to find the n-sided polygon of largest area inscribed in the unit circle.
While it is natural to prove that all vertices of a global optimum reside on the unit circle, the problem is formulated so that every vertex i is at distance rhoi from the center, and at angle thetai. We would expect that the local optimum found has all rho's are equal to 1. The example file contains instructions for drawing the resulting polygon using matplotlib.
The objective function is the total area of the polygon. Considering the segment S[i] joining the center to the i-th vertex and A(i,j) the area of the triangle defined by the two segments S[i] and S[j], the objective function is A(0,1) + A(1,2) + ... + A(N-1,0), where A(i,j) = 1/2 * rhoi * rhoj * sin (thetai - thetaj). We first define the set Vertices as the set of integers from 0 to n-1.
rho = [xp.var (lb = 1e-5, ub = 1.0) for i in Vertices]
theta = [xp.var (lb = -math.pi, ub = math.pi) for i in Vertices]
p = xp.problem ()
p.addVariable (rho, theta)
p.setObjective (
0.5*(xp.Sum (rho[i]*rho[i-1]*xp.sin (theta[i]-theta[i-1]) for i in Vertices if i != 0)
+ rho[0]*rho[N-1]*xp.sin (theta[0]-theta[N-1])), sense = xp.maximize)We establish that the angles must be increasing in order to obtain a sensible solution:
p.addConstraint (theta[i] >= theta[i-1] + 1e-4 for i in Vertices if i != 0)
Note also that we enforce that the angles be different as otherwise they might form a local optimum where all of them are equal.
Solving the n-queens problem
In chess, the queen can move in all directions (even diagonally) and travel any distance. The problem of the n queens consists in placing n queens on an n×n chessboard so that none of them can be eaten in one move.
We first create a dictionary of variables, mapping each cell of the chessboard to one variable so that we can refer to it later. All variables are clearly binary as they indicate whether a given cell has a queen or not.
n = 10 # the size of the chessboard
N = range (n)
x = {(i,j): xp.var (vartype = xp.binary, name='q{0}_{1}'.format (i,j))
for i in N for j in N}
vertical = [xp.Sum (x[i,j] for i in N) <= 1 for j in N]
horizontal = [xp.Sum (x[i,j] for j in N) <= 1 for i in N]
diagonal1 = [xp.Sum (x[k-j,j] for j in range (max(0,k-n+1), min(k+1,n))) <= 1
for k in range (1,2*n-2)]
diagonal2 = [xp.Sum (x[k+j,j] for j in range (max(0,-k), min(n-k,n))) <= 1
for k in range (2-n,n-1)]
p = xp.problem()
p.addVariable (x)
p.addConstraint (vertical, horizontal, diagonal1, diagonal2)
# Objective, to be maximized: number of queens on the chessboard
p.setObjective (xp.Sum (x), sense = xp.maximize)
p.solve ()As a rudimentary form of visualization, we print the solution on the chessboard with different symbols for variables at one or zero.
for i in N:
for j in N:
if p.getSolution (x[i,j]) == 1:
print ('@', sep='', end='')
else:
print ('.', sep='', end='')
print ('')Solving Sudoku problems
The well-known Sudoku puzzles ask one to place numbers from 1 to 9 into a 9×9 grid such that no number repeats in any row, in any column, and in any 3x3 sub-grid. For a more general version of the game, replace 3 with q and 9 with q2.
We model this problem as an assignment problem where certain conditions must be met for all numbers in the columns, rows, and sub-grids.
These subgrids are lists of tuples with the coordinates of each subgrid. In a 9×9 sudoku, for instance, subgrids[0,1] has the 9 elements in the middle top square.
The input is a starting grid where the unknown numbers are replaced by zero. The example file contains a relatively hard 9×9 sudoku, which we show below, and also a 16×16 variant of the same game.
q = 3
starting_grid = \
[[8,0,0,0,0,0,0,0,0],
[0,0,3,6,0,0,0,0,0],
[0,7,0,0,9,0,2,0,0],
[0,5,0,0,0,7,0,0,0],
[0,0,0,0,4,5,7,0,0],
[0,0,0,1,0,0,0,3,0],
[0,0,1,0,0,0,0,6,8],
[0,0,8,5,0,0,0,1,0],
[0,9,0,0,0,0,4,0,0]]
n = q**2 # the size must be the square of the size of the subgrids
N = range (n)
x = {(i,j,k): xp.var (vartype = xp.binary, name='x{0}_{1}_{2}'.format(i,j,k))
for i in N for j in N for k in N}
# define all q^2 subgrids
subgrids = {(h,l): [(i,j) for i in range (q*h, q*h + q)
for j in range (q*l, q*l + q)]
for h in range (q) for l in range (q)}
vertical = [xp.Sum (x[i,j,k] for i in N) == 1 for j in N for k in N]
horizontal = [xp.Sum (x[i,j,k] for j in N) == 1 for i in N for k in N]
subgrid = [xp.Sum (x[i,j,k] for (i,j) in subgrids[h,l]) == 1
for (h,l) in subgrids.keys() for k in N]
# Assign exactly one number to each cell
assign = [xp.Sum (x[i,j,k] for k in N) == 1 for i in N for j in N]Then we fix those variables that are non-zero in the input grid. We don't need an objective function as this is a feasibility problem. After computing the solution, we print it to the screen.
init = [x[i,j,k] == 1 for k in N for i in N for j in N
if starting_grid[i][j] == k+1]
p = xp.problem()
p.addVariable (x)
p.addConstraint (vertical, horizontal, subgrid, assign, init)
p.solve ()
print ('Solution:')
for i in N:
for j in N:
l = [k for k in N if p.getSolution (x[i,j,k]) >= 0.5]
assert (len (l) == 1)
print ('{0:2d}'.format (1 + l[0]), end = '', sep='')
print ('')Examples using NumPy
Using NumPy multidimensional arrays to create variables
Use NumPy arrays for creating a 3-dimensional array of variables, then use it to create a mode.
S1 = range (2)
S2 = range (3)
S3 = range (4)
m = xp.problem ()
h = np.array ([[[xp.var (vartype = xp.binary)
for i in S1]
for j in S2]
for k in S3])
m.addVariable (h)
m.setObjective (h[0][0][0] * h[0][0][0] +
h[1][0][0] * h[0][0][0] +
h[1][0][0] * h[1][0][0] +
xp.Sum (h[i][j][k] for i in S3 for j in S2 for k in S1))
cons00 = - h[0][0][0] ** 2 +
xp.Sum (i * j * k * h[i][j][k]for i in S3 for j in S2 for k in S1) >= 11
m.addConstraint (cons00)
m.solve ()The final part of the code retrieves the matrix representation of the quadratic part of the only constraint.
mstart1=[]
mclind1=[]
dqe1=[]
m.getqrowqmatrix (cons00, mstart1, mclind1, dqe1, 29, h[0][0][0], h[3][2][1])
print ("row 0:", mstart1, mclind1, dqe1)Using the dot product to create arrays of expressions
Here we use NumPy arrays to print the product of a matrix by a random vector, and the xpress.Dot function on a matrix and a vector. Note that the NumPy dot operator works perfectly fine here, but should be avoided for reasons of performance, especially when handling large arrays where at least one contains optimization variables or expressions.
x = np.array ([xp.var() for i in range(5)]) p = xp.problem () p.addVariable (x) p.addConstraint (xp.Sum (x) >= 2) p.setObjective (xp.Sum (x[i]**2 for i in range (5))) p.solve() A = np.array (range(30)).reshape(6,5) # A is a 6x5 matrix sol = np.array (p.getSolution ()) # a vector of size 5 columns = A*sol # not a matrix-vector product! v = np.dot (A,sol) # an array: matrix-vector product A*sol w = xp.Dot (A,x) # an array of expressions print (v,w)
Using the Dot product to create constraints and quadratic functions
This is an example of a problem formulation that uses the xpress.Dot operator to formulate constraints in a concise fashion. Note that the NumPy dot operator is not suitable here as the result is an expression in the Xpress variables.
A = np.random.random(30).reshape(6,5) # A is a 6x5 matrix Q = np.random.random(25).reshape(5,5) # Q is a 5x5 matrix x = np.array ([xp.var() for i in range(5)]) # vector of variables x0 = np.random.random (5) # random vector Q += 4 * np.eye (5) # add 5 * the identity matrix Lin_sys = xp.Dot (A,x) <= np.array([3,4,1,4,8,7]) # 6 constraints (rows of A) Conv_c = xp.Dot (x,Q,x) <= 1 # one quadratic constraint p = xp.problem () p.addVariable (x) p.addConstraint (Lin_sys, Conv_c) p.setObjective (xp.Dot (x-x0, x-x0)) # minimize distance from x0 p.solve ()
Using NumPy to create quadratic optimization problems
This example creates and solves a simple quadratic optimization problem. Given an n×n matrix Q and a point x0, minimize the quadratic function xT (Q + n3 I) x subject to the linear system (x - x0)T Q + e = 0, where e is the vector of all ones, the inequalities Q x ≥ 0, and nonnegativity on all variables. Report solution if available.
n = 10
Q = np.arange (1, n**2 + 1).reshape (n, n)
x = np.array ([xp.var () for i in range (n)])
x0 = np.random.random (n)
p = xp.problem ()
p.addVariable (x)
c1 = xp.Dot ((x - x0), Q) + 1 == 0
c2 = xp.Dot (Q, x) >= 0
p.addConstraint (c1,c2)
p.setObjective (xp.Dot (x, Q + N**3 * np.eye (N), x))
p.solve ('')
print ("nrows, ncols:", p.attributes.rows, p.attributes.cols)
print ("solution:", p.getSolution ())
p.write ("test5-qp", "lp")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 + cardWe 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] = newnodeWe 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 subtreesQuery 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 x02 - 2 x0 x3 + x13 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 |
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 NoneThe 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 0We 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