## Monday, August 12, 2024

### 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 firstoption 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 minminvs(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    rhsc     -3   -2    1   -1  r1     1    0.5            1            10r2     2    1   -1   -1        1        12 r3                    1            1     2; * extract data from problemParameters  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 basisbasis(i) = yes;nonbasic(k) = not basis(k); * initial:* structural variables  = 0* slacks are equal to rhsx(i) = rhs(i); * initial : unit or all-slack basisB(i,i) = 1; * initially the unit matrixinvB(i,i) = 1; * initialbasmap(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
```