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

[Albert Chan]

I added another option to simplex_le, to handle integer programming problems.
Default is maximize (no '=' constraint), and variables not needed to be integers.

simplex_le(a, [k=inf], [z=[]])

Code:

// s.t.  x,  y ≥ 0
//      2x +4y ≥ 22
//      3x +2y ≥ 20
//      4x +5y ≥ 40
// min  6x +5y = z

> m:=[-[2,4,22],-[3,2,20],-[4,5,40],[6,5,0]]
> simplex_le(m, -inf)

[320/7 ... [20/7,40/7,46/7,0,0]]

What if (x,y) required to be non-negative integer? (variables x = 1st, y = 2nd)

> simplex_le(m, -inf, [1,2])

[47 ... [2,7,10,0,3]]

min(z) is higher as expected, solution is x=2, y=7

I am posting only simplex_le wrapper, the rest is from previous post.
With this update, all test_cases can be done with simplex_le()

Code:

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

 WHILE 1 DO     // Gomory's plane cutting algorithm
  IF len(x:=a[0])==0 THEN break END // no solution
  IF len(x[1])>1 THEN break END     // inf solution
  FOR I FROM 1 to len(z) DO
   IF type(x[z[I]])!=DOM_INT THEN break END
  END;
  IF I>len(z) THEN break END        // one solution
// Detetmine constraint with largest frac_part
  a,p,I := tail(a);
  d:=dim(a);
  b:=col(a,d[2]);
  b-=floor(b);
  b[d[1]]:=0;
  x:=a[POS(b, max(b))];
// Create new constraint
  x:=floor(x)-x;
  x[d[2]+1]:=x[d[2]];
  x[d[2]]:=1;
  p[d[1]]:=d[2]; // pivot on 1
// solve with new constraint
  ADDCOL(a,b-b,d[2]);
  ADDROW(a,x,d[1]);
  a:=simplex_core(a,p,v,0,0);
  a[4]+=I;       // pivot1 calls
 END;

 IF k>0 THEN a[1]:=-a[1] END     // undo flip
 RETURN a;
END;

[Albert Chan]

I encountered a bug trying this problem:

Code:

x = (x1,x2,x3) all non-negative integers

2 x1 + 3 x2 + 4 x3 <= 6
4 x1 + 2 x2 + 5 x3 <= 5
22 x1 + 20 x2 + 40 x3 = P, to maximize.

It hang the machine.
Further investigation, it did solve the problem, we just don't know it.
Solution is infinite, if we ignore integer requirements.

x = (0,2,0) --> max(P) = 40
x = (0,0,1) --> max(P) = 40

Current implementation use 0 ≤ t ≤ 1, to denote infinite solutions.
But this is messy to test, perhaps x(t=0) == x(t=1), assuming no numerical issues.

I think better solution is to have x returned as matrix.
Test for infinite solution is simply len(x[1])>1

I have updated version 5 and latest simplex_le() code, here is the result:

> simplex_le([[2,3,4,6],[4,2,5,5],[22,20,40,0]],∞,[1,2,3])

\([40 \;\cdots \left(\begin{array}{cc}
0 & 0 \\
2 & 0 \\
0 & 1 \\
0 & 2 \\
1 & 0
\end{array}\right) ]\)

[ftneek]

The changes all seem good. I'm adding them to the source code.

Just so it's clear, simplex_core now requires canonical form to use while simplex2 allows for standard form input, correct? so simplex2==current simplex_core in source code, and changing the call to simplex2 + updating the output for the current tests should cause them to pass.

I thought you said in one comment that simplex_int would use the simplex() wrapper but I don't think I see that. That would require the equality constraints to be listed first right?

Last, how come you added the plane cutting algorithm again inside simplex_le? Can you not just use simplex_int inside simplex_le?

Albert Chan Wrote:In simplex_core(), it have 2 branch return no solution (P is total pivot1() calls)

RETURN ["No Solution: empty feasible region (w≠0)",a,bv,P,[]];
...
RETURN ["No Solution",a,bv,P,[]];

What is the difference?

The w!=0 check happens when the system has artificial variables and they have been zero'd out of the artificial w row.

The other return statement you are referencing happens during an iteration of dual simplex, and I think should be analogous to the check that happens in step 3. It is possible that the return statements could use a better description.

