Simplex Algorithm for Solving Linear Programming Problems

Solving Linear Programming Problems with the Simplex Algorithm

In 1947, George Dantzig introduced a method for solving linear programming problems, known today as the simplex method. This 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 to a linear programming problem is always a basic feasible solution. The idea of the simplex method is to find different basic feasible solutions under different basis vectors and then identify the optimal solution. Geometrically, this involves transitioning from one vertex to another until the optimal vertex is found.

Thus, the algorithm can be divided into three subproblems: 1. How to transition from one basic feasible solution to another. 2. How to determine which vertex to transition to. 3. When to stop transitioning, i.e., how to judge if the current basic feasible solution is optimal.

Transitioning Between Vertices

The standard form of a linear programming problem is:

$$\begin{align*} minimize\quad c^{T}x\\ subject\ to\quad Ax=b\\ x\ge0 \end{align*}$$

Considering the equation $Ax=b$, it can be expanded into the following canonical form:

$$\begin{align*} x_{1}+y_{1\ m+1}x_{m+1}+...+y_{1n}x_{n}=y_{10}\\ x_{2}+y_{2\ m+1}x_{m+1}+...+y_{2n}x_{n}=y_{20}\\ \vdots\\ x_{m}+y_{m\ m+1}x_{m+1}+...+y_{mn}x_{n}=y_{m0} \end{align*}$$

This system of equations can be transformed into $[I_{m},Y_{m,n-m}]x=y_{0}$. This form of the equation $Ax=b$ is called the canonical form. The solutions of the canonical form and the original system are the same. In the canonical expression, the variables corresponding to the basis vectors are basic variables, and the others are non-basic variables. In the equation $[I_{m},Y_{m,n-m}]x=y_{0}$, $x_{1},x_{2},...,x_{m}$ are basic variables, while the others are non-basic variables. Considering the augmented matrix in canonical form $[I_{m},Y_{m,n-m},y_{0}]$, the elements in the last column are the coordinates of vector $b$ with respect to the basis {$a_{1},...,a_{m}$}.

Now consider updating the augmented matrix, i.e., replacing a basic variable with a non-basic variable to find a new canonical expression for the basic variables. For example, replace the basic variable $a_{p},1\le p\le m$ with the non-basic variable $a_{q},m<q\le n$. In the original matrix, $a_{q}$ can be expressed as:

$$a_{q}=\sum_{i=1}^{m}y_{iq}a_{i}=\sum_{i=1\\i\ne p}^{m}y_{iq}a_{i}+y_{pq}a_{p}$$

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

$$a_{p}=\frac{1}{y_{pq}}a_{q}-\sum_{i=1\\i\ne p}^{m}\frac{y_{iq}}{y_{pq}}a_{i}$$

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

$$a_{j}=y_{1j}a_{1}+y_{2j}a_{2}+...+y_{mj}a_{m}$$

Substituting the expression for $a_{p}$, we get:

$$a_{j}=\sum_{i=1\\i\ne p}^{m}(y_{ij}-\frac {y_{pj}}{y_{pq}}y_{iq})a_{i}+\frac{y_{pj}}{y_{pq}}a_{q}$$

