The Simplex Algorithm for Solving Linear Programming Problems

In 1947, Dantzig proposed a method for solving linear programming problems, now known as the simplex method. It is a concise and efficient algorithm, hailed as one of the ten algorithms with the greatest impact on scientific development and engineering practice in the 20th century. As mentioned earlier, the optimal solution of a linear programming problem must be a basic feasible solution. The idea of the simplex method is to compute different basic feasible solutions under different basis vectors, and then find the optimal solution. From a geometric perspective, this is the process of moving from one extreme point to another until the optimal extreme point is found. In that case, the algorithm can be divided into three subproblems: 1. How to move from one basic feasible solution to another. 2. How to determine which extreme point to move to. 3. When to stop the moving operation, i.e., how to determine whether the current basic feasible solution is optimal.

Moving Between Extreme Points

The standard form of a linear program is

Consider the equation $Ax=b$, expanded into the following canonical system of equations

This system can be converted into $[I_{m},Y_{m,n-m}]x=y_{0}$. A system of equations $Ax=b$ in this form is called the canonical form. The canonical form of a system has the same solutions as the original system. In the canonical expression, the variables corresponding to the basis column vectors are the basic variables, and the others are the nonbasic variables. That is, in the equation $[I_{m},Y_{m,n-m}]x=y_{0}$, $x_{1},x_{2}...,x_{m}$ are the basic variables and the rest are nonbasic variables. Consider the canonical augmented matrix $[I_{m},Y_{m,n-m},y_{0}]$; the elements of its last column are the coordinates of the vector $b$ with respect to the basis {$a_{1},...,a_{m}$}.

Now consider updating the augmented matrix, i.e., replacing some basic variable with some nonbasic variable and finding the canonical expression corresponding to the new basis. For example, replace the basic variable $a_{p},1\le p\le m$ with the nonbasic variable $a_{q},m<q\le n$. In the original matrix, $a_{q}$ can be expressed as

Note that the vector set {$a_{1},...,a_{p-1},a_{q},a_{p+1},...,a_{m}$} is linearly independent and can form a basis only when $y_{pq}\ne 0$. From the above equation we can solve for

Using the original augmented matrix, any vector $a_{j}(m<j\le n)$ can be expressed as

Substituting the expression for $a_{p}$ yields

Letting $y_{ij}^{'}$ denote the elements of the new augmented matrix, the above gives

The transformation of the matrix using the above formulas is called the pivot transformation about the element $(p,q)$. The above gives the formula for moving from one extreme point to another.

Determining the Extreme Point to Move To

Determining the extreme point to move to means transforming the vector set into another coordinate system. In detail, this means taking one of the $m$ vectors out of the vector set and selecting one vector from the other $n-m$ vectors to add to the basis vector set, forming a new basis. We denote this as taking out the vector $a_{p}$ (leaving the basis) and adding the vector $a_{q}$ (entering the basis). This problem can again be split into two parts, namely determining $p$ and $q$ separately. For determining $q$, first look at the pivot transformation formula below

Now consider another way to turn a nonbasic variable $a_{q}$ into a basis vector. Express the vector $a_{q}$ using the current basis vectors

Multiplying both sides by $\epsilon>0$ yields

Substituting the basic solution $x=[y_{10},...,y_{m0},0,...,0]^T$ into the equation $Ax=b$ yields

Combining the two equations above gives

That is, the vector $[y_{10}-\epsilon y_{1q},...,y_{m0}-\epsilon y_{mq},0,...,\epsilon,...,0]$ is a solution of the equation $Ax=b$. As $\epsilon$ continually increases, when a 0 first appears among the first $m$ elements of the vector, we obtain a basic feasible solution. The objective function value corresponding to the transformed basic feasible solution is

where $z_{0}=c_{1}y_{10}+...+c_{m}y_{m0}$. Letting $z_{q}=c_{1}y_{1q}+...+c_{m}y_{mq}$, we obtain $z=z_{0}+(c_{q}-z_{q})\epsilon$. When $z-z_{0}=(c_{q}-z_{q})\epsilon<0$, the objective function value corresponding to the basic feasible solution has decreased, meaning that adding this vector to the basis vector set is beneficial.

So we can compute $c_{q}-z_{q}$ for each $q,m<q\le n$. Letting $r_{i}=0(i=1,...,m),r_{i}=c_{i}-z_{i}(i=m+1,...,n)$, we call $r_{i}$ the $i$-th reduced cost coefficient, also known as the reduced cost (testing number). It can be proven that a basic feasible solution is optimal if and only if all corresponding reduced costs are nonnegative.

That is, we can determine whether the current solution is optimal by computing the reduced costs. When there are negative numbers among the reduced costs, the current solution is not optimal, and the objective function value can be decreased by adding the $q$-th vector to the basis vector set. If there are multiple negative numbers, we choose the $q$-th vector with the smallest reduced cost to add to the basis vector set.

The above determined the method for choosing $q$. Now let us introduce how to determine the leaving vector $a_{p}$. As shown above, for the vector $[y_{10}-\epsilon y_{1q},...,y_{m0}-\epsilon y_{mq},0,...,\epsilon,...,0]$, as $\epsilon$ continually increases, among the first $m$ elements the $i$-th element $y_{i0}-\epsilon y_{iq}$ will be the first to equal 0. Then we choose $a_{i}$ as the leaving vector $a_{p}$, i.e., $p=arg\ min_{i}\{y_{i0}/y_{iq}:y_{iq}>0\}$. Note that if no $y_{iq}>0$ exists, then the problem has an unbounded solution, so the computation stops. Additionally, if multiple 0 elements appear simultaneously, then a degenerate basic feasible solution is obtained; in this case we choose the smallest value of $p$.

Determining the Stopping Condition

The above introduced how to choose the extreme point to move to. When determining $q$, we compute the reduced costs; if all reduced costs are nonnegative, then the current basic feasible solution is optimal and the computation stops. Moreover, when determining $p$, if as $\epsilon$ continually increases the first $m$ elements of the vector also increase, i.e., $y_{iq}<0$ for $0<i\le m$, then the problem has an unbounded solution, and the computation also stops.

The Simplex Algorithm

The above introduced the methods for handling the several subproblems of the simplex algorithm. Now let us summarize the simplex algorithm. 1. Construct the canonical augmented matrix from the initial basic feasible solution. 2. Compute the reduced costs of the nonbasic variables. 3. If $r_{j}\le 0$ for all $j$, then stop the computation; the current basic feasible solution is optimal. Otherwise, proceed to the next step. 4. Among the reduced costs less than 0, choose the $q$ with the smallest reduced cost. 5. If no $y_{iq}>0$ exists, then stop the computation; the problem has an unbounded solution. Otherwise, compute $p=arg\ min_{i}\{y_{i0}/y_{iq}:y_{iq>0}\}$. If multiple indices $i$ satisfying the condition are obtained, let $p$ be the smallest index value. 6. Using the element $(p,q)$ as the pivot element, perform the pivot transformation and update the canonical augmented matrix. 7. Go to step 2.

In addition, the simplex algorithm also has a matrix expression. To address the difficulties encountered in finding an initial basic feasible solution, the two-phase simplex method was proposed; to reduce the amount of computation, the revised simplex method was proposed; and so on.

To be continued…