From my understanding here is what can happen:
X=optimal solution to primal <--> Y=optimal solution to dual
other possibilities:
Max is unbounded above and Dual Min with no feasible solutions
Min is unbounded below and Dual Max with no feasible solutions
Both with no feasible solutions

So I think the w!=0 check might correspond to the 3rd 'other possibility'

[Albert Chan]

(11-23-2023 01:23 AM)ftneek Wrote:  Just so it's clear, simplex_core now requires canonical form to use while simplex2 allows for standard form input, correct? so simplex2==current simplex_core in source code, and changing the call to simplex2 + updating the output for the current tests should cause them to pass.

Yes to all.

Quote:I thought you said in one comment that simplex_int would use the simplex() wrapper but I don't think I see that.

That was in version 4, but switched back in version 5, matching exactly your simplex_int().
The reason is again to simplify updating your current tests.

My preference is simplex_int() using simplex(a, [art=0]) wrapper.
It is faster (simplex wrapper already confirmed unit columns, thus no pivot1 calls),
and less error prone (I had machine lock-up many times, inputting wrong bv values),
and less typing to enter problem.

Quote:Last, how come you added the plane cutting algorithm again inside simplex_le?
Can you not just use simplex_int inside simplex_le?

simplex_int() is not that user friendly, with more arguments required than core!
I moved simplex_int details into simplex_le(), and not use simplex_int() anymore.

2nd reason is efficiency.

Unlike simplex() wrapper, which check for unit columns for pivots, simplex_le does not need all that.
It does not even need to check, because it build the identity matrix!
The matrix is in canonical form, and does not need to go thru simplex2() again.

Also, simplex_le() does not use art/ign variables, code should be much cleaner coded here.
Also, with '=' constraints, matrix produced (less slack vaiables, 0 art variables) is smaller than simplex_int().

Most important of all, less typing to input problem!
No more filling out identity matrix, and counting where the 1's are, for bv.
Especially so after I had to pull artrow outside of simplex_core.

// lpsolve(2x+5y,[3x-y=1,x-y<=5],assume=nonnegint)

> simplex_int([[3,-1,0,1,1],[1,-1,1,0,5],[2,5,0,0,0],[0,0,0,1,0]],{4,3},4,1,0,{1,2})
> simplex_le([[3,-1,1],[1,-1,5],[2,5,0]],-1,[1,2])

[ftneek]

I uploaded new source files with quite a few improvements. The biggest being that simplex() does not require you to manually enter any slack or artificial variables, and it is able to handle maximization problems, integer restrictions, binary restrictions, and unrestricted variables. Also simplex_int() accepts programs in canonical form, to match the change to simplex_core().

Binary example: test16()
> simplex([[5,11,13,16,9,32,32],[7,19,20,23,15,49,62],[9,25,26,29,16,63,77],[9,28,30,35,15,72,88],[-1,1,0,0,0,0,0],[0,0,0,1,1,1,1],[0,0,-1,-1,-1,-1,-1],[10,18,22,25,12,12,0]],0,[],lpmax,[1,2,3,4,5,6])
> [53, ...]

Next I would like to work on simplifying the arguments, so that you can omit an argument if it is not needed. For example currently if there is an unrestricted variable you have to input 6 parameters even if some aren't used or the default value is used. In that sense I think I would like to be a bit more like lpsolve().

Happy holidays!

[Albert Chan]

(12-22-2023 07:38 AM)ftneek Wrote:  I uploaded new source files with quite a few improvements. The biggest being that simplex() does not require you to manually enter any slack or artificial variables, and it is able to handle maximization problems, integer restrictions, binary restrictions, and unrestricted variables. Also simplex_int() accepts programs in canonical form, to match the change to simplex_core().

Wow! That is quite an update.

Unfortunately, this update crash HP Prime emulator 2.1.14181 (2018 10 16)

Start with a hard reset, then close emulator
copy OP *.hpprgm to "HP Prime\Calculators\prime" directory.
start emulator, Cas> test_simplex()      → {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}
restart emulator                                  → CRASHED

Problem seems to arise from simplex.hppgrm, not the others.
And, it is not because of binary nature of .hppgrm

