The Simplex Algorithm

The simplex algorithm finds the optimal solution of a LP problem by an iterative process that traverses along a sequence of edges of the polytopic feasible region, starting at the origin and through a sequence of vertices ${\bf x}$ with progressively greater objective value $f({\bf x})$, until eventually reaching the optimal solution. By doing so, it avoids checking exhaustively all vertices of the feasible region for optimality.

We first consider the equality constraint ${\bf Ax}=[{\bf A}_{M\times N}\;\vert\;{\bf I}_{M\times M}] {\bf x}={\bf b}$ of the standard LP problem. This is an under-determined linear system of $N+M$ variables but only $M$ equations, only $M$ of its $N+M$ columns of the coefficient matrix ${\bf A}$ are independent (assuming rank$({\bf A})=M$). Initially we choose the $M$ columns of the identity coefficient matrix ${\bf I}$ as the independent columns, but subsequently we can choose any other $M$ columns of ${\bf A}$ as the independent columns and convert them into a standard basis vector by Gauss-Jordan elimination, together with the corresponding variables in ${\bf x}$. The resulting equation ${\bf Ax}={\bf b}$ remains equivalent to the original one, i.e., the $M$ equality constraints are always preserved and never violated in the process.

For convenience of discussion and without loss of generality, we could reorder (not actually carried out) the $N+M$ columns in ${\bf A}$, together with the corresponding variables in ${\bf x}$, so that the constraining equation would always takes the following form:

$${\bf A}{\bf x}=[{\bf A}_n\;\;{\bf I}] \left[\begin{array}{c}{\bf x}_n\\ {\bf x}_b\end{array}\right] ={\bf A}_n{\bf x}_n+{\bf I}\;{\bf x}_b={\bf b}$$ (229)

Now the variables in the M-D vector ${\bf x}_b$ corresponding to the $M$ independent columns in ${\bf I}$ are the basic variables, while the N-D vector ${\bf x}_n$ corresponding to the remaining $N$ dependent columns, now denoted by ${\bf A}_n$, are the non-basic variables. This equation will always hold if ${\bf x}_n={\bf0}$ and ${\bf x}_b={\bf b}$. Such a solution is called a basic solution of the linear system ${\bf A}{\bf x}={\bf b}$:

$${\bf x}=\left[\begin{array}{c}{\bf x}_n\\ {\bf x}_b\end{array}\right] =\left[\begin{array}{c}{\bf0}\\ {\bf b}\end{array}\right]$$ (230)

Initially, ${\bf A}_n$ is the coefficient matrix in the inequality constraint ${\bf Ax}\le {\bf b}$, and the corresponding non-basic variables ${\bf x}_n={\bf0}$ is actually the origin. But through the iteration, we will keep converting some other columns of ${\bf A}$ into standard basis vectors $\{{\bf e}_1,\cdots,{\bf e}_M\}$, the column vectors of ${\bf I}$, by Gauss-Jordan elimination. The variables in ${\bf x}$ corresponding to the new standard basis vectors become the basic variables, while those corresponding to columns that become non-standard basis vectors become non-basis variables.

The basic solutions corresponding the $C_{N+M}^N$ ways to choose $M$ of the $N+M$ columns of ${\bf A}$ as independent are actually the intersections formed by any $N$ of the $N+M$ constraining hyper-planes in the N-D space. As noted before, only a subset of these $C_{N+M}^N$ basic solutions satisfy the $M$ constraints, and they are called basic feasible solutions (BFS). The goal of the iterative process of the simplex algorithm is to find the optimal basic feasible solution that maximizes $f({\bf x})$ without exhausting all $C_{N+M}^N$ possibilities. This is done by selecting the columns in such a way that the value of the objective function will always be maximally increased, until eventually we find the optimal solution ${\bf x}^*$ at one of the vertex of the polytopic feasible region.

Specially, the implementation of the simplex algorithm is based on a tableau with $N+M+1$ columns and $M+1$ rows, initialized as below:

