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 with progressively greater objective
value
, 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
of the standard LP problem. This is an under-determined linear system
of
variables but only
equations, only
of its
columns of the coefficient matrix
are independent (assuming
rank
). Initially we choose the
columns of the
identity coefficient matrix
as the independent columns, but
subsequently we can choose any other
columns of
as the
independent columns and convert them into a standard basis vector by
Gauss-Jordan elimination,
together with the corresponding variables in
. The resulting
equation
remains equivalent to the original one,
i.e., the
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 columns in
,
together with the corresponding variables in
, so that the
constraining equation would always takes the following form:
(229) |
(230) |
The basic solutions corresponding the ways to choose
of the
columns of
as independent are actually the
intersections formed by any
of the
constraining hyper-planes
in the N-D space. As noted before, only a subset of these
basic solutions satisfy the
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
without exhausting all
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
at one of the vertex of the polytopic feasible region.
Specially, the implementation of the simplex algorithm is based on a
tableau with columns and
rows, initialized as below:
(231) |
At this initial stage of the iteration, the original variables
are the non-basic variables, while the
slack variables
are the basic variables
and their coefficients form an
identity matrix
,
composed of the standard basis vectors. The corresponding feasible basic
solution is simply
(232) |
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
in the last row will be
maximally increased, while all constraints remain satisfied. Here are the
steps in each iteration:
This section is based on the maximization of
. We select
in the jth column of the tableau if
is most negative, i.e.,
is most heavily weighted by
, so that
it will increase
more than any other
.
This selection is based on the constraints imposed on the selected
entering variable . In general, the restriction on
set by
the kth constraint
is
when
. If
, then
the ith constraint is most restrictive on
, we will therefore
select the corresponding basic variable
as the leaving variable,
i.e., it becomes a non-basic variable to be set to zero, so that
can be maximally increased without violating the constraints. If all
for all
, variable
is not bounded.
To convert the entering variable to a basic variable to replace the
leaving variable
, we need to turn the corresponding jth column into a
standard basis vector
. This is realized by pivoting on
in the following steps:
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
and
, the
constraint equations
always holds.
This iterative process keeps replacing the slack variables
(the basic variables initially) by the original
variables
(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 standard basis vectors
(
)
are the final basic variables that take the values in
in the
right-most column of the tableau. They form the optimal basic solution,
with the maximum of the objective function
given by the last
element also in the last column. The remaining
variables corresponding
to non-standard columns are non-basic variables that take the value zero.
When
(more constraints than variables), some of the slack variables
may remain in the basic variable group taking non-zero values; when
,
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:
In this initial state, the basic variables are ,
,
and
, the non-basic variables are
. The
corresponding basic feasible solution is at the origin.
The entering variable becomes a basic variable, and
the leaving variable
becomes a non-basic variable,
replaced by new basic variable
. The corresponding basic
feasible solution is at
, with
,
,
. The objective function value is
.
The entering variable becomes a basic variable, and
the leaving variable
becomes a non-basic variable,
replaced by the new basic variable
. The corresponding
basic feasible solution is at
, with
. The objective function value is
. Now that both
and
have
beome basic variables, the optimal solution has been obtained.
This problem has
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
at the origin with
, through
the vertex at
and
with
, to the optimal solution
at
and
with
.
The Matlab code for implementing the simplex algorithm is
listed below. The function takes input including matrix
as well as the vectors
and
, and generates
the optimal solution
, the corresponding maximum value
of the ojbective function
, together
with the values of the slack vaerible
.
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