Simplex method

The simplex method is the standard algorithm to solve LP problems and is based on the following observations:

• An LP problem is a convex problem; hence, if we find a locally optimal solution, we have also found a global optimizer.
• An LP problem is a concave problem, too; hence, we know that we may restrict the search for the optimal solution to the boundary of the feasible set, which is a polyhedron. Actually, it can be shown that there is an optimal solution corresponding to a vertex, or extreme point, of the polyhedron.
• So, we need a clever way to explore extreme points of the feasible set. Geometrically, we may imagine moving from a vertex to a neighboring one, trying to improve the objective function. When there is no neighboring vertex improving the objective, we have found a local optimizer, which is also global.

It is useful to visualize this process, by referring to Fig. 1.3. Imagine starting from a trivially feasible solution, the origin M₀. We may improve this production plan by moving along the edges of the polyhedron; one possible path is (M₀, M₁, M₂, M₃); another path is (M₀, M₄, M₃). Both paths lead to the optimizer, even though one seems preferable, as it visits less vertices. In large-scale problems, only a small subset of vertices is actually visited to find an optimizer.

Now, in order to obtain a working algorithm, we should translate this geometric intuition into algebraic terms. The first step is transforming the LP problem in standard form,

where , , , and . Clearly, the problem makes sense only if the matrix A has less rows than columns, i.e., if m < n. If so, the set of equality constraints, regarded as a system of linear equations, has an infinite set of solutions, which may be considered as ways of expressing the right-hand side vector b as linear combinations of columns of A. Let us express A in terms of its column vectors a_j, j = 1, …, n:

where . There are n columns, but a subset B of m columns suffices to express b as follows:

To be precise, we should make sure that this subset of columns is a basis, i.e., that they are linearly independent; in what follows, we will cut a few corners and assume that this is the case. A solution of this system, in which n − m variables are set to zero, and only m are allowed to assume a nonzero value is a basic solution. Hence, we can partition the vector x into two subvectors: the subvector  of the basic variables and the subvector  of the nonbasic variables. Using a suitable permutation of the variable indices, we may rewrite the system of linear equations

as

where  is nonsingular and . However, in the LP model we have a nonnegativity restriction on variables. A basic solution where x_B ≥ 0 is called a basic feasible solution. The simplex algorithm relies on a fundamental result, that we state loosely and without proof: Basic feasible solutions correspond to extreme points of the polyhedral feasible set of the LP problem.

Solving an LP amounts to finding a way to express b as a least-cost linear combination of at most m columns of A, with nonnegative coefficients. Assume that we have a basic feasible solution x; we will consider later how to obtain an initial basic feasible solution. If x is basic feasible, it may be written as

where

The objective function value corresponding to x is

Now we should look for neighboring vertices improving this value. Neighboring vertices may be obtained by swapping a column in the basis with a column outside the basis. This means that one nonbasic variable is brought into the basis, and one basic variable leaves the basis.

To assess the potential benefit of introducing a nonbasic variable into the basis, we should express the objective function in terms of nonbasic variables. To this aim, we rewrite the objective function in (12.82), making its dependence on nonbasic variables explicit. Using Eq. (12.81), we may express the basic variables as

Then we rewrite the objective function in terms of nonbasic variables only

where

The quantities  are called reduced costs, as they measure the marginal variation of the objective function with respect to the nonbasic variables. If , it is not possible to improve the objective function; in such a case, bringing any nonbasic variable into the basis at some positive value cannot reduce the overall cost. Therefore, the current basis is optimal if . If, on the contrary, there exists a q ∈ N such that , it is possible to improve the objective function by bringing x_q into the basis. A simple strategy is to choose q such that

This selection does not necessarily result in the best performance of the algorithm; we should consider not only the rate of change in the objective function, but also the value attained by the new basic variable. Furthermore, it may happen that the entering variable is stuck to zero and does not change the value of the objective. In such a case, there is danger of cycling on a set of bases; ways to overcome this difficulty are well explained in the literature.

When x_q is brought into the basis, a basic variable must “leave” the basis in order to maintain Ax = b. To spot the leaving variable, we can reason as follows. Given the current basis, we can use it to express both b and the column a_q corresponding to the entering variable:

