Linear Programming
Introduction
Linear programming problems are optimization problems where the objective function and constraints are all linear.
Mathematica has a collection of algorithms for solving linear optimization problems with real variables, accessed via
LinearProgramming,
FindMinimum,
FindMaximum,
NMinimize,
NMaximize,
Minimize, and
Maximize.
LinearProgramming gives direct access to linear programming algorithms, provides the most flexibility for specifying the methods used, and is the most efficient for largescale problems.
FindMinimum,
FindMaximum,
NMinimize,
NMaximize,
Minimize, and
Maximize are convenient for solving linear programming problems in equation and inequality form.
This solves a linear programming problem
Min x + 2 y
s.t. 5 x + y = 7
x + y ≥ 26
x ≥ 3, y ≥ 4
using Minimize.
Out[1]=  

The LinearProgramming Function
LinearProgramming is the main function for linear programming with the most flexibility for specifying the methods used, and is the most efficient for largescale problems.
The following options are accepted.
Options for LinearProgramming.
The
Method option specifies the algorithm used to solve the linear programming problem. Possible values are
Automatic,
"Simplex",
"RevisedSimplex", and
"InteriorPoint". The default is
Automatic, which automatically chooses from the other methods based on the problem size and precision.
The
Tolerance option specifies the convergence tolerance.
Examples
Difference between Interior Point and Simplex and/or Revised Simplex
The simplex and revised simplex algorithms solve a linear programming problem by moving along the edges of the polytope defined by the constraints, from vertices to vertices with successively smaller values of the objective function, until the minimum is reached.
Mathematica's implementation of these algorithm uses dense linear algebra. A unique feature of the implementation is that it is possible to solve exact/extended precision problems. Therefore these methods are suitable for smallsized problems for which nonmachinenumber results are needed, or a solution on the vertex is desirable.
Interior point algorithms for linear programming, loosely speaking, iterate from the interior of the polytope defined by the constraints. They get closer to the solution very quickly, but unlike the simplex/revised simplex algorithms, do not find the solution exactly.
Mathematica's implementation of an interior point algorithm uses machine precision sparse linear algebra. Therefore for largescale machineprecision linear programming problems, the interior point method is more efficient and should be used.
This solves a linear programming problem that has multiple solutions (any point that lies on the line segment between {1, 0} and {1, 0} is a solution); the interior point algorithm gives a solution that lies in the middle of the solution set.
Out[6]=  

Using Simplex or RevisedSimplex, a solution at the boundary of the solution set is given.
Out[7]=  

Finding Dual Variables
Given the general linear programming problem
It is useful to know solutions for both for some applications.
The relationship between the solutions of the primal and dual problems is given by the following table.
 
feasible  feasible 
unbounded  infeasible or unbounded 
infeasible  unbounded or infeasible 
When both problems are feasible, then the optimal values of (P) and (D) are the same, and the following complementary conditions hold for the primal solution
x, and dual solution
y, z.
DualLinearProgramming returns a list
{x, y, z, w}.
This solves the primal problem
Min 3 x_{1}+4 x_{2 }
s.t. x_{1}+2 x_{2}≥ 5
1≤x_{1}≤4, 1≤x_{2}≤4,
as well as the dual problem
Max 5y_{1}+ z_{1}+z_{2}  4 w_{1}  4 w_{2 }
s.t. y_{1}+ z_{1} w_{1}= 3
2 y_{1} + z_{2}  w_{2} = 4
y_{1}, z_{1}, z_{2}, w_{1}, w_{2}≥ 0
Out[14]=  

This confirms that the primal and dual give the same objective value.
Out[15]=  
Out[16]=  

The dual of the constraint is y = {2.}, which means that for one unit of increase in the righthand side of the constraint, there will be 2 units of increase in the objective. This can be confirmed by perturbing the righthand side of the constraint by 0.001.
Out[17]=  

Indeed the objective value increases by twice that amount.
Out[18]=  

Dealing with Infeasibility and Unboundedness in the Interior Point Method
The primaldual interior point method is designed for feasible problems; for infeasible/unbounded problems it will diverge, and with the current implementation, it is difficult to find out if the divergence is due to infeasibility, or unboundedness.
A heuristic catches infeasible/unbounded problems and issues a suitable message.
Out[19]=  

Sometimes the heuristic cannot tell with certainty if a problem is infeasible or unbounded.
Out[20]=  

Using the Simplex method as suggested by the message shows that the problem is unbounded.
Out[21]=  

The Method Options of "InteriorPoint"
"TreatDenseColumn" is a method option of
"InteriorPoint" that decides if dense columns are to be treated separately. Dense columns are columns of the constraint matrix that have many nonzero elements. By default, this method option has the value
Automatic, and dense columns are treated separately.
Large problems that contain dense columns typically benefit from dense column treatment.
Out[98]=  
Out[99]=  

Importing Large Datasets and Solving LargeScale Problems
A commonly used format for documenting linear programming problems is the Mathematical Programming System (MPS) format. This is a text format consisting of a number of sections.
Importing MPS Formatted Files in Equation Form
Mathematica is able to import
MPS formatted files. By default,
Import of
MPS data returns a linear programming problem in equation form, which can then be solved using
Minimize or
NMinimize.
This solves the linear programming problem specified by MPS file "afiro.mps".
Out[25]=  
Out[26]=  

LargeScale Problems: Importing in Matrix and Vector Form
For largescale problems, it is more efficient to import the
MPS data file and represent the linear programming using matrices and vectors, then solve using
LinearProgramming.
This shows that for MPS formatted data, the following 3 elements can be imported.
Out[101]=  

