Simplex Algorithm

[Albert Chan]

(11-21-2023 05:58 PM)Albert Chan Wrote:  With vars moved to the end, user see [min(c) ... vars]

We might as well put more information in returned list now.
I tracked number of pivot1 calls inside simplex_core(), simplex algorithm is very efficient.

simplex test_cases pivot1 calls.

 1: 3+3 = 6      // use simplex2 wrapper
 2: 3+3 = 6
 3: 3+2 = 5
 4: 2+3 = 5
 5: 3+3 = 6
 6: 2+2 = 4
 7: 2+2 = 4      // simplex_int use simplex2 once
 8: 2+3 = 5
 9: 2+11 = 13
10: 1+1 = 2
11: 2+10 = 12

1st number is simplex2 len(bv) pivot1 calls, to make all pivots=1.
2nd number is simplex_core() total pivot1 calls, to reach solution.

Here is version 5, added all posted ideas since 4

Code:

#cas
simplex(a):=
BEGIN
LOCAL c,v,colJ,I,J,bv,n,(art:=0);
IF dim(dim(a))!=2 THEN a,art:=a END
n:=art;
c,v:=dim(a).-1;  /* constraints, variables */
bv:=makelist(c); /* pivots */
FOR J FROM v DOWNTO 1 DO
 IF n>0 AND J<=v-art THEN error(1,a) END
 colJ:=col(a,J)[1..c];
 IF (I:=contains(colJ,1))==0 THEN continue END
 IF bv[I] or l1norm(colJ[I]:=0) THEN continue END
 bv[I]:=J;
 IF c > art-(n-=1) THEN continue END
 RETURN simplex_core(addartrow(a,art),bv,v,art,0)
END;
error(2,a)  /* too few pivots */
END;

simplex2(a,bv,v,art,ign):=
BEGIN
 LOCAL I;
 FOR I FROM 1 to len(bv) DO a:=pivot1(a,I,bv[I]) END
 RETURN simplex_core(a,bv,v,art,ign)
END;

simplex2(a,bv,v,art,ign):=
BEGIN
 LOCAL I;
 FOR I FROM 1 to len(bv) DO a:=pivot1(a,I,bv[I]) END
 RETURN simplex_core(a,bv,v,art,ign)
END;

simplex_core(a,bv,v,art,ign):=
BEGIN
LOCAL a2,b,c,d,bmin,cmin,n,l1,l2,I,J,P:=0;
d:=dim(a);
c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));
cmin:=MIN(c);
n:=1+(art>0);
b:=SUPPRESS(col(a,0),d(1)-n+1,d(1));
bmin:=MIN(b);

// Step 2: If all ci,bi≥0 go to Step 6
WHILE cmin<0 OR bmin<0 DO

// Step 3: For any ci < 0 check for unbounded solution
 FOR I FROM 1 TO (d(2)-1-ign) DO
  IF c(I) < 0 AND MAX(col(a,I))<=0 THEN
// UNBOUNDED SOLUTION
   RETURN [-inf,a,bv,P,[]];
  END;
 END;

// Step 4: Determine bv to enter and exit
 J:=POS(c,cmin); // xin
 I:=minRatio(col(a,J),col(a,0),n); // xout

// Step 5: Canonical reduction with new bv
 IF bv[I]!=J THEN
  bv[I]:=J;
  a:=pivot1(a,I,J); P+=1;
  c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));
  cmin:=MIN(c);
 END;

 IF cmin>=0 THEN
  IF art>0 THEN
   IF a(d(1),d(2))!=0 THEN
// w≠0? NO SOLUTION / EMPTY FEASIBLE REGION
    RETURN ["Bad Row "+d(1),a,bv,P,[]];
   END;
   a:=delrows(a,d[1]);
   d[1]-=(n:=1);
   ign:=art;
   art:=0;
   c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));
   cmin:=MIN(c);
  ELSE
   b:=SUPPRESS(col(a,0),d(1)-n+1,d(1));
   bmin:=MIN(b);
// IF ∃bi<0, use Dual Simplex
   IF bmin<0 THEN
    FOR I FROM 1 TO len(b) DO
// For any bi < 0 check for unbounded/no solution
     IF b(I)<0 AND MIN(SUPPRESS(row(a,I),d(2)-ign,d(2)))>=0 THEN
      RETURN ["Bad Row "+I,a,bv,P,[]];
     END;
    END;
    I:=POS(b,bmin);
    J:=minRatio(-a[I],c,0); /* max ratio */
    bv[I]:=J;
    a:=pivot1(a,I,J); P+=1;
    c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));
    cmin:=MIN(c);
   END;
  END;
 END;
END;

// Step 6: ci,bi≥0 optimum achieved
l1:=sol(a,v);

a2:=a;
c:=SUPPRESS(row(a2,d(1)),d(2)-ign,d(2));
FOR J FROM 1 TO len(c) DO
 IF c(J) or contains(bv,J) THEN continue END
 IF (b:=abs(trn(col(a2,J))))==1 or b==0 THEN continue END
 I:=minRatio(col(a2,J),col(a2,0),1);
 a2:=pivot1(a2,I,J); P+=1;
 l2:=sol(a2,v);
// IFINITE SOLUTIONS
 IF l1==l2 THEN BREAK END
 RETURN [-a[0][0],a,bv,P,transpose(l1,l2)];
END;

// UNIQUE SOLUTION
RETURN [-a[0][0],a,bv,P,l1];
END;