• The first $M$ rows are for the coefficients in the equality constraints $\sum_{j=1}^N a_{ij}x_j+s_i=b_i$ ( $i=1,\cdots,M$). An additional $(M+1)st$ row at the bottom is for the coefficients in the objective function $f({\bf x})=\sum_{j=1}^N c_jx_j$, which is zero initially at the origin $x_1=\cdots=x_N=0$, i.e., $z-c_1x_1-\cdots-c_Nx_N=0$.
• The first $N+M$ columns contain the coefficients for the $N$ original variables $\{x_1,\cdots,x_N\}$ and the $M$ slack variables $\{s_1,\cdots,s_M\}$, and an additional $(N+M+1)st$ column is for the constants ${\bf b}=[b_1,\cdots,b_M]^T$ on the right-hand side of the $M$ constraint equations. The last element of the column is the right-hand side of the objective equation $z=f({\bf x})$, which is initially zero.

$$\begin{tabular}{c\vert\vert cccc\vert cccc\vert\vert c}\hline bas... ...-c_2$\ & $\cdots$\ & $-c_M$\ & 0 & 0 & $\cdots$\ & 0 & 0\\ \hline \end{tabular}$$ (231)

At this initial stage of the iteration, the $N$ original variables ${\bf x}_n=[x_1,\cdots,x_N]^T$ are the non-basic variables, while the $M$ slack variables ${\bf x}_b=[s_1,\cdots,s_M]^T$ are the basic variables and their coefficients form an $M\times M$ identity matrix ${\bf I}$, composed of the standard basis vectors. The corresponding feasible basic solution is simply

$${\bf x}_b=[s_1,\cdots,s_M]^T={\bf b},\;\;\;\;\; {\bf x}_n=[x_1,\cdots,x_N]^T={\bf0}$$ (232)

that satisfies all the constraints ${\bf Ax}={\bf A}_n{\bf x}_n+{\bf I}{\bf x}_b={\bf b}$.

In each of the subsequent iterations, we will select one of the non-basic variables, called the entering variable, to replace one of the basic variables, called the leaving variable, in such a way that the value of the objective function value $z=f({\bf x})$ in the last row will be maximally increased, while all constraints remain satisfied. Here are the steps in each iteration:

• Selection of the entering variable

  This section is based on the maximization of $z=f({\bf x})$. We select $x_j$ in the jth column of the tableau if $-c_j$ is most negative, i.e., $x_j$ is most heavily weighted by $c_j=\max\{c_1,\cdots,c_N\}$, so that it will increase $z=\sum_{j=1}^N c_jx_j$ more than any other $x_k\;\;(k\ne j)$.

• Selection of the leaving variable

  This selection is based on the constraints imposed on the selected entering variable $x_j$. In general, the restriction on $x_j$ set by the kth constraint $\sum_{j=1}^N a_{kj}x_j\le b_k$ is $x_j\le b_k/a_{kj}$ when $a_{kj}>0$. If $b_i/a_{ij}=\min\{b_k/a_{1j},\cdots,b_k/a_{Mj}\}$, then the ith constraint is most restrictive on $x_j$, we will therefore select the corresponding basic variable $s_i$ as the leaving variable, i.e., it becomes a non-basic variable to be set to zero, so that $x_j$ can be maximally increased without violating the constraints. If all $a_{kj}<0$ for all $k=1,\cdots,M$, variable $x_j$ is not bounded.

• Gauss-Jordan elimination based on Pivoting

  To convert the entering variable $x_j$ to a basic variable to replace the leaving variable $s_i$, we need to turn the corresponding jth column into a standard basis vector ${\bf e}_i$. This is realized by pivoting on $a_{ij}$ in the following steps:

  1. Divide the pivot row ${\bf r}_i$ by $a_{ij}$: ${\bf r}_i\leftarrow{\bf r}_i/a_{ij}$. Now $a_{ij}=1$.
  2. Subtract $a_{kj}{\bf r}_i$ from the kth row: ${\bf r}_k\leftarrow{\bf r}_k-a_{kj}{\bf r}_i$ so that $a_{kj}=0$ for all $k=1,\cdots,M+1$, $k\ne i$, including the last row ${\bf r}_{M+1}$.
  Now the jth column corresponding to the entering variable $x_j$ becomes a standard basis vector ${\bf e}_i$, i.e., $x_j$ becomes a basic variable that takes the value of $b_i$, and the column corresponding to the leaving variable $s_i$ is no longer a standard basis vector, and $s_i=0$ as it is now a non-basic variable.

