linearProgramming

Solves a linear programming problem.

Synopsis

linearProgramming (a, b, c)

Required Arguments

float a[[]] (Input)
Array of size m × n containing a matrix with coefficients of the m constraints.
float b[] (Input)
Array with m components containing the right-hand side of the constraints; if there are limits on both sides of the constraints, then b contains the lower limit of the constraints.
float c[] (Input)
Array with n components containing the coefficients of the objective function.

Return Value

The solution x of the linear programming problem. If no solution can be computed, then None is returned.

Optional Arguments

upperLimit, double[] (Input)
Array with m components containing the upper limit of the constraints that have both the lower and the upper bounds. If no such constraint exists, then upperLimit is not needed.
constrType, int[] (Input)

Array with m components indicating the types of general constraints in the matrix a. Let \(r_i=a_{i1} x_1+\ldots+a_{in} x_n\). Then, the value of constrType[i] signifies the following:

constrType[i] Constraint
0 \(r_i=b_i\)
1 \(r_i\leq upperLimit_i\)
2 \(r_i\geq b_i\)
3 \(b_i\leq r_i\leq upperLimit_i\)
4 Ignore this constraint

Default: constrType = 0

lowerBound, double[] (Input)

Array with n components containing the lower bound on the variables. If there is no lower bound on a variable, then \(10^{30}\) should be set as the lower bound.

Default: lowerBound = 0

upperBound, double[] (Input)

Array with n components containing the upper bound on the variables. If there is no upper bound on a variable, then \(-10^{30}\) should be set as the upper bound.

Default: no upper bound

refinement (Input)

The coefficient matrices and other data are saved at the beginning of the computation. When finished this data together with the solution obtained is checked for consistency. If the discrepancy is too large, the solution process is restarted using the problem data just after processing the equalities, but with the final x values and final active set.

Default: Refinement is not performed.

extendedRefinement (Input)

This is similar to refinement, except it iterates until there is a sign that no further progress is possible (recommended if all the accuracy possible is desired).

Default: Extended refinement is not performed.

obj (Output)
Optimal value of the objective function.
dual (Output)
An array with m components containing the dual solution. On return, the necessary space is allocated by linearProgramming.

Description

The function linearProgramming uses an active set strategy to solve linear programming problems, i.e., problems of the form

\[\begin{split}\begin{array}{llc} \min\limits_{x \in R^n} c^T x & \text{ subject to } & b_l \leq A_x \leq b_u \\ & & x_l \leq x \leq x_u \\ \end{array}\end{split}\]

where c is the objective coefficient vector, A is the coefficient matrix, and the vectors \(b_l\), \(b_u\), \(x_l\), and \(x_u\) are the lower and upper bounds on the constraints and the variables, respectively.

Examples

Example 1

The linear programming problem in the standard form

\[\begin{split}\text{subject to } \begin{array}{l} \min f(x) = -x_1 -3_{x_2} \\ \begin{array}{lrrrrr} x_1 & +x_2 & +x_3 & & & & =1.5 \\ x_1 & +x_2 & & -x_4 & & & =0.5\\ x_1 & & & & +x_5 & & =1.0 \\ & x_2 & & & & +x_6 & =1.0 \\ \end{array} \\ x_i \geq 0, \text{ for } i = 1, \ldots, 6 \end{array}\end{split}\]

is solved.

from numpy import *
from pyimsl.math.linearProgramming import linearProgramming
from pyimsl.math.writeMatrix import writeMatrix


m = 4
n = 6
a = array([[1.0, 1.0, 1.0, 0.0, 0.0, 0.0],
           [1.0, 1.0, 0.0, -1.0, 0.0, 0.0],
           [1.0, 0.0, 0.0, 0.0, 1.0, 0.0],
           [0.0, 1.0, 0.0, 0.0, 0.0, 1.0]])
b = (1.5, 0.5, 1.0, 1.0)
c = (-1.0, -3.0, 0.0, 0.0, 0.0, 0.0)

# Solve the LP problem
x = linearProgramming(a, b, c)

# Print x
writeMatrix("x", x)