This imports the problem "ganges", with 1309 constraints and 1681 variables, in a form suitable for LinearProgramming. 
This solves the problem and finds the optimal value.
Out[104]=  

The "ConstraintMatrix" specification can be used to get the sparse constraint matrix only.
Out[105]=  

Free Formatted MPS Files
Standard
MPS formatted files use a fixed format, where each field occupies a strictly fixed character position. However some modeling systems output
MPS files with a free format, where fields are positioned freely. For such files, the option
"FreeFormat">True can be specified for
Import.
This string describes an MPS file in free format. 
This gets a temporary file name, and exports the string to the file. 
This imports the file, using the "FreeFormat">True option.
Out[126]=  

Linear Programming Test Problems
Through the
ExampleData function, all
NetLib linear programming test problems can be accessed.
This finds all problems in the Netlib set.
Out[12]=  

This shows other properties that can be imported for the "afiro" problem.
Out[10]=  

This imports "afiro" in equation form.
Out[11]=  

Application Examples of Linear Programming
L1Norm Minimization
It is possible to solve an
l_{1} minimization problem
by turning the system into a linear programming problem
This defines a function for solving an l_{1} minimization problem. 
The following is an overdetermined linear system.
Design of an Optimal Anchor
The example is adopted from
[2]. The aim is to design an anchor that uses as little material as possible to support a load.
This problem can be modeled by discretizing and simulating it using nodes and links. The modeling process is illustrated using the following figure. Here a grid of 7×10 nodes is generated. Each node is then connected by a link to all other nodes that are of Manhattan distance of less than or equal to three. The three red nodes are assumed to be fixed to the wall, while on all other nodes, compression and tension forces must balance.
Each link represents a rigid rod that has a thickness, with its weight proportional to the force on it and its length. The aim is to minimize the total material used, which is
Hence mathematically this is a linearly constrained minimization problem, with objective function a sum of absolute values of linear functions.
The absolute values
force*link_length in the objective function can be replaced by breaking down
force into a combination of compression and tension forces, with each nonnegative. Thus assume
E is the set of links,
V the set of nodes,
l_{ij} the length of link between nodes
i and
j,
c_{ij}and
t_{ij} the compression and tension forces on the link; then the above model can be converted to a linear programming problem
The following sets up the model, solves it, and plots the result; it is based on an AMPL model [2]. 
This solves the problem by placing 30 nodes in the horizontal and vertical directions.
Out[6]=  

If, however, the anchor is fixed not on the wall, but on some points in space, notice how the results resemble the shape of some leaves. Perhaps the structure of leaves is optimized in the process of evolution.
Out[11]=  

Algorithms for Linear Programming
Simplex and Revised Simplex Algorithms
The simplex and revised simplex algorithms solve linear programming problems by constructing a feasible solution at a vertex of the polytope defined by the constraints, and then moving along the edges of the polytope to vertices with successively smaller values of the objective function until the minimum is reached.
Although the sparse implementation of simplex and revised algorithms are quite efficient in practice, and are guaranteed to find the global optimum, they have a poor worstcase behavior: it is possible to construct a linear programming problem for which the simplex or revised simplex method takes a number of steps exponential in the problem size.
Mathematica implements simplex and revised simplex algorithms using dense linear algebra. The unique feature of this implementation is that it is possible to solve exact/extended precision problems. Therefore these methods are more suitable for smallsized problems for which nonmachine number results are needed.
This sets up a random linear programming problem with 20 constraints and 200 variables. 
This solves the problem. Typically, for a linear programming problem with many more variables than constraints, the revised simplex algorithm is faster. On the other hand, if there are many more constraints than variables, the simplex algorithm is faster. 
If only machinenumber results are desired, then the problem should be converted to machine numbers, and the interior point algorithm should be used. 
Interior Point Algorithm
Although the simplex and revised simplex algorithms can be quite efficient on average, they have a poor worstcase behavior. It is possible to construct a linear programming problem for which the simplex or revised simplex methods take a number of steps exponential in the problem size. The interior point algorithm, however, has been proven to converge in a number of steps that are polynomial in the problem size.
Furthermore, the
Mathematica simplex and revised simplex implementation use dense linear algebra, while its interior point implementation uses machinenumber sparse linear algebra. Therefore for largescale, machinenumber linear programming problems, the interior point method is more efficient and should be used.
Interior Point Formulation
Consider the standardized linear programming problem
where
c, xR^{n},
AR^{m×n},
bR^{m}. This problem can be solved using a barrier function formulation to deal with the positive constraints
The firstorder necessary condition for the above problem gives
Let
X denote the diagonal matrix made of the vector
x, and
z=t X^{1} e.
This is a set of
2m+n linear/nonlinear equations with constraints. It can be solved using Newton's method
One way to solve this linear system is to use Gaussian elimination to simplify the matrix into block triangular form.
To solve this block triangular matrix, the socalled normal system can be solved, with the matrix in this normal system
This matrix is positive definite, but becomes very illconditioned as the solution is approached. Thus numerical techniques are used to stabilize the solution process, and typically the interior point method can only be expected to solve the problem to a tolerance of about
, with tolerance explained below.
Mathematica uses Mehrotra's predictorcorrector scheme [1].
Convergence Tolerance
General Linear Programming problems are first converted to the standard form
with the corresponding dual
The convergence criterion for the interior point algorithm is
with the tolerance set, by default, to
.
References
[1] Vanderbei R. Linear Programming: Foundations and Extensions. SpringerVerlag (2001)
[2] Mehrotra S. "On the Implementation of a PrimalDual Interior Point Method" SIAM Journal on Optimization 2 (1992): 575601.