As Gauss-Jordan elimination converts the linear equations to a set of equivalent equations, the constraints remain satisfied through out the process. Although the membership of the basic and non-basic variable groups keeps changing in the iteration, so long as ${\bf x}_b={\bf b}$ and ${\bf x}_n={\bf0}$, the $M$ constraint equations ${\bf Ax}={\bf A}_n{\bf x}_n+{\bf I}{\bf x}_b={\bf b}$ always holds.

This iterative process keeps replacing the slack variables $\{s_1,\cdots,s_M\}$ (the basic variables initially) by the original variables $\{x_1,\cdots,x_N\}$ (the non-basic variables initially), one at a time, until all original variables become basic variables and all elements in the last row are non-negative.

The final result can be read out directly from the tableau. The variables corresponding to the $M$ standard basis vectors ${\bf e}_i$ ( $i=1,\cdots,M$) are the final basic variables that take the values in ${\bf b}$ in the right-most column of the tableau. They form the optimal basic solution, with the maximum of the objective function $z=f({\bf x})$ given by the last element also in the last column. The remaining $n$ variables corresponding to non-standard columns are non-basic variables that take the value zero. When $M>N$ (more constraints than variables), some of the slack variables may remain in the basic variable group taking non-zero values; when $M<N$, some of the original variable may be in the non-basic group taking the value zero.

Example:

We re-solve the LP problem considered in the previous examples, now in the standard form:

$$\begin{tabular}{ll} max: & $f({\bf x})=2x_1+3x_2$\ \\ s. t.: & $... ...; x_2\ge 0,\; s_1\ge 0,\; s_2\ge 0,\; s_3\ge0 \end{array}\right.$ \end{tabular}$$

The standard form is further converted to a tableau as shown below. The left-most column indicates the basic variables, the next $N=2$ columns are for the variables $x$ and $y$,the next $M=3$ columns are for the slack variables $u$, $v$, and $w$, the right most column is for the constants $b_1$, $b_2$, and $b_3$ on the right-hand side of the equations.

$$\begin{tabular}{c\vert cc\vert ccc\vert c} \hline & $x_1$\ & $x_2... ...& 0 & 0 & 1 & 40 \\ \hline $z$\ &-2 &-3 & 0 & 0 & 0 & 0 \\ \hline \end{tabular}$$

In this initial state, the basic variables are $s_1=18$, $s_2=60$, and $s_3=40$, the non-basic variables are $x_1=x_2=0$. The corresponding basic feasible solution is at the origin.

• Select column $j=2$ as the pivot column, as $-c_2=-3<-c_1$ is the most negative coefficient.
• Select ${\bf r}_3$ as the pivot row, as $b_3/a_{32}=40/5=8=\min\{18/1=18,\;60/5=12,\;40/5=8\}$.
• Divide pivot row ${\bf r}_i={\bf r}_3=[2,\;5,\;0,\;0,\;1,\;40]$ by the pivot element $a_{ij}={32}=5$ to get ${\bf r}_i={\bf r}_3=[0.4,\;1,\;0,\;0,\;0.2,\;8]$.
• Subtract $a_{kj}{\bf r}_k=a_{k2}{\bf r}_k$ from row ${\bf r}_k$, so that $a_{kj}=a_{k2}=0,\;(k=1,\cdots,M+1=4,\,k\ne i)$, including the last row ${\bf r}_4$.

$$\begin{tabular}{c\vert cc\vert ccc\vert c} \hline & $x_1$\ & $x_2... ...0 & 0 & 0.2& 8 \\ \hline $z$\ &-0.8& 0 & 0 & 0 &0.6& 24 \\ \hline \end{tabular}$$