Start with a hard reset
SHIFT PROGRAM NEW, name=simplex, copy actual #cas ... #end code
Cas>lpmax --> ∞
restart emulator ... failed
restart emulator, Cas> lpmax --> lpmax
restart emulator, Cas> lpmax --> ∞
restart emulator ... failed
restart emulator, Cas> lpmax --> lpmax
restart emulator, Cas> lpmax --> ∞
...

[ftneek]

(12-22-2023 04:07 PM)Albert Chan Wrote:  Wow! That is quite an update.

Unfortunately, this update crash HP Prime emulator 2.1.14181 (2018 10 16)

Strange, though I also encountered an issue.

When I took a break yesterday and came back to work on simplex() I noticed the file no longer compiled on my virtual calculator (2020 01 16), it gave a syntax error on a part of the code I didn't change earlier. I uploaded the same files to my physical calculator (2023 04 13), simplex.hpprgm compiled and test_simplex() passed. On the virtual calculator I clicked Calculator > Reset and transferred the files again, they complied successfully.

I'm not sure what caused it but the reset seems to have resolved the syntax error that occurred on my virtual calculator.
Edit: but now I see what you mean about crashing the virtual calculator... I'm not sure what could be causing it, but I am open to solutions.

[Albert Chan]

(12-22-2023 04:07 PM)Albert Chan Wrote:  Problem seems to arise from simplex.hppgrm, not the others.

Problem is simplex() first CASE ... END, to parse user input.
Removing first CASE ... END solved crashing issue.

I simplified a bit, with extended simplex_le(a, [dir=inf], [integers]) pattern.
Slight difference, simplex() wrapper default is minimize.

Code:

// simplex(a,[dir],[integers],[binary],[unrestricted])
simplex(a):=
BEGIN
LOCAL (binary:=0),(integers:=0),(unrestricted:=0),(art:=0),(dir:=-inf);
LOCAL b,c,d,v,bv,n,I,J;
IF dim(d:=dim(a))!=2 THEN
 IF d>5 THEN error(a) END;
 a,dir,integers,binary,unrestricted := a;
 IF (art:=abs(dir))==inf THEN art:=0 END;
 d:=dim(a);
END;

IF unrestricted!=0 THEN
 FOR I FROM len(unrestricted) DOWNTO 1 DO
  b:=col(a,unrestricted[I]);
  a:=ADDCOL(a,-b,unrestricted[I]+1);
  n:=contains(integers,unrestricted[I]);
  IF n THEN
   integers:=REDIM(integers,len(integers)+1);
   FOR J FROM len(integers) DOWNTO n+1 DO integers[J]:=integers[J-1]+1 END;
  END;
 END;
d[2]+=len(unrestricted);
END;

c:=row(a,d[1]);
a:=delrow(a,d[1]);

IF binary!=0 THEN
 v:=[];
 FOR I FROM 1 TO len(binary) DO
  n:=makelist(d[2]);
  n[binary[I]]:=1;
  n[0]:=1;
  v:=append(v,n);
 END;
 a:=append(a,op(v));
 d[1]+=len(binary);
 integers:=UNION(integers,binary);
END;

b:=col(a,d[2]);
a:=delcol(a,d[2]);
d:=d.-1;

CASE
 IF art==0 OR art==d[1] THEN a:=extend(a,IDENMAT(d[1])) END;
 DEFAULT
  v:=append(makemat(0,art,d[1]-art),op(IDENMAT(d[1]-art)));
  v:=extend(v,append(IDENMAT(art),op(makemat(0,d[1]-art,art))));
  a:=extend(a,v);
END;
a:=ADDCOL(a,b,d[2]+d[1]+1);
c:=REDIM(c,len(c)+d[1]);
c[0]:=c[d[2]+1];
c[d[2]+1]:=0;
a:=append(a,c);

n:=art;
c,v:=dim(a).-1;  /* constraints, variables */
bv:=range(v-art+1,v+1);
bv:=CONCAT(bv,range(d[2]+1,v-art+1));
IF dir>0 THEN a[0]:=-a[0] END;
IF integers != 0 THEN
 c:=simplex_int(addartrow(a,art),bv,v,art,0,integers)
