cabeggar

Linear Programming with CPLEX

Posted on |

If we are trying to solve optimization problems, sometimes we will end up with a mathematical programming which requires modelling the problem and solve it with some solvers. For example, when we are computing VM placement in data centers, we can use this kind of optimization to minimize usage of computer resources and network utilization. IBM ILOG CPLEX Optimizer is a powerful solver that can help us easily solve this kind of problems.

CPLEX includes CPLEX studio, an interactive command line tool and APIs for different languages like C++, Java and Python. In this blog, I’ll go through CPLEX’s example for Python API to solve a linear programming problem, which is a type of mathematical problems.

Problem Statement

The problem we want to solve is as follows.

1
2
3
4
5
6
Maximize			x1 + 2*x2 + 3*x3
subject to			-x1 + x2 + x3 <= 20
					x1 - 3*x2 + x3 <= 30
with these bounds	0 <= x1 <= 40
					0 <= x2 <= infinity
					0 <= x3 <= infinity

Prameters

For the previous examples, some global variables can be written as follows:

1
2
3
4
5
6
7
8
9
10
11
import cplex
from cplex.exceptions import CplexError
import sys
# data common to all populatedby functions
my_obj = [1.0, 2.0, 3.0]
my_ub = [40.0, cplex.infinity, cplex.infinity]
my_colnames = ["x1", "x2", "x3"]
my_rhs = [20.0, 30.0]
my_rownames = ["c1", "c2"]
my_sense = "LL" # "L" for "less-than", 2 "L" for 2 contraints

Populate by rows

The most common way to model this problems is to set constraint as rows.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def populatebyrow(prob):
    prob.objective.set_sense(prob.objective.sense.maximize)
    # since lower bounds are all 0.0 (the default), lb is omitted here
    prob.variables.add(obj = my_obj, ub = my_ub, names = my_colnames)
    # can query variables like the following bounds and names:
    # lbs is a list of all the lower bounds
    lbs = prob.variables.get_lower_bounds()
    # ub1 is just the first lower bound
    ub1 = prob.variables.get_upper_bounds(0)
    # name is ["x1", "x3"]
    names = prob.variables.get_names([0, 2])
    rows = [[[0, 1, 2], [-1.0, 1.0, 1.0]],
            [[0, 1, 2], [1.0, -3.0, 1,0]]]
    # We can also write as this
    # rows = [[[0, "x2", "x3"], [-1.0, 1.0, 1.0]],
    #         [["x1", 1, 2], [1.0, -3.0, 1,0]]]
    prob.linear_constraints.add(lin_expr = rows,
                                senses = my_sense,
                                rhs = my_rhs,
                                names = my_rownames)
    # because there are two arguments, they are taken to specify a range
    # thus, cols is the entire constraint matrix as a list of column vectors
    cols = prob.variables.get_cols("x1", "x3")

We first set our target to maximize our objective. Then we add our objectives, upper bound of variables (lower bound set to 0 as default) and name of columns (variables). Then we set the constrainst by rows, first adding indexes or name of variables, and set coefficient secondly. At last we add our constraints into the problem solver.

Populate by columns

In some situations, it may be better to populate by columns.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def populatebycolumn(prob):
    prob.objective.set_sense(prob.objective.sense.maximize)
    prob.linear_constraints.add(rhs = my_rhs,
                                senses = my_sense,
                                names = my_rownames)
    """
    populate by columns
    constraints: -x1+x2+x3<=20
                 x1-3x2+x3<=30
    converted into: -c1+c2
                    c1-3c2
                    c1+c2
    """
    c = [[[0, 1], [-1.0, 1.0]],
         [["c1", 1], [1.0, -3.0]],
         [[0, "c2"], [1.0, 1.0]]]
    prob.variables.add(obj = my_obj,
                       ub = my_ub,
                       names = my_colnames,
                       columns = c)

It’s like we tranverse the matrix in previous method and we set “c” as basic element in constraints. We also add constraints before adding variables in this method. Whether to populate by rows or columns should be based on number of constraints and number of variables.

Populate by non-zero

At last, there’s a third method to populate the problem with non-zero elements in the matrix.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def populatebynonzero(prob):
    prob.objective.set_sense(prob.objective.sense.maximize)
    prob.linear_constraints.add(rhs = my_rhs,
                                senses = my_sense,
                                names = my_rownames)
    prob.variables.add(obj = my_obj,
                       ub = my_ub,
                       names = my_colnames)
    """
    The following is a sparse matrix
    number in rows and cols represents location in matrix
    vals represent values at specific location
    """
    rows = [0, 0, 0, 1, 1, 1]
    cols = [0, 1, 2, 0, 1, 2]
    vals = [-1.0, 1.0, 1.0, 1.0, -3.0, 1.0]
    prob.linear_constraints.set_coefficients(zip(rows, cols, vals))
    # can also change one coefficient at a time
    # prob.linear_constraints.set_coefficients(1, 1, -3.0)
    # or pass in a list of triples
    # prob.linear_constraints.set_coefficients([(0, 1, 1.0), (1, 1, -3.0)])

This method is good at handle sparse matrix. In this example, the matrix is like this:

1
2
-1	1	1
1	-3	1

The value in rows and cols composite the location of a value in corresponding order in vals.

Final

There are still some final codes which is easy to understand.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def lpex1(pop_method):
    try:
        my_prob = cplex.Cplex()
        if pop_method == "r":
            handle = populatebyrow(my_prob)
        if pop_method == "c":
            handle = populatebycolumn(my_prob)
        if pop_method == "n":
            handle = populatebynonzero(my_prob)
        my_prob.solve()
    except CplexError, exc:
        print exc
        return
    numrows = my_prob.linear_constraints.get_num()
    numcols = my_prob.variables.get_num()
    print
    # solution.get_status() returns an integer code
    print "Solution status = ", my_prob.solution.get_status(), ":",
    # the following line prints the corresponding string
    print my_prob.solution.status[my_prob.solution.get_status()]
    print "Solution value = ", my_prob.solution.get_objective_value()
    slack = my_prob.solution.get_linear_slacks()
    pi = my_prob.solution.get_dual_values()
    x = my_prob.solution.get_values()
    dj = my_prob.solution.get_reduced_costs()
    for i in range(numrows):
        print "Row %d:  Slack = %10f    Pi = %10f" % (i, slack[i], pi[i])
    for j in range(numcols):
        print "Column %d:   Value = %10f    Reduced cost = %10f" % (j, x[j], dj[j])
    my_prob.write("lpex1.lp")
if __name__ == "__main__":
    if len(sys.argv) != 2 or sys.argv[1] not in ["-r", "-c", "-n"]:
        print "Usage: lpex1.py -X"
        print "   where X is one of the following options"
        print "      r          generate problem by row"
        print "      c          generate problem by column"
        print "      n          generate problem by nonzero"
        print " Exiting..."
        sys.exit(-1)
    lpex1(sys.argv[1][1])
else:
    prompt = """Enter the letter indicating how the problem data should be populated:
    r : populate by rows
    c : populate by columns
    n : populate by nonzeros\n ? > """
    r = 'r'
    c = 'c'
    n = 'n'
    lpex1(input(prompt))