Simplex Algorithm
11-21-2023, 11:20 PM (This post was last modified: 11-23-2023 07:48 PM by Albert Chan.)
Post: #41
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
(11-21-2023 05:58 PM)Albert Chan Wrote:  With vars moved to the end, user see [min(c) ... vars]

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.

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]]

[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.

11-22-2023, 06:44 PM (This post was last modified: 11-24-2023 12:40 PM by Albert Chan.)
Post: #42
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
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;
11-22-2023, 10:10 PM
Post: #43
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
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) ]$$
11-23-2023, 01:23 AM (This post was last modified: 11-23-2023 01:53 AM by ftneek.)
Post: #44
 ftneek Member Posts: 74 Joined: Oct 2022
RE: Simplex Algorithm
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'

- neek
11-23-2023, 06:35 AM
Post: #45
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
(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])
12-22-2023, 07:38 AM (This post was last modified: 12-22-2023 08:08 AM by ftneek.)
Post: #46
 ftneek Member Posts: 74 Joined: Oct 2022
RE: Simplex Algorithm
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!

- neek
12-22-2023, 04:07 PM
Post: #47
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
(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)

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

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 --> ∞
...
12-22-2023, 07:28 PM (This post was last modified: 12-22-2023 11:36 PM by ftneek.)
Post: #48
 ftneek Member Posts: 74 Joined: Oct 2022
RE: Simplex Algorithm
(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.

- neek
12-23-2023, 04:44 AM (This post was last modified: 12-23-2023 02:21 PM by Albert Chan.)
Post: #49
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
(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) ]$$
12-23-2023, 07:46 AM
Post: #50
 ftneek Member Posts: 74 Joined: Oct 2022
RE: Simplex Algorithm
(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

- neek
12-23-2023, 10:23 AM
Post: #51
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
(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.
12-23-2023, 09:30 PM
Post: #52
 ftneek Member Posts: 74 Joined: Oct 2022
RE: Simplex Algorithm
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;

- neek
12-24-2023, 03:46 PM
Post: #53
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
(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!
12-24-2023, 07:32 PM (This post was last modified: 12-24-2023 07:40 PM by ftneek.)
Post: #54
 ftneek Member Posts: 74 Joined: Oct 2022
RE: Simplex Algorithm
(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)?

- neek
12-24-2023, 08:05 PM
Post: #55
 Albert Chan Senior Member Posts: 2,516 Joined: Jul 2018
RE: Simplex Algorithm
(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;
01-07-2024, 02:40 AM (This post was last modified: 01-07-2024 06:13 AM by ftneek.)
Post: #56
 ftneek Member Posts: 74 Joined: Oct 2022
RE: Simplex Algorithm
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.

- neek
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)