sol(a,v):=
BEGIN
 LOCAL x,b,c, j,k;
 x:=makelist(v); /* all zeroes */
 b:=col(a,0);    /* last column */
 FOR j FROM 1 TO v DO
  IF abs(trn(c:=col(a,j)))!=1 THEN continue END
  IF (k:=contains(c,1))==0 THEN continue END
  IF b[k] THEN x[j]:=b[k]; b[k]:=0 END
 END;
 RETURN x;
END;

minRatio(a,b,n):=
 BEGIN
 LOCAL i,ratio, (j:=0),(minratio:=inf);
 FOR i FROM 1 TO len(b)-n DO
  IF a[i]<=0 THEN continue END
  ratio:=b[i]/a[i];
  IF minratio>ratio THEN minratio:=ratio; j:=i END
 END;
 RETURN j
END;

// first attempt of Gomory's plane cutting algorithm
simplex_int(a,bv,v,art,ign,integers):=
BEGIN
 LOCAL b,d,x,n;

// Step 1: ignore integer restriction, apply simplex
 a:=simplex2(a,bv,v,art,ign);
 ign:=art;

 WHILE 1 DO
  IF len(x:=a[0])==0 THEN return a END // no solution
  IF len(x[1])>1 THEN return a END     // inf solution
  
  FOR I FROM 1 to len(integers) DO
   IF type(x[integers[I]])!=DOM_INT THEN break END
  END;
  IF I>len(integers) THEN return a END // all integers

// Step 2: Detetmine constraint with largest fractional part
  a,bv,n := a[2..4];
  d:=dim(a);
  b:=col(a,d[2]);
  x := b[1 .. d[1]-1];
  x -= floor(x);
  x := a[POS(x, max(x))];

// Create new constraint
  x:=floor(x)-x;
  I:=len(x)-ign;
  x:=extend(x[1..I-1], [1], x[I..0]);
  bv:=append(bv,d(2)-ign); /* pivot on 1 */
  
// Step 3: solve with new constraint
  REDIM(a,{d(1),d(2)-1});
  ADDCOL(a,b-b,d(2)-ign);
  ADDCOL(a,b,d(2)+1-ign);
  ADDROW(a,x,d(1));
  a:=simplex_core(a,bv,v,0,ign);
  a[4]+=n;  // pivot1 calls
 END;
END;

simplex_le(a) :=
BEGIN
 LOCAL b, p, v,(p:=[]),(d:=inf);
 IF dim(dim(a))!=2 THEN a,d:=a END
 IF abs(d)<len(a) THEN a,p:=remove_eq_c(a,d) END
 b := dim(a);
 a := transpose(extend(a, identity(b[1])));
 a[0] := a[b[2]];
 a := transpose(delrows(a, b[2] .. b[2]+len(p)));
 IF d==0 THEN return a END
 v := sum(b)-2 - len(p);         /* variables */
 p := extend(p,range(b[2],v+1)); /* pivots */
 IF d>0 THEN a[0] := -a[0] END   /* maximize */
 a := simplex_core(a,p,v,0,0);
 IF d>0 THEN a[1] := -a[1] END   /* undo flip */
 return a;
END;

remove_eq_c(a, n) :=
BEGIN
 LOCAL I,J,c := len(a[1]);
 LOCAL p := makelist(n:=abs(n)); /* max pivots */
 FOR I FROM n DOWNTO 1 DO
  p[I] := J := pivot_index(suppress(a[I],c));
  IF J==0 and a[I][c] THEN error(3,a) END
  a := J ? pivot1(a,I,J) : delrows(a,I);
 END
 return a, remove(0,p);
END;

pivot_index(a) :=
BEGIN
 LOCAL b:=contains((a:=abs(a)),1);
 IF b THEN RETURN b END
 IF (b:=max(a)) THEN b:=contains(a,b) END
 RETURN b;
END;

pivot1(m,I,J):= pivot((m[I]:=m[I]/m[I][J]),I,J);
artrow(col,art):=append(makelist(k->k>col-art,2,col),0);

addartrow(a,art):=  /* art vars are valid pivots */
BEGIN
 local c:=len(a[1]);
 IF art<=0 THEN RETURN a END
 append(a,artrow(c,art)-sum(a[POS(col(a,J),1)],J,c-art,c-1))  
END;
#end

simplex_core() is fast, but assume all supplied pivots are good (unit column).
It now returns [min(c), matrix, bv, pivot1_calls, vars]
On HP Prime, most likely, is showed as [min(c) ... vars]

simplex2() is the safer simplex_core(), go thru all the pivots, to make sure all 1's.

simplex() is also very safe, it confirmed a unit column before using it as pivot.
And because it only use true pivot, it does not need pivot1() calls.

Update: added convenience function addartrow(a,art)
It is not simply adding artrow to a, but also to make art variables available to use as pivots.
Note: if art <= 0 it just return a. i.e. no art row is added.

Example:
> m:=[[1,2,0,1,1,0,0,20],[2,1,1,0,0,1,0,10],[-1,4,-2,3,0,0,1,40],[1,4,3,2,0,0,0,0]]
> simplex_core(addartrow(m,3),[5,6,7],7,3,0)

[35 ... [5,0,0,15,0,0,0]]

This is equivalent to simplex2 with plain artrow added, but faster, without pivot1() calls.

> simplex2(append(m,artrow(8,3)),[5,6,7],7,3,0)

[35 ... [5,0,0,15,0,0,0]]

My previous version hide this detail, with artrow added inside simplex_core.
But, to make core restartable at any stage, we need to expose this artrow.

Update2: replaced error msg "No solutions" to more informative "Bad Row " + row_number