This computes 

min  c x
       A x = b
       lb <= x <= ub
       
through the Simplex method.

How to use this notebook

First you need input the problem data in the section Problem data. You are required to provide the matrix A, the vectors lb (lower bounds), ub (upper bounds), c (costs), b (rhs). Other quantities will be computed automatically, like m (# rows) and n (# columns).

Then you are ready to solve Phase One. To this aim you have to set up the data for Phase One in the section Setting Phase One.
In this section you needn't provide data. Everything is computed automatically.

Now you may iterate the simplex steps to find an initial feasible solution. You do this in the section Iteration. (see instructions within the section).

Once a feasible solution is found the data must be reset to solve Phase Two.
This is done in Resetting for Phase Two. Follow carefully instructions within the section.

At this point you go back to Iteration and carry out the simplex steps until you find a solution (the problem cannot be unbounded because of the bounds on the variables).

Problem data

A manual instance

A = {{ 1,-3,2, 1,-1,1,2, 1, 0,1,0,0, 0},
      { 2, 1,5, 2, 1,0,0, 0, 1,0,1,0, 0},
      { 1, 0,2,-1, 0,1,0, 2, 3,0,0,1, 0},
      {-2,-1,3, 0, 0,0,0,-1,-2,0,0,0,-1}};

m=Length[A];n=Length[A[[1]]];

If a bound is infinity do not use Infinity. Use instead a large number. It will be treated as a normal bound.

inf =100000;

lb ={0,0,0,0,0,0,0,0,0,0,0,0,0};
ub={5,5,4,3,3,5,6,2,4,inf,inf,inf,inf};

c = {1, -3 , -1, 4,-3, 10,-6,-4,-1,0,0,0,0};

b = {8,5,4,-6};

If[Length[b]!=m, Print["Check b !"]];
If[Length[c]!=n, Print["Check c !"]];
If[Length[lb]!=n, Print["Check lb !"]];
If[Length[ub]!=n, Print["Check ub !"]];

A relaxed TSP instance

The arc costs are uniformly randomly  generated
between 1 and 99.

inf=100000;

nn=5;na=nn(nn-1)/2;

SeedRandom[10];

bb=Table[Random[Integer,{1,99}],{na}]

costs=Table[0,{nn+1},{nn+1}];
k=1;
Do[costs[[1,i+1]]=i;costs[[i+1,1]]=i,{i,nn}];
Do[costs[[i+1,j+1]]=bb[[k]];
   costs[[j+1,i+1]]=bb[[k++]],
{i,nn-1},{j,i+1,nn}];
MatrixForm[costs]

n=na;m=nn;
A=Table[0,{i,nn},{j,na}];
k=1;
Do[A[[i,k]]=1;A[[j,k]]=1;
   k++;
   ,{i,nn-1},{j,i+1,nn}];
AT=Transpose[A];
c=bb;
b=Table[2,{nn}];
lb=Table[0,{na}];
ub=Table[1,{na}];

MatrixForm[A]

Setting Phase One

Off[General::spell1];

AAA=A;
cc=c;
nn=n;

x0=lb

res=b-A.x0

sign=Table[If[res[[i]]<0,-1,1],{i,m}]

AT=Join[Transpose[A],DiagonalMatrix[sign]];

A=Transpose[AT];

c=Table[If[i<=n,0,1],{i,n+m}]

lb=Join[lb,Table[0,{m}]];

ub=Join[ub,Table[inf,{m}]];

n=n+m;

base = Table[i,{i,nn+1,n}]

eta1=Table[i,{i,nn}]

eta2={};

resume={};

Resetting for Phase Two

first compute this

base

value

artinbase=Select[base,#>nn &]

If different from {} and value=0 then there are degenerate artificial vars in base.
Do the following for each one of them

artinbase

artificial index to be pushed out

indexout=artinbase[[1]];

artinbase=Drop[artinbase,1];

baseout=Position[base,indexout][[1,1]]

index to be brought in (new matrix must be non singular and only indices with non zero values in the following vector must be considered)

WAi=(W.A)[[baseout]]

notzero=Select[WAi, #!=0 &]

indexin=Position[WAi,notzero[[1]]][[1,1]]

base updating

base[[baseout]]=indexin;

AppendTo[eta1,indexout]

eta1=Complement[eta1,{indexin}]

If there are more artificial vars in base redo the computation

Now do the following and go to Iteration

eta1=Sort[eta1]

eta1=Drop[eta1,-m]

A=AAA;

c=cc;
n=nn;

lb=Drop[lb,-m];ub=Drop[ub,-m];AT=Transpose[A];

Iteration

Current base

{base,eta1,eta2}

Computation of reduced costs

B is the basis matrix

B = Transpose[Transpose[A][[ base]]];
MatrixForm[B];

W = Inverse[B];
MatrixForm[W];

cb = c[[ base]]

Which[eta1=={}
     ,x = W. (b-Transpose[AT[[eta2]]].ub[[eta2]]);
      value=cb.x+c[[eta2]].ub[[eta2]]
     ,eta2=={}
     ,x = W. (b-Transpose[AT[[eta1]]].lb[[eta1]]);
      value=cb.x+c[[eta1]].lb[[eta1]]
     ,True
     ,x = W. (b-Transpose[AT[[eta2]]].ub[[eta2]]
               -Transpose[AT[[eta1]]].lb[[eta1]]);
      value=cb.x+
         c[[eta2]].ub[[eta2]]+c[[eta1]].lb[[eta1]]
      ];

current x in base and current value

{x,N[value]}

AppendTo[resume,{{base,eta1,eta2},value}];

dual variables

u = cb . W

reduced costs

cr = c - u . A

Selection of entering column

enteringset1=Select[
             Transpose[{N[cr],Range[n]}],
             (#[[1]] < 0)&&  (MemberQ[eta1,#[[2]]])&];

enteringset2=Select[
             Transpose[{N[cr],Range[n]}],
             (#[[1]] > 0)&&  (MemberQ[eta2,#[[2]]])&];

If[enteringset1=={} && enteringset2=={},
   Print["optimal base !"]]

variables at lower bounds eligible to enter the basis:
first row    -> reduced costs
second row -> variable indices

If[enteringset1 != {},TableForm[Transpose[enteringset1]]]

variables at upper bounds eligible to enter the basis:

If[enteringset2 != {},TableForm[Transpose[enteringset2]]]

Computation of leaving index and base updating

assign to p the value of the index to enter the base. The rest is computed automatically

p =6

If[MemberQ[eta2,p]
  ,etafrom=2;
   ppos=Flatten[Position[eta2,p]][[1]]
  ,etafrom=1;
   ppos=Flatten[Position[eta1,p]][[1]]]

Block[{yy},
    y = W . Transpose[A][[ p ]];
    If[etafrom==2,yy=-y;y=yy;];
    y
]

x

Selection of variable to leave the basis

Block[{i,ratio,min},
      min = 1000000;
      Do[Which[y[[i]]>0
               ,ratio = (x[[i]]-lb[[base[[i]]]])/y[[i]];
                If[ratio < min,
                   min = ratio; q = i;etato=1]
               ,y[[i]]<0
               ,ratio = (x[[i]]-ub[[base[[i]]]])/y[[i]];
                If[ratio < min,
                   min = ratio; q = i;etato=2]
               ];
         ,{i,m}     
        ];
      If[ub[[p]]-lb[[p]]<min,
          min=ub[[p]]-lb[[p]];etato=0];q
]

etato

indexout=base[[q]]

Which[etato==1
     ,base[[q]]=p;
      If[etafrom==1
        ,eta1[[ppos]]=indexout
        ,eta2=Drop[eta2,{ppos,ppos}];
         eta1=Append[eta1,indexout]
        ]
     ,etato==2
     ,base[[q]]=p;
      If[etafrom==2
        ,eta2[[ppos]]=indexout
        ,eta1=Drop[eta1,{ppos,ppos}];
         eta2=Append[eta2,indexout]
        ]
     ,etato==0
     ,If[etafrom==1
        ,eta1=Drop[eta1,{ppos,ppos}];
         eta2=Append[eta2,p]
        ,eta2=Drop[eta2,{ppos,ppos}];
         eta1=Append[eta1,p]
        ]
     ];

Final solution

solution=Table[0,{n}];
Do[solution[[base[[i]]]]=x[[i]],{i,m}];
Do[solution[[eta1[[i]]]]=lb[[eta1[[i]]]]
            ,{i,Length[eta1]}];
Do[solution[[eta2[[i]]]]=ub[[eta2[[i]]]]
            ,{i,Length[eta2]}];
solution

value=solution.c

resume