The entering variable $x_j=x_2$ becomes a basic variable, and the leaving variable $s_i=s_3$ becomes a non-basic variable, replaced by new basic variable $x_2$. The corresponding basic feasible solution is at $x_1=0,\;x_2=b_3=8$, with $s_1=10$, $s_2=20$, $s_3=0$. The objective function value is $z={\bf c}^T{\bf x}=24$.

• Select column $j=1$ as the pivot column, as $-c_1=-0.8$ is the most (only) negative coefficient in the last row.
• Select row $i=2$ as the pivot row, as the ratio $b_2/a_{21}=20/4=5 =\min\{10/1.6=6.25,\;20/4=5,\;8/0.4=20\}$.
• Divide pivot row ${\bf r}_i={\bf r}_2=[4,\;0,\;0,\;1,\;-1,\;20]$ by the pivot element $a_{ij}=a_{21}=4$ to get ${\bf r}_i={\bf r}_2=[1,\;0,\;0,\;0.25,\;-0.25,\;5]$.
• Subtract $a_{kj}{\bf r}_k=a_{k1}{\bf r}_k$ from row ${\bf r}_k$, so that $a_{kj}=a_{k1}=0, (k=1,\cdots,M+1=4)$, including the last row ${\bf r}_4$.

$$\begin{tabular}{c\vert cc\vert ccc\vert c} \hline & $x_1$\ & $x_2... ...1 & 0.3 & 6 \\ \hline $z$\ & 0 & 0 & 0 & 0.2 & 0.4 & 28 \\ \hline \end{tabular}$$

The entering variable $x_1$ becomes a basic variable, and the leaving variable $s_1$ becomes a non-basic variable, replaced by the new basic variable $x_1$. The corresponding basic feasible solution is at $x_1=b_2=5,\;x_2=b_3=6$, with $s_1=2,\;s_2=s_3=0$. The objective function value is $z={\bf c}^T{\bf x}=28$. Now that both $x_1$ and $x_2$ have beome basic variables, the optimal solution has been obtained.

This problem has $C_{N+M}^N=C_5^2=10$ basic solutions, corresponding to the same number of intersections, out of which five are feasible, as previously obtained. The simplex method finds three of these feasible solutions, starting from $x_1=x_2=0$ at the origin with $z=0$, through the vertex at $x_1=0$ and $x_2=8$ with $z=24$, to the optimal solution at $x_1=5$ and $x_2=6$ with $z=28$.

The Matlab code for implementing the simplex algorithm is listed below. The function takes input including matrix ${\bf A}$ as well as the vectors ${\bf b}$ and ${\bf c}$, and generates the optimal solution ${\bf x}^*$, the corresponding maximum value of the ojbective function $f({\bf x})={\bf c}^T{\bf x}^*$, together with the values of the slack vaerible ${\bf s}$.

function [x s z] = MySimplex(A,b,c)
    [m n]=size(A)                       % coefficient matrix for variables x
    x=zeros(n,1);                       % initial zero values for x
    s=zeros(m,1);
    A=[A eye(m) b; -c zeros(1,m+1)];    % initialization of tableau 
          
    [cmin, pc]=min( A(m+1,1:n+m) );     % pc, index of the first pivot column 
    it=0;
    while cmin<0                        % iteration until all c's in last row are zero
        it=it+1;
        rmin=9e9;
        for i=1:m                       % find pivot row
            if A(i,n+m+1)~=0 & A(i,pc)>0
                w=A(i,n+m+1)./A(i,pc);
                if w<rmin
                    pr=i;               % pr, index of the pivot row
                    rmin=w;
                end
            end
        end
        p=A(pr,pc);                     % get the pivot element
        A(pr,:)=A(pr,:)/p;              % modify pivot row
        for i=1:m+1                     % modify all k+1 rows
            if i~=pr
                A(i,:)=A(i,:)-A(pr,:)*A(i,pc);
            end
        end
        z=A(m+1,n+m+1);                 % maximum value of objective function so far
        
        [cmin, pc]=min( A(m+1,1:n+m) ); % pivot column of next iteration
    end
end