Examples of use
Topics covered in this chapter:
- Creating simple problems
- Modeling examples
- Examples using NumPy
- Advanced examples: callbacks and problem querying, modifying, and analysis
- Translated Mosel examples
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.
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 list 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'], # rowtype
[-2.4, -3, 4, 5], # rhs
None, # rng
[3,4,5], # objcoef
[0,2,4,8], # start
None, # collen
[0,1,2,3,0,1,2,3], # rowind
[1,2,2,1,1,3,3,1], # rowcoef
[-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.optimize() 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 = p.addVariable()
p.setObjective(x**2 + 2*x + 444)
p.optimize()
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'], # rowtype
[-2.4, -3, 4, 5], # rhs
None, # rng
[3,4,5], # objcoef
[0,2,4,8], # start
None, # collen
[0,1,2,3,0,1,2,3], # rowind
[1,1,1,1,1,1,1,1], # rowcoef
[-1,-1,-1], # lb
[3,5,8], # ub
[0,0,0,1,1,2], # objqcol1
[0,1,2,1,2,2], # objqcol2
[2,1,1,2,1,2], # objqcoef
[2,3], # qrowind
[2,3], # nrowqcoefs
[1,2,0,0,2], # rowqcol1
[1,2,0,2,2], # rowqcol2
[3,4,1,1,1], # rowqcoef
['I','S'], # coltype
[0,1], # entind
[0,2], # limit
['1','1'], # settype
[0,2,4], # setstart
[0,1,0,2], # setind
[1.1,1.2,1.3,1.4]) # refval
p.optimize() Modeling examples
A simple model
This example demonstrates how variables and constraints, or lists/arrays thereof, can be added into a problem. The script then prints the solution and attributes of the problem.
N = 4
S = range(N)
m = xp.problem()
v = [m.addVariable(name="y{0}".format(i), lb=0, ub=2*N) for i in S]
v1 = m.addVariable(name="v1", lb=0, ub=10, threshold=5, vartype=xp.continuous)
v2 = m.addVariable(name="v2", lb=1, ub=7, threshold=3, vartype=xp.semicontinuous)
vb = m.addVariable(name="vb", vartype=xp.binary)
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)
solvestatus, solstatus = m.optimize()
print("solve status: ", solvestatus.name)
print("solution status: ", solstatus.name)
print("solution:", m.getSolution()) Using IIS to investigate an infeasible problem
The problem modeled below is infeasible,
import xpress as xp
minf = xp.problem("ex-infeas")
x0 = minf.addVariable()
x1 = minf.addVariable()
x2 = minf.addVariable(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.addConstraint(c1,c2,c3,c4)
minf.optimize()
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 list of N=5 variables and sets constraints and objective. We define all constraints and the objective function using a Python aggregate type.
import xpress as xp
N = 5
S = range(N)
m = xp.problem("problem 1")
v = [m.addVariable(name="y{0}".format(i)) for i in S]
print("variable:", 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.optimize()
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.
import xpress as xp
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]
p = xp.problem("knapsack")
x = p.addVariables(5, vartype=xp.binary)
profit = xp.Sum(value[i] * x[i] for i in S)
p.addConstraint(xp.Sum(weight[i] * x[i] for i in S) <= 130)
p.setObjective(profit, sense=xp.maximize)
p.optimize() Note that the same result could have been achieved using NumPy arrays and the Xpress module's dot product as follows:
import xpress as xp
import numpy as np
value = np.array([102, 512, 218, 332, 41])
weight = np.array([ 21, 98, 44, 59, 9])
p = xp.problem("knapsack")
x = p.addVariables(5, vartype=xp.binary)
profit = xp.Dot(value, x)
p.addVariable(x)
p.addConstraint(xp.Dot(weight, x) <= 130)
p.setObjective(profit, sense=xp.maximize)
p.optimize() A Min-cost-flow problem using NumPy
This example solves a min-cost-flow problem using NumPy and the incidence matrix of the graph.
import numpy as np import xpress as xp # 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
We then generate the incidence matrix by creating a NumPy matrix with n rows and m columns, such that each column, which corresponds to an arc (i,j), has a -1 at row i and a 1 at row j.
# Generate incidence matrix: begin with a NxM zero matrix
A = np.zeros((n,m))
# Then for each column i of the matrix, add a -1 in correspondence to
# the tail of the arc and a 1 for the head of the arc. Because Python
# uses 0-indexing, the row of A should be the node index minus one.
for i, edge in enumerate(E):
A[edge[0] - 1][i] = -1
A[edge[1] - 1][i] = 1 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])
p = xp.problem('network flow')
flow = p.addVariables(m) # flow variables declared on arcs
p.addConstraint(xp.Dot(A, flow) == - demand)
p.setObjective(xp.Dot(cost, flow))
p.optimize()
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.
import xpress as xp
a,b = 1,100
p = xp.problem()
x = p.addVariable(lb=-xp.infinity)
y = p.addVariable(lb=-xp.infinity)
p.setObjective((a-x)**2 + b*(y-x**2)**2)
p.controls.xslp_solver = 0 # solve it with SLP, not Knitro
solvestatus, solstatus = p.optimize()
print("solve status: ", solvestatus.name)
print("solution status: ", solstatus.name)
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.
p = xp.problem()
rho = [p.addVariable(name='rho_{}'.format(i), lb=1e-5, ub=1.0) for i in Vertices]
theta = [p.addVariable(name='theta_{}'.format(i), lb=-math.pi, ub=math.pi)
for i in Vertices]
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 2D array 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)
p = xp.problem()
# Create a 2D numpy array of (i,j) variables and link them to problem p
x = p.addVariables(N, N, vartype=xp.binary, name='q')
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.addConstraint(vertical, horizontal, diagonal1, diagonal2)
# Objective, to be maximized: number of queens on the chessboard
p.setObjective(xp.Sum(x), sense=xp.maximize)
p.optimize() 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)
p = xp.problem()
x = p.addVariables(N, N, N, vartype=xp.binary)
# 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.addConstraint(vertical, horizontal, subgrid, assign, init)
p.optimize()
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([[[m.addVariable(vartype=xp.binary)
for i in S1]
for j in S2]
for k in S3], dtype=xp.npvar)
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.optimize() 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.
p = xp.problem() x = p.addVariables(5) p.addConstraint(xp.Sum(x) >= 2) p.setObjective(xp.Sum(x[i]**2 for i in range(5))) p.optimize() 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 p = xp.problem() # Add a NumPy array of variables x = p.addVariables(5) 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.addConstraint(Lin_sys, Conv_c) p.setObjective(xp.Dot(x-x0, x-x0)) # minimize distance from x0 p.optimize()
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
p = xp.problem()
Q = np.arange(1, N**2 + 1).reshape(N, N)
x = p.addVariables(N)
x0 = np.random.random(N)
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.optimize('')
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 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.optimize()
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.
p = xp.problem()
x = p.addVariable()
y = p.addVariable()
cons1 = x + y >= 2
upperlim = 2*x + y <= 3
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 = [p.addVariable(vartype=xp.binary) for i in S]
# 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 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 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 |
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],
start=[0, 3, 9],
colind=[x[0],x[1],x[2], x[0],x[1],x[2],x[3],x[4], 5],
rowcoef=[1,2,3,4,5,6,7,8,-3],
names=['newcon1', 'newcon2'])
p.optimize()
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.optimize()
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.optimize()
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.optimize()
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.getCallbackSolution()
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 = prob.getCallbackSolution()
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)
p = xp.problem()
x = {(i, j): p.addVariable(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.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.timelimit = -20 # 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.optimize()
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 xixj is assigned to a new variable wij so that every occurrence of that product in the problem is replaced by wij. Assuming li and ui are the lower and upper bound on xi, respectively, we add the linear McCormick inequalities:
- wij ≥ lj xi + li xj - lj li
- wij ≥ uj xi + ui xj - uj ui
- wij ≤ lj xi + ui xj - lj ui
- wij ≤ uj xi + li xj - uj li
The bounds on the new auxiliary variable wij are a function of the bounds on xi and xj.
Below is the code that takes care of reformulating the problem. We first have to identify all terms xixj and create a dictionary linking each pair (i,j) to an auxiliary variable wij. 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(rowind=None, colind=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 wij, 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] = p.addVariable(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)
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 wij = xixj.
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 wij are satisfied, and if not it creates a branching object on either xi or xj. 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 wii = xi2, 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 wij = xi xj 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:
sol = prob.getCallbackSolution()
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):
try:
sol = prob.getCallbackSolution()
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 wij, possible tighter lower or upper bounds can be deduced on xi and/or xj.
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 for term x[i]^2: from new bounds on x[i],
# infer new bound for x[i]^2.
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.
Translated Mosel examples
The subdirectory modeling_examples of the Python examples directory contains a few examples from the Mosel distribution that were adapted to the Xpress Python interface:
- blend.py, blend2.py: variants of an oil blending optimization model;
- burglari.py, burglar.py, burglarl.py, burglar_rec.py: several variants of the knapsack problem
- catenary.py: optimization model for finding the position of all elements of a hanging chain
- chess.py, chess2.py: two variants on the simple problem of production management;
- coco.py: Multiperiod production planning problem;
- complex_test.py: an example of complex numbers (a native type in Python
- fstns.py: the problem of firestation siting;
- date_test.py: an example of dates using the datatime module;
- pplan.py: a production planning example;
- trans.py: a transportation problem.
© 2001-2025 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.