Using $y_{ij}^{'}$ to denote the elements in the new augmented matrix, we have:

$$\begin{align*} y_{ij}^{'}=y_{ij}-\frac{y_{pj}}{y_{pq}}y_{iq},i\ne p\\ y_{pj}^{'}=\frac{y_{pj}}{y_{pq}} \end{align*}$$

The transformation of the matrix using the above formulas is called a pivot transformation with respect to the element $(p,q)$. This is the formula for transitioning from one vertex to another.

Determining the Transition Vertex

Determining the transition vertex involves converting the set of vectors to another coordinate system. Specifically, it means removing one of the m vectors from the set and adding one of the n-m vectors to form a new basis. We denote this as removing vector $a_{p}$ (leaving the basis) and adding vector $a_{q}$ (entering the basis). This problem can be divided into two parts: determining $p$ and $q$.

To determine $q$, consider the pivot transformation formula:

$$\begin{align*} y_{ij}^{'}=y_{ij}-\frac{y_{pj}}{y_{pq}}y_{iq},i\ne p\\ y_{pj}^{'}=\frac{y_{pj}}{y_{pq}} \end{align*}$$

Another method to turn a non-basic variable $a_{q}$ into a basis vector is to express vector $a_{q}$ using the current basis vectors:

$$a_{q}=y_{1q}a_{1}+y_{2q}a_{2}+...+y_{mq}a_{m}$$

Multiplying both sides of the above equation by $\epsilon>0$ gives:

$$\epsilon a_{q}=\epsilon y_{1q}a_{1}+\epsilon y_{2q}a_{2}+...+\epsilon y_{mq}a_{m}$$

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

$$y_{10}a_{1}+...+y_{m0}a_{m}=b$$

Combining the above two equations, we get:

$$(y_{10}-\epsilon y_{1q})a_{1}+(y_{20}-\epsilon y_{2q})a_{2}+...+(y_{m0}-\epsilon y_{mq})a_{m}+\epsilon a_{q}=b$$

The vector $[y_{10}-\epsilon y_{1q},...,y_{m0}-\epsilon y_{mq},0,...,\epsilon,...,0]$ is a solution to the equation $Ax=b$. As $\epsilon$ increases, the first m elements of the vector will eventually have a zero, resulting in a basic feasible solution. The objective function value corresponding to the transformed basic feasible solution is:

$$\begin{align*} z=c_{1}(y_{10}-y_{1q}\epsilon)+...+c_{m}(y_{m0}-y_{mq}\epsilon)+c_{q}\epsilon\\ =z_{0}+[c_{q}-(c_{1}y_{1q}+...+c_{m}y_{mq})]\epsilon \end{align*}$$

Where $z_{0}=c_{1}y_{10}+...+c_{m}y_{m0}$. Let $z_{q}=c_{1}y_{1q}+...+c_{m}y_{mq}$, then $z=z_{0}+(c_{q}-z_{q})\epsilon$. If $z-z_{0}=(c_{q}-z_{q})\epsilon<0$, it indicates that the objective function value corresponding to the basic feasible solution has decreased, meaning that adding this vector to the basis is beneficial.

Therefore, for each $q,m<q\le n$, calculate $c_{q}-z_{q}$. Let $r_{i}=0(i=1,...,m),r_{i}=c_{i}-z_{i}(i=m+1,...,n)$, where $r_{i}$ is the reduced cost coefficient or test number. It can be proven that a basic feasible solution is optimal if and only if all the corresponding test numbers are non-negative.

Thus, by calculating the test numbers, we can determine if the current solution is optimal. If there are negative test numbers, the current solution is not optimal, and the objective function value can be reduced by adding the $q$-th vector to the basis. If there are multiple negative numbers, choose the vector with the smallest test number to add to the basis.

Having determined the method for choosing $q$, let’s discuss how to determine the exiting basis vector $a_{p}$. As mentioned earlier, for the vector $[y_{10}-\epsilon y_{1q},...,y_{m0}-\epsilon y_{mq},0,...,\epsilon,...,0]$, as $\epsilon$ increases, the first m elements will have an element $y_{i0}-\epsilon y_{iq}$ that first equals zero. We choose $a_{i}$ as the exiting basis vector $a_{p}$, i.e., $p=arg\ min_{i}\{y_{i0}/y_{iq}:y_{iq}>0\}$. Note that if there is no $y_{iq}>0$, the problem has an unbounded solution, and the computation stops. If multiple zero elements appear simultaneously, a degenerate basic feasible solution is obtained, and the smallest $p$ value is chosen.

Determining the Stopping Condition

The method for choosing the transition vertex was introduced earlier. When determining $q$, calculate the test numbers. If all test numbers are non-negative, the current basic feasible solution is optimal, and the computation stops. Additionally, when determining $p$, if as $\epsilon$ increases, the first m elements of the vector also increase, i.e., for $0<i\le m,y_{iq}<0$, the problem has an unbounded solution, and the computation stops.

Simplex Algorithm

The methods for handling the subproblems of the simplex algorithm were introduced earlier. Here is a summary of the simplex algorithm:

1. Construct the augmented matrix in canonical form based on the initial basic feasible solution. 2. Calculate the test numbers for the non-basic variables. 3. If $r_{j}\le 0$ for all $j$, stop the computation; the current basic feasible solution is optimal. Otherwise, proceed to the next step. 4. Choose the $q$ with the smallest test number among the negative test numbers. 5. If there is no $y_{iq}>0$, stop the computation; the problem has an unbounded solution. Otherwise, calculate $p=arg\ min_{i}\{y_{i0}/y_{iq}:y_{iq>0}\}$. If multiple indices $i$ satisfy the condition, let $p$ be the smallest index. 6. Perform a pivot transformation using element $(p,q)$ as the pivot element, updating the augmented matrix in canonical form. 7. Return to step 2.

Additionally, for the simplex algorithm, there is a matrix expression. To address the difficulties in finding the initial basic feasible solution, the two-phase simplex method was proposed. To reduce computation, the revised simplex method was also introduced.

To be continued…