Revised Simplex LP Solver written in GAMS

I am teaching some GAMS classes, and a question arose: "How does the Simplex method work?" It's not easy to answer in a few sentences, but I want to touch upon the concept of a basis anyway. Once you have a good intuition of what a basis is, a simple Simplex method is not so far-fetched. I find the tableau presentation somewhat confusing and far removed from what actual Simplex solvers do. I strongly prefer the Revised Simplex Method in matrix notation. 

Minor rant: I just don't understand the appeal of the tableau method. It looks to me like an invention for torturing undergrad students. Most of all, it is not very structure-revealing; it does not help you understand the underlying concepts. But about 100% of the LP textbooks insist we should learn that first.

As a gimmick, I implemented a simplified version in the GAMS language. This reminds me that someone spent the effort writing a Basic interpreter in TeX [1]. This is probably just as useful.

References

1. Andrew Marc Greene, BASIX -An Interpreter Written in TeX.  TUGboat, Volume 11 (1990), No. 3-Proceedings of the 1990 Annual Meeting. https://tug.org/tugboat/tb11-3/tb29greene.pdf.

Appendix

Revised Simplex in GAMS

$onText       Revised Simplex Method in GAMS         basic version       Assumptions:        minimization        positive variables x>=0        <= constraints        x=0 is feasible             $offText   *----------------------------------------------------------- * scalar LP *-----------------------------------------------------------   positive variables x1,x2,x3,x4; variable z; equations   obj   e1   e2   e3 ;   obj..  z =e= 3*x1+2*x2-x3+x4; e1..   2*x1+x2 =l= 20; e2..   x1+x2-x3-x4 =l= 12; e3..   x4 =l= 2;   model m /all/; solve m maximizing z using lp;   *----------------------------------------------------------- * example of Singleton Set * this concept is used later on *----------------------------------------------------------- * in case of multiple assignment, pick first option strictSingleton=0; set s /a,b,c,d/; alias (s,ss); parameter vs(s) /a 1, b 2, c 3, d 1/; singleton set minvs(s); * implementation of arg min minvs(s) = smin(ss,vs(ss)) = vs(s); display "singleton set",minvs;     *----------------------------------------------------------- * simplex method (simplified) * * problem setup *-----------------------------------------------------------   Set   k 'columns (including slacks)' /x1*x4,r1*r3/   i(k) 'rows' /r1*r3/   j(k) 'columns' /x1*x4/ ;   * specify problem in matrix format * add slacks (forming EQ constraints) * convert to minimization * we assume all variables (including slacks are >= 0) Table problem(*,*) 'includes cost (c) and rhs (b)' *     ---structural----   --slacks--           x1   x2   x3   x4   r1  r2  r3    rhs c     -3   -2    1   -1   r1     1    0.5            1            10 r2     2    1   -1   -1        1        12  r3                    1            1     2 ;   * extract data from problem Parameters   c(k) 'cost coefficients'   rhs(i) 'right-hand side'   A(i,k) 'LP matrix' ; c(k) = problem('c',k); rhs(i) = problem(i,'rhs'); A(i,k) = problem(i,k); display c,rhs,A;     *----------------------------------------------------------- * initial basis *-----------------------------------------------------------   alias (i,ii);   set    dummy     'for ordering of display output' /init/    iter      'max number of iterations' /iter1*iter10/    basis(k)  'number of elements = number of rows'    nonbasic(k) 'non-basic variables'    basmap(i,k) 'basis mapping i-th basis var = k-th problem var' ;   parameter    x(k)     'all levels'    B(i,ii)  'Basis matrix'    invB(i,ii) 'basis inverse'    xb(i)    'values of basic variables'    cb(i)    'mapped cost' ;     * feasible starting basis basis(i) = yes; nonbasic(k) = not basis(k);   * initial: * structural variables  = 0 * slacks are equal to rhs x(i) = rhs(i);   * initial : unit or all-slack basis B(i,i) = 1;   * initially the unit matrix invB(i,i) = 1;   * initial basmap(i,i) = yes;   xb(i) = sum(basmap(i,k),x(k)); cb(i) = sum(basmap(i,k),c(k))   parameter trace(*,*); trace('init',k)=x(k); trace('init','obj')=sum(i,cb(i)*xb(i));   *----------------------------------------------------------- * iterating *-----------------------------------------------------------     parameters    itno     'iteration number'    sigma(k) 'reduced cost'    minsigma 'min reduced cost of nonbasics'    opttol   'optimality tolerance' /1e-7/    transf(i) 'transformed column invB*A(entering)'    step(i)   'steps in ratio test'    minstep   'most limiting step'    objval    'objective value' ;   singleton Sets    enter(k) 'entering variable (non-basic -> basic)'    leave(i) 'leaving variable (basic -> non-basic)' ;   loop(iter,      * calculate the reduced cost * for basic variables: 0 * for nonbasics j: c(j) - cb*invB*A_j    sigma(basis) = 0;    sigma(nonbasic) = c(nonbasic)-sum(i,sum(ii,cb(ii)*invB(ii,i))*A(i,nonbasic));      * find most negative reduced cost    minsigma = smin(nonbasic,sigma(nonbasic));    if(minsigma >= -opttol,       display trace,"optimal";       break;    );    * entering Variable * note that enter() is a singleton set, * so we find one single entering variable    enter(nonbasic) = sigma(nonbasic) = minsigma;    * leaving variable    transf(i) = sum(ii, invB(i,ii)*A(ii,enter));    if(smax(i,transf(i))<=0,       display trace,"unbounded";       break;    );    step(i) = (xb(i)/transf(i))$(transf(i)>0);    minstep = smin(i$(transf(i)>0),step(i));    leave(i) = step(i) = minstep;      * basis change bookkeeping    basis(enter) = yes;    nonbasic(enter) = no;    loop(basmap(leave,k),       basis(k) = no;       nonbasic(k) = yes;    );    basmap(leave,k) = no;    basmap(leave,enter) = yes;      loop(basmap(ii,k),       B(i,ii) = A(i,k);    );    * invert from scratch (real simplex solvers are smarter)    invB(i,ii)=0;    executeTool.checkErrorLevel 'linalg.invert i B invB';       cb(i) = sum(basmap(i,k),c(k));    xb(i) = sum(ii,invB(i,ii)*rhs(ii));    x(k) = sum(basmap(i,k),xb(i));      objval = sum(i, cb(i)*xb(i));      trace(iter,k)=x(k);    trace(iter,'obj')=sum(i,cb(i)*xb(i));   );

 

The trace log illustrates nicely that we always have m=3 basics:

 

----    171 PARAMETER trace  

               x1          x2          x3          x4          r1          r2          r3         obj

init                                                       10.000      12.000       2.000
iter1       6.000                                           4.000                   2.000     -18.000
iter2       7.000                               2.000       3.000                             -23.000
iter3                  14.000                   2.000       3.000                             -30.000
iter4                  20.000       6.000       2.000                                         -36.000