where B(i) is the index of the ith basic variable (i = 1, …, m), a_B(i) is the corresponding column, and

If we multiply Eq. (12.87) by a number θ and subtract it from Eq. (12.86), we obtain

From Eq. (12.89) we see that θ is the value of the entering variable in the new solution and the values of the current basic variables are affected in a way depending on the sign of d_i. If d_i ≤ 0, x_B(i) remains nonnegative when x_q increases. But if there is an index i such that d_i > 0, then we cannot increase x_q at will, since there is a limit value for which a currently basic variable becomes zero. This limit value is attained by the entering variable x_q, and the first current basic variable that gets zero leaves the basis:

If d ≤ 0, there is no limit on the increase of x_q, and the optimal solution is unbounded.

In order to start the iterations, an initial basis is needed. One possibility is to introduce a set of auxiliary artificial variables z in the constraints:

The artificial variables can be regarded as residuals, in the same vein as residuals in linear regression. Assume also that the equations have been rearranged in such a way that b ≥ 0. Clearly, a basic feasible solution of the system (12.91) where z = 0 is also a basic feasible solution for the original system Ax = b. In order to find such a solution, we can introduce and minimize an auxiliary function  as follows

using the simplex method itself. Finding an initial basic feasible solution for this artificial problem is trivial: z = b. If the optimal value of (12.92) is * = 0, we have found a starting point for the original problem; otherwise, the original problem is infeasible.

The reader is urged to keep in mind that what we have stated here is only the principle of the simplex method. Many things may go wrong with such a naive idea, and considerable work is needed to make it robust and efficient. Indeed, even though the simplex method dates back to 1947, it is still being improved today. We leave these refinement of the simplex method to the specialized literature and illustrate its application to the familiar optimal mix problem.³²

Example 12.23 A first step, which is actually carried out by good solvers, is to preprocess the formulation eliminating redundant constraints. This results in the streamlined problem

We should rewrite it in standard form, by introducing three slack variables s₁,s₂,s₃ ≥ 0:³³

In matrix terms, we have:

Finding an initial feasible solution is easy; as a starting basis, we consider {s₁,s₂,s₃}, which corresponds to a production plan where x₁ = x₂ = 0, i.e., the origin M₀ in Fig. 1.3. The corresponding basis matrix is A_B = I, i.e., a 3×3 identity matrix. Using the notation we introduced, we have

Since the cost coefficients of the basic variables are zero, application of Eq. (12.84) yields the reduced costs [−45, −60]. To improve the mix, according to Eq. (12.85) we should bring x₂ into the basis. Geometrically, we move vertically along an edge of the polyhedron, from M₀ to M₁. By the way, this is sensible as item 2 has the largest contribution to profit, but we observe that this is not necessarily the best choice, as we start moving along the longer of the two paths that lead to the optimal solution M₃; in fact, selecting the nonbasic variable with the most negative reduced cost is not the strategy of choice of modern solvers. Now we should apply Eq. (12.90) to find the variable leaving the basis. Since the matrix of basic columns is the identity matrix, applying (12.88) is easy:

and

We see that the slack variable s₃ leaves the basis, as x₂ reaches its upper bound.

Now, let us avoid terribly boring calculations and cut a long story short. Repeating these steps, we would bring x₁ into the basis, making one capacity constraint binding by driving to zero the corresponding slack variable. We know that both capacity constraints are binding at the optimal solution M₃. Let us just check that this is indeed the optimal solution. The extreme point M₃ corresponds to the basis {x₁, x₂, x₃}; therefore, we have

Finding the corresponding basic feasible solution requires solving a system of linear equations. Formally:

Now we need the reduced costs of s₁ and s₂

As both reduced costs are positive, we have indeed found the optimal basis. We note that in practice no matrix inversion is required; we may just solve systems of linear equations using by the Gaussian elimination process of Section 3.2.2.

In closing this example, the very careful reader might notice that the reduced costs above are exactly the shadow prices of capacity constraints that we have seen in Example 12.21. Indeed, there is a connection, which we cannot pursue here; we just notice that the simplex method yields shadow prices as a very useful by-product.

Fig. 12.20 Search tree for a pure binary programming problem.