ELSE
 c:=simplex_core(addartrow(a,art),bv,v,art,0);
END;
IF dir>0 THEN c[1]:=-c[1] END;
RETURN c
END;

(11-22-2023 06:44 PM)Albert Chan Wrote:  > m:=[-[2,4,22],-[3,2,20],-[4,5,40],[6,5,0]]
> simplex_le(m, -inf)

[320/7 ... [20/7,40/7,46/7,0,0]]

What if (x,y) required to be non-negative integer? (variables x = 1st, y = 2nd)

> simplex_le(m, -inf, [1,2])

[47 ... [2,7,10,0,3]]

min(z) is higher as expected, solution is x=2, y=7

> m:=[-[2,4,22],-[3,2,20],-[4,5,40],[6,5,0]]
> simplex(m); /* simplex default = minimize */

\([\frac{320}{7},\left(\begin{array}{cccccc}
0 & 0 & 1 & \frac{6}{7} & \frac{-8}{7} & \frac{46}{7} \\
1 & 0 & 0 & \frac{-5}{7} & \frac{2}{7} & \frac{20}{7} \\
0 & 1 & 0 & \frac{4}{7} & \frac{-3}{7} & \frac{40}{7} \\
0 & 0 & 0 & \frac{10}{7} & \frac{3}{7} & \frac{-320}{7}
\end{array}\right) ,[3,1,2],2,\left(\begin{array}{c}
\frac{20}{7} \\
\frac{40}{7} \\
\frac{46}{7} \\
0 \\
0
\end{array}\right) ]\)

> simplex(m, -inf, [1,2])

\([47,\left(\begin{array}{ccccccc}
0 & 0 & 1 & 2 & 0 & -4 & 10 \\
1 & 0 & 0 & -1 & 0 & 1 & 2 \\
0 & 1 & 0 & 1 & 0 & \frac{-3}{2} & 7 \\
0 & 0 & 0 & 1 & 1 & \frac{-7}{2} & 3 \\
0 & 0 & 0 & 1 & 0 & \frac{3}{2} & -47
\end{array}\right) ,[3,1,2,5],3,\left(\begin{array}{c}
2 \\
7 \\
10 \\
0 \\
3
\end{array}\right) ]\)

Note that variables list now returned as a matrix.

Also, simplex() solved equality constraints using art variables, simplex_le() by variable eliminations.

test01() example:

> simplex([[1,2,0,1,20],[2,1,1,0,10],[-1,4,-2,3,40],[1,4,3,2,0]],-3)

\([35,\left(\begin{array}{cccccccc}
1 & \frac{1}{2} & \frac{1}{2} & 0 & \frac{3}{4} & 0 & \frac{-1}{4} & 5 \\
0 & 0 & 0 & 0 & \frac{-3}{2} & 1 & \frac{1}{2} & 0 \\
0 & \frac{3}{2} & \frac{-1}{2} & 1 & \frac{1}{4} & 0 & \frac{1}{4} & 15 \\
0 & \frac{1}{2} & \frac{7}{2} & 0 & \frac{-5}{4} & 0 & \frac{-1}{4} & -35
\end{array}\right) ,[1,6,4],3,\left(\begin{array}{c}
5 \\
0 \\
0 \\
15 \\
0 \\
0 \\
0
\end{array}\right) ]\)

> simplex_le([[1,2,0,1,20],[2,1,1,0,10],[-1,4,-2,3,40],[1,4,3,2,0]],-3)

\([35,\left(\begin{array}{ccccc}
0 & \frac{3}{2} & \frac{-1}{2} & 1 & 15 \\
1 & \frac{1}{2} & \frac{1}{2} & 0 & 5 \\
0 & \frac{1}{2} & \frac{7}{2} & 0 & -35
\end{array}\right) ,[4,1],3,\left(\begin{array}{c}
5 \\
0 \\
0 \\
15
\end{array}\right) ]\)

[ftneek]

(12-23-2023 04:44 AM)Albert Chan Wrote:  Problem is simplex() first CASE ... END, to parse user input.
Removing first CASE ... END solved crashing issue.

Thanks for catching the source of the crash/reset, I implemented the change and uploaded new files.

(12-23-2023 04:44 AM)Albert Chan Wrote:  Also, simplex() solved equality constraints using art variables, simplex_le() by variable eliminations.