Output

 
                                      x
          1            2            3            4            5            6
        0.5          1.0          0.0          1.0          0.5          0.0

Example 2

This example demonstrates how the function readMps can be used together with linearProgramming to solve a linear programming problem defined in an MPS file. The MPS file used in this example is an uncompressed version of the file ‘afiro’, available from http://www.netlib.org/lp/data/. This example also demonstrates the use of the optional argument refinement to activate iterative refinement in linearProgramming.

from __future__ import print_function
import os
from numpy import *
from pyimsl.math.linearProgramming import linearProgramming
from pyimsl.math.readMps import readMps
from pyimsl.math.writeMatrix import writeMatrix


def A(i, j):
    return a[(i) * problem.ncolumns + j]


# Read the MPS file
filename = os.getcwd() + os.sep + "afiro.mps"
problem = readMps(filename)

# Setup constraint type array
irtype = empty(problem[0].nrows)
for i in range(0, problem[0].nrows):
    irtype[i] = 3

# Setup the constraint matrix
a = zeros((problem[0].nrows, problem[0].ncolumns), dtype=double)
for k in range(0, problem[0].nonzeros):
    i = problem[0].constraint[k].row
    j = problem[0].constraint[k].col
    a[i, j] = problem[0].constraint[k].val

# Setup constraint bounds
bl = empty((problem[0].nrows), dtype=double)
bu = empty((problem[0].nrows), dtype=double)
for i in range(0, problem[0].nrows):
    bl[i] = problem[0].lower_range[i]
    bu[i] = problem[0].upper_range[i]

# Setup variable bounds.  Be sure to account for
# how unbounded variables should be set.
xlb = empty((problem[0].ncolumns), dtype=double)
xub = empty((problem[0].ncolumns), dtype=double)
for i in range(0, problem[0].ncolumns):
    if (problem[0].lower_bound[i] == problem[0].negative_infinity):
        xlb[i] = 1.0e30
    else:
        xlb[i] = problem[0].lower_bound[i]
    if (problem[0].upper_bound[i] == problem[0].positive_infinity):
        xub[i] = -1.0e30
    else:
        xub[i] = problem[0].upper_bound[i]

inputObjective = empty((problem[0].ncolumns), dtype=double)
for i in range(0, problem[0].ncolumns):
    inputObjective[i] = problem[0].objective[i]

# Solve the LP problem
objective = []
x = linearProgramming(a, bl, inputObjective,
                      upperLimit=bu,
                      constrType=irtype,
                      lowerBound=xlb,
                      upperBound=xub,
                      refinement=True,
                      obj=objective)

# Output results
print("Problem Name: ", problem[0].name.decode("utf-8"))
print("Objective: %e" % (objective[0]))
writeMatrix("Solution", x, column=True)

Output

Problem Name:  AFIRO    
Objective: -4.647531e+02
 
   Solution
 1         80.0
 2         25.5
 3         54.5
 4         84.8
 5         57.9
 6          0.0
 7          0.0
 8          0.0
 9          0.0
10          0.0
11          0.0
12          0.0
13         18.2
14         39.7
15         61.3
16        500.0
17        475.9
18         24.1
19          0.0
20        215.0
21        363.9
22          0.0
23          0.0
24          0.0
25          0.0
26          0.0
27          0.0
28          0.0
29        339.9
30         20.1
31        156.5
32          0.0

Note Errors

IMSL_MULTIPLE_SOLUTIONS Multiple solutions giving essentially the same minimum exist.

Warning Errors

IMSL_SOME_CONSTRAINTS_DISCARDED Some constraints were ignored or discarded because they were too linearly dependent on other active constraints.
IMSL_ALL_CONSTR_NOT_SATISFIED All constraints are not satisfied. If a feasible solution is possible then try using refinement by supplying optional argument refinement.
IMSL_CYCLING_OCCURRING The algorithm appears to be cycling. Using refinement may help.

Fatal Errors

IMSL_PROB_UNBOUNDED The problem is unbounded.
IMSL_PIVOT_NOT_FOUND An acceptable pivot could not be found.