Linear Programming
The goal of linear programming is to optimize a linear function subject to linear constraints. In its most general form, a linear program is an optimization problem
Maximize | Σ ci xi |
|---|---|
subject to | lhsi ≤ Aix ≤ rhsi |
li ≤ xi ≤ ui |
A typical linear program can have many thousands of variables and several thousand constraints.
Variables
Linear programs use a specialized kind of decision variable, LinearProgramVariable. Since the objective function is linear, it is specified fully by the coefficient of each variable, usually called the cost. It can be set or retrieved through the Cost property.
By convention, linear program variables have a lower bound of zero by default.
Linear programs may have variables that are restricted to be integers. Integer variables can be created using the model's AddIntegerVariable method. If the variable can only take on the values 0 or 1, the AddBinaryVariable method provides a shortcut. Existing variables can be restricted to integer values by setting the IsInteger property to true.
The LinearConstraint class is used to represent constraints in a linear program. Every constraint must have a unique name, which can be accessed through the Name property.
A constraint is a linear combination of variables that is required to lie between specified boundaries. The LowerBound and UpperBound properties define these boundaries. When a constraint is not bounded from above, its upper bound is Double.PositiveInfinity. When a constraint is not bounded from below, its lower bound is Double.NegativeInfinity.
Variables and constraints can be accessed through the linear program's Variables and Constraints collections. These collections are indexed by position and by name.
Defining linear programs
A linear program is represented by the LinearProgram class. There are three ways to construct LinearProgram objects.
Defining linear programs using constructors
The quickest way to define a linear program is by passing the vectors and matrices that appear in the definition above to a constructor. There are three constructors that can be used in this way.
The first constructor takes three arguments. It defines a linear program in standard form. This is a linear program where variables are constrained to be positive, and all constraints are inequalities with an upper bound. The first parameter is a Vector<T> containing the coefficients of the objective function. It corresponds to the vector c in the definition. The second argument is a Matrix<T> that contains the coefficients of the inequalities. It corresponds to the matrix A in the definition. The third argument is a Vector<T> containing the right-hand-sides of the inequalities. The code below creates the following a linear program:
s.t. | x1 + 2x2 ≤ 5 |
5x1 + 3x2 ≤ 11 |
Vector<double> c = Vector.Create(1.0, 1.0);
Matrix<double> A = Matrix.CopyFrom(
new [] { 1.0, 5.0, 2.0, 3.0 },
2, 2, MatrixElementOrder.RowMajor);
Vector<double> b = Vector.Create(5.0, 11.0);
LinearProgram lp1 = new LinearProgram(c, A, b);The second constructor adds a fourth parameter. If present, this parameter indicates the number of equality constraints. The equality constraints appear above the inequality constraints.
The third constructor is the most flexible. It has six parameters in all. The first two are once again the cost coefficients and the coefficient matrix. The third and fourth parameters are vectors that contain the lower and upper bounds of the constraints. Each type of constraint can be expressed using a lower and upper bound, as shown in the table below.
Constraint | Lower bound | Upper bound |
|---|---|---|
c ≤ rhs | -infinity | rhs |
c ≥ rhs | rhs | +infinity |
c = rhs | rhs | rhs |
lhs ≤ c ≤ rhs | lhs | rhs |
The fifth and sixth parameters are vectors that contain the lower and upper bound for the variables. The different types of constraints on variables can be expressed the same way as for constraints. The following code sample builds the following linear program:
s.t. | 0.5 ≤ x1 + x2 ≤ 1.5 |
0 ≤ x1 ≤ 1 | |
0 ≤ x2 ≤ 1 |
Vector<double> c = Vector.Create(-1.0, -3.0);
Matrix<double> A = Matrix.CopyFrom(new double[] { 1.0, 1.0 },
1, 2, MatrixElementOrder.RowMajor);
Vector<double> lc = Vector.Create(0.5);
Vector<double> uc = Vector.Create(1.5);
Vector<double> lx = Vector.Create(0.0, 0.0);
Vector<double> ux = Vector.Create(1.0, 1.0);
LinearProgram lp2 = new LinearProgram(c, A, lc, uc, lx, ux);Building linear programs from scratch
The LinearProgram class also has a default constructor that creates an empty linear program. You can then call the linear program's AddVariable and AddLinearConstraint methods.
The AddVariable method has two overloads. The first takes two arguments. The first is the name of the variable. The second is the cost associated with the variable. The second overload has two additional parameters that specify the lower and upper bound of the variable. If no values are specified, a lower bound of zero and an upper bound of infinity are assumed.
The AddLinearConstraint method has four overloads. The first takes four parameters. The first is the name of the constraint. The second is a Vector<T> that contains the coefficients of the variables. The third parameter is a ConstraintType value that specifies the type of constraint. The possible values are listed in the table below. The last parameter is the right-hand side of the constraint.
Value | Description |
|---|---|
Equal | c = rhs |
GreaterThanOrEqual | c ≥ rhs |
LessThanOrEqual | c ≤ rhs |
The second overload also takes four parameters. The first two are the same as before. The third and fourth parameters are the lower and upper bound of the constraint. The different types of constraints can be specified according to the examples in Table 1.
The third and fourth overloads are similar to the first two. The only difference is that the coefficients are specified using a Double array instead of a vector.
The following example creates a linear program that is identical to the last example:
LinearProgram lp3 = new LinearProgram(c, A, lc, uc, lx, ux);
lp3.AddVariable("X1", -1, 0.0, 1.0);
lp3.AddVariable("X2", -3, 0.0, 1.0);
lp3.AddLinearConstraint("C1", new double[] { 1.0, 1.0 }, 0.5, 1.5);Although not strictly necessary, it is advisable to create the variables first, followed by the constraints. If the vector containing the coefficients for a constraint is longer than the number of variables that have been defined, additional variables are added automatically. The kth variable is given the name x<g>.
Whenever possible, constraints on variables should be preferred over explicit linear program constraints. The time to solve a linear program depends most heavily on the number of constraints.
Importing linear programs: the MPS format
Use the MpsReader class to read a LinearProgram from a file in MPS format. The format was named after an early linear programming system from IBM. It has since become a de facto standard for defining linear programs in text format in commercial linear programming systems. There are many online resources that describe the format in detail.
The MpsReader class has one static method, Read, which is overloaded. The first method takes a string containing the path to an MPS file and returns a LinearProgram object that represents the linear program in the file. The second method does the same but gets its input from a System.IO.StreamReader.
Several extensions of the MPS format have been defined over the years. The current implementation of MpsReader doesn't support these extensions. A System.FormatException is thrown when an unrecognized extension is found.
Default Method and Model-Specific Options
Continuous linear programs use the revised dual simplex method with LUSOL factorization and Bartels-Golub updates. The updates reuse most of the current factorization after a pivot. The default also combines useful bound flips within a step. When eight consecutive searches find no useful flips, it uses the Harris ratio test for 64 selections before trying again.
Presolve removes unnecessary constraints, substitutes suitable variables and combines exact duplicates before solving. All seven additional reductions run by default in the existing fixed-point loop. Scaling then applies to the reduced model. The returned solution, objective and dual variables are reconstructed in the original model's units and order.
// Minimize -3*x-2*y, with x+y <= 4, 0 <= x <= 2 and 0 <= y <= 3.
LinearProgram model = new LinearProgram();
model.AddVariable("x", -3, 0, 2);
model.AddVariable("y", -2, 0, 3);
model.AddLinearConstraint("capacity", Vector.Create(1.0, 1), ConstraintType.LessThanOrEqual, 4);
Vector<double> solution = model.Solve();
// The default uses LUSOL, adaptive bound flipping and all seven reductions.
if (model.Status != OptimizationModelStatus.Optimal) throw new InvalidOperationException("LP did not converge.");
// The solution is (2,2), with objective -10.Use SimplexOptions to change a supported setting. The defaults are a useful starting point; a smaller reduced problem or fewer iterations does not always mean a faster solve. Compare alternatives on representative models from your application.
For characterized large sparse network models, PartialSteepestEdge examines a random-start half of the current candidate rows, with at least 32 candidates. It performs a full scan before concluding that no useful candidate remains. This can reduce pricing work, but can increase iterations on degenerate or poorly conditioned models. Keep full SteepestEdge pricing for those models. The following sparse example shows explicit partial pricing; it does not install an automatic model classifier.
const int n = 1000;
int[] rows = new int[2 * n], columns = new int[2 * n];
double[] entries = new double[2 * n];
for (int j = 0; j < n; j++)
{
rows[2 * j] = j; columns[2 * j] = j; entries[2 * j] = 2;
rows[2 * j + 1] = (j + 1) % n; columns[2 * j + 1] = j; entries[2 * j + 1] = 1;
}
// Minimize 3*sum(x), with 2*x[i]+x[i-1] >= 3 and x >= 0.
LinearProgram model = new LinearProgram(Vector.FromFunction(n, i => 3.0),
Matrix.CreateSparse(n, n, rows, columns, entries), Vector.FromFunction(n, i => 3.0),
Vector.FromFunction(n, i => double.PositiveInfinity), Vector.Zeros<double>(n),
Vector.FromFunction(n, i => double.PositiveInfinity));
SimplexOptions options = new SimplexOptions
{ UsePartialPricing = true, RandomSeed = 12345 };
// Partial pricing is an explicit choice for a characterized sparse family.
// Compare it with the default full pricing on your own models.
Vector<double> solution = model.Solve(options, null);
if (model.Status != OptimizationModelStatus.Optimal) throw new InvalidOperationException("LP did not converge.");
// Every solution component is 1. Summing the constraints proves the minimum is 3*n.PresolveReductions is a set of flags. To reduce allocation on a model family where column merging is costly, disable DuplicateColumns while keeping the other reductions. This can require more pivots. None retains the elementary presolve rules. BasisUpdateMethod can explicitly select EtaUpdates, and setting UseBoundFlippingRatioTest to false selects Harris without bound flipping. UsePartialPricing enables partial scanning with dual steepest-edge weights. The BFRT cooldown iteration counts remain configurable, with 8/64 recommended as the general starting point. These settings choose supported methods, without enabling the lab's perturbation or refinement experiments.
Options are copied at the beginning of a solve. You can reuse and edit the supplied options between solves. RandomSeed makes the perturbation and independent partial-pricing streams reproducible; equal-score choices in parallel full pricing can still vary. Automatic integer programs retain their existing sparse-LU, Harris and elementary-presolve policy. Explicit LUSOL, partial-pricing and column-indexing choices apply to continuous LPs.
Solving linear programs
Once a linear program has been defined, solving it is very straightforward. The Solve method does all the work.
The revised dual simplex method is used. Running time depends on the number of variables, constraints and nonzero coefficients, and on the number of pivots needed. The options above let you choose a supported method for your model family.
The Solve method returns a Vector<T> containing the optimal solution. You should inspect the Status property to make sure an optimal solution was indeed found. This property is of type OptimizationModelStatus and can take on the following values:
Value | Description |
|---|---|
Unknown | Nothing is known about the solution. This is the status returned before the SolveModel method is called. |
Optimal | The linear program was solved successfully and the solution that was returned is the optimal solution. |
Infeasible | The linear program does not have a solution because some of the constraints conflict with each other. |
Unbounded | The linear program does not have a finite solution. The cost function can be made arbitrarily small. |
Two other status values, Feasible and DualFeasible are used while solving the linear program. Finally, the OptimalValue property returns the value of the cost function at the solution. The GetDualSolution method returns a Vector<T> containing the solution to the dual problem.
The code below solves the linear program defined earlier:
Vector<double> x = lp1.Solve();
Vector<double> y = lp1.GetDualSolution();
Console.WriteLine("Status: {0}", lp1.Status);
Console.WriteLine("Primal: {0:F1}", x);
Console.WriteLine("Dual: {0:F1}", y);
Console.WriteLine("Optimal value: {0:F1}", lp1.OptimalValue);