test01() example:

> simplex([[1,2,0,1,20],[2,1,1,0,10],[-1,4,-2,3,40],[1,4,3,2,0]],-3)
.
.
.
> simplex_le([[1,2,0,1,20],[2,1,1,0,10],[-1,4,-2,3,40],[1,4,3,2,0]],-3)

Is simplex_le able to solve test14? It is problem 2d in section 6.3 of the referenced textbook.

> simplex([[1,-1,0,12],[0,2,3,28],[2,-4,1,0]],2,[1,2,3])
> [28, ... , trn([14,2,8,...])]

> simplex_le([[1,-1,0,12],[0,2,3,28],[2,-4,1,0]],2,[1,2,3])
> infinte loop

[Albert Chan]

(12-23-2023 07:46 AM)ftneek Wrote:  Is simplex_le able to solve test14? It is problem 2d in section 6.3 of the referenced textbook.

Yes. We need a patch for variable list now returned as matrix.
simplex_le() also patched to counted *all* pivot1() used, matching simplex() behavior.

> simplex([[1,-1,0,12],[0,2,3,28],[2,-4,1,0]],2,[1,2,3])

\([28,\left(\begin{array}{cccccccc}
1 & 0 & 0 & 0 & -3 & 1 & 2 & 14 \\
0 & 0 & 1 & 0 & 2 & 0 & -1 & 8 \\
0 & 1 & 0 & 0 & -3 & 0 & 2 & 2 \\
0 & 0 & 0 & 1 & -2 & 0 & 1 & 1 \\
0 & 0 & 0 & 0 & 8 & 2 & -5 & 28
\end{array}\right) ,[1,3,2,4],4,\left(\begin{array}{c}
14 \\
2 \\
8 \\
1 \\
0
\end{array}\right) ]\)

> simplex_le([[1,-1,0,12],[0,2,3,28],[2,-4,1,0]],2,[1,2,3])

\([28,\left(\begin{array}{cccccc}
1 & 0 & 0 & 0 & -3 & 14 \\
0 & 0 & 1 & 0 & 2 & 8 \\
0 & 1 & 0 & 0 & -3 & 2 \\
0 & 0 & 0 & 1 & -2 & 1 \\
0 & 0 & 0 & 0 & 8 & 28
\end{array}\right) ,[1,3,2,4],4,\left(\begin{array}{c}
14 \\
2 \\
8
\end{array}\right) ]\)

I wil sent you updated code via PM, for consideration.

[ftneek]

I incorporated your changes and uploaded new files. Also now simplex_core returns ∅ for infeasible model with no solutions rather than "No solution" or the empty feasible region string, they all mean the same thing. If you want to avoid negating for infeasible max in simplex_le, the check is just

Code:

IF dir>0 AND c[1]!=∅ THEN c[1]:=-c[1] END;

[Albert Chan]

(11-22-2023 10:10 PM)Albert Chan Wrote:  > simplex_le([[2,3,4,6],[4,2,5,5],[22,20,40,0]],∞,[1,2,3])

\([40 \;\cdots \left(\begin{array}{cc}
0 & 0 \\
2 & 0 \\
0 & 1 \\
0 & 2 \\
1 & 0
\end{array}\right) ]\)

simplex_le() simply return infinite solutions (integers or not).
It is just luck that above variables happened to be integers.

This fixed it.
simplex_le(), before patch

Code:

 WHILE 1 DO     // Gomory's plane cutting algorithm
  IF len(x:=a[0])==0 THEN break END // no solution
  IF len(x[1])>1 THEN break END     // inf solution
  FOR I FROM 1 to len(z) DO
   IF type(x[z[I]][1])!=DOM_INT THEN break END
  END;
  IF I>len(z) THEN break END        // one solution

simplex_le(), after patch

Code:

 WHILE len(x:=a[0]) DO  // Gomory's plane cutting algorithm
  FOR b FROM 1 to len(x[1]) DO
   FOR I FROM 1 to len(z) DO
    IF type(x[z[I]][b])!=DOM_INT THEN break END
   END;
   IF I>len(z) THEN break END // solution
  END;
  IF I>len(z) THEN break END  // solution

