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)² + b (y-x²)².

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 rho_i from the center, and at angle theta_i. 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 * rho_i * rho_j * sin (theta_i - theta_j). 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 q².

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('')