break 2 unable to get out of WHILE loop (why?), thus the ugly multiple breaks.
If there is a better way, don't hesitate to post!

[ftneek]

(12-24-2023 03:46 PM)Albert Chan Wrote:  break 2 unable to get out of WHILE loop (why?), thus the ugly multiple breaks.
If there is a better way, don't hesitate to post!

Your first for loop seems to iterate from 1 to len(x[1]), we only need to check the indices from the integers list, anything else does not matter if it is integer or not.

Also that break statement is nested inside a while and 2 for loops, so you may need to use break 3. Or you are checking I>len(z) twice, should the second one be b>len(z)?

[Albert Chan]

(12-24-2023 07:32 PM)ftneek Wrote:  Your first for loop seems to iterate over all x[I] (from 1 to len(x[1])), we only need to check the indices from the integers list, anything else does not matter if it is integer or not.

It was simply simplex_le original code, nested with b loop to handle possibly 2 edges.

1st b loop select column of x = 1 .. len(x[1]) = 1 .. (1 or 2)
2nd I loop select rows from integer list z, *not* all rows of x

Quote:break statement is nested inside a while and 2 for loops, so you may need to use break 3.
Or you are checking I>len(z) twice, should the second one be b>len(z)?

I mean to break out from b loop and WHILE loop.

But, with break n bug, this does not work.
It does not matter if you change 2 to 3, or even break 99

Code:

 WHILE len(x:=a[0]) DO  // Gomory's plane cutting algorithm
  FOR b FROM 1 to len(x[1]) DO
   FOR I FROM 1 to len(z) DO
    IF type(x[z[I]][b])!=DOM_INT THEN break END
   END;
   IF I>len(z) THEN break 2 END // 2 get ignored!
  END;

[ftneek]

Uploaded latest changes which returns "∅" instead of ∅ for infeasible model and switches pivoting rule to Bland's rule when a repeated basis is detected.

Appendix B: An example of cycling / test18()
> simplex([[1,0,0,1/4,-8,-1,9,0],[0,1,0,1/2,-12,-1/2,3,0],[0,0,1,0,0,1,0,1],[0,0,0,-3/4,20,-1/2,6,0]],-3)

[-5/4,..., transpose([3/4,0,0,1,0,1,0,0,0,0])]

I tried to verify the solution using other solvers (xcas lpsolve, Wolfram alpha..) but had no luck so far. One thing I noticed is that the appendix says the repeated basis is detected on the 6th iteration, but according to the P counter it was detected on the 8th pivot1 call (see the commented out return on line 185), 2 tableaus passed the initial tableau shown.

The hashin file uses the new program structure (part PPL and part python) to access the python hash function. The downside is that this causes the outdated connectivity kit on my main computer (mac) to crash upon trying to open the file, but it is simple enough it can probably be split up or included in the simplex file to still work on older software. Hopefully updates will come soon for outdated software, and the new program structure can actually be used.

[ftneek]

Updated Nov 17, 2024:
The hashin file has been replaced with basisID. There is no longer a dependency on micro python or the new program structure, so this code should work fine on older VCs and CKs. The basisID file adds the commands basis_to_id and id_to_basis. I also packaged all the files as a zip.

Previously in hashin, I relied on the micro python hash function to detect a repeated basis. I recently learned of a method to map each basis to an id and that the basis can be recovered from the id, which has a few applications for linear programming, including to detect cycling. Let me know if any issues are detected.

Example:
basis=[3,4,5]
total # of variables=5
3 variables in basis

basis_to_id([3,4,5],5) -> 9
id_to_basis(9,5,3) -> [3,4,5]

[ftneek]

Nov 26, 2024
Most importantly, there was a bug fix affecting some integer problems with = constraints.

Besides that, infeasible models will now return a localized string depending on the current language that translates to "No solutions" (thanks to the STRINGFROMID command). NO_SOLUTIONS() will return the string for the current language.

Additionally, several tests were updated for one reason or another (bug fix, incorrect input, internal adjustments).

Lastly, all PPL commands (besides STRINGFROMID) have been migrated to equivalent CAS commands. Next, I would like to make adjustments to add support for 0 based indexing, which I'm hoping would allow the same code to work on the HP Prime as well as inside Xcas.