Simplex Algorithm
11-20-2023, 06:12 PM (This post was last modified: 11-21-2023 08:08 PM by Albert Chan.)
Post: #36
 Albert Chan Senior Member Posts: 2,637 Joined: Jul 2018
RE: Simplex Algorithm
(11-19-2023 04:58 PM)Albert Chan Wrote:  I am not convinced that -1 slack variable will cause problem.

I am now convinced that it will not cause problem, even if scaled.
This is because minRatio and maxRatio calculated in completely opposite way.
minRatio(a,b) skipped over a≤0, while maxRatio skipped over a≥0
As long as pivot is non-zero, we are safe, the mistaken pivot will be resolved.

Anyway, we can always grab a +1 pivot, there is no reason to worry about -1
≤ constraint: slack variable of 1
= constraint: art variable of 1
≥ constraint: slack of -1, art of +1

I also noticed user cannot just grab any cell as pivot, it may produce wrong result.
Thus, a wrapper is very important, disallowed randomly picking pivots.
If we do want (i,j) as pivot, do a:=pivot1(a,i,j) before picking the next one.

I am reverting back to the fast method, without (again) canonicalize all pivots.
And, no need to fill artrow, even in simplex_core(), it handles this internally

simplex(m, [art=0]) ≡ simplex_core(m, pivots, variables, art, 0)

Code:
IF art THEN /* add artrow */  a:=append(a,artrow((c:=len(a[1])),art));  a[0]-=sum(a[POS(bv,J)],J,c-art,c-1); END;

With the same logic, simplex_int() now use simplex wrapper.
test_case_7(), before and after. Less typing, less errors.

< 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_int([[3,-1,0,1,1],[1,-1,1,0,5],[2,5,0,0,0]],1,{1,2})

Or course, the tests now need adjustments.
I am including both code and tests. test_simplex() now run 20% faster.

> test_simplex()      → {1,1,1,1,1,1,1,1,1,1,1}

Lets call this version 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(a,bv,v,art,0) END; error(2,a)  /* too few pivots */ END; simplex_core(a,bv,v,art,ign):= BEGIN LOCAL a2,b,c,d,bmin,cmin,n,l1,l2,I,J; IF art THEN /* add artrow */  a:=append(a,artrow((c:=len(a[1])),art));  a[0]-=sum(a[POS(bv,J)],J,c-art,c-1); END; 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];   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);   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 ["No Solution: empty feasible region (w≠0)",[],a,bv];    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 ["No Solution",[],a,bv];      END;     END;     I:=POS(b,bmin);     J:=minRatio(-a[I],c,0); /* max ratio */     bv[I]:=J;     a:=pivot1(a,I,J);     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);  l2:=sol(a2,v); // IFINITE SOLUTIONS  IF l1==l2 THEN BREAK END  RETURN [-a[0][0],l2-(l2-l1)*t,a,bv]; END; // UNIQUE SOLUTION RETURN [-a[0][0],l1,a,bv]; 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,art,integers):= BEGIN  LOCAL b,d,x,m,I,bv,ign,v:=len(a[1]); // Step 1: ignore integer restriction, apply simplex  a:=simplex(a,art);  ign:=art;  WHILE 1 DO   IF len(x:=a[2])==0 THEN return a END /* no 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 := a[3..4];   d:=dim(a);   b:=col(a,d[2]);   m := b[1 .. d[1]-1];   m -= floor(m);   m := a[POS(m, max(m))]; // Create new constraint   m:=floor(m)-m;   I:=len(m)-ign;   m:=extend(m[1..I-1], [1], m[I..0]);   REDIM(a,{d(1),d(2)-1});   ADDCOL(a,b-b,d(2)-ign);   ADDCOL(a,b,d(2)+1-ign);   ADDROW(a,m,d(1));   bv:=append(bv,d(2)-ign); /* pivot on 1 */ // Step 3: solve with new constraint   a:=simplex_core(a,bv,v,0,ign);  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,j,k):= pivot((m[j]:=m[j]/m[j][k]),j,k); artrow(col,art):=append(makelist(k->k>col-art,2,col),0); #end

Code:
#cas test_simplex():= BEGIN RETURN {test_case_1(),test_case_2(),test_case_3(),test_case_4(),test_case_5(), test_case_6(),test_case_7(),test_case_8(),test_case_9(),test_case_10(), test_case_11()}; END; test_case_1():= BEGIN // Example 3.7.2 // artificial variables // s.t.  x1,  x2,  x3,  x4 ≥ 0 //       x1 +2x3      + x4 = 20 //      2x1 + x2 + x3      = 10 //      -x1 +4x2 -2x3 +3x4 = 40 // min   x1 +4x2 +3x3 +2x4 = z LOCAL expected,result; expected:=[35,[5,0,0,15,0,0,0],[[1,1/2,1/2,0,3/4,0,-1/4,5],[0,0,0,0,-3/2,1,1/2,0],[0,3/2,-1/2,1,1/4,0,1/4,15],[0,1/2,7/2,0,-5/4,0,-1/4,-35]],{1,6,4}]; result:=simplex_core([[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]],{5,6,7},7,3,0); RETURN sum(abs(expected - result)) == 0; END; test_case_2():= BEGIN // -slack variables, artificial variables // s.t.  x,  y ≥ 0 //      2x +4y ≥ 22 //      3x +2y ≥ 20 //      4x +5y ≥ 40 // min  6x +5y = z LOCAL expected,result; expected:=[320/7,[20/7,40/7,46/7,0,0,0,0,0],[[0,1,0,4/7,-3/7,0,-4/7,3/7,40/7],[1,0,0,-5/7,2/7,0,5/7,-2/7,20/7],[0,0,1,6/7,-8/7,-1,-6/7,8/7,46/7],[0,0,0,10/7,3/7,0,-10/7,-3/7,-320/7]],{2,1,3}]; result:=simplex_core([[2,4,-1,0,0,1,0,0,22],[3,2,0,-1,0,0,1,0,20],[4,5,0,0,-1,0,0,1,40],[6,5,0,0,0,0,0,0,0]],{6,7,8},8,3,0); RETURN sum(abs(expected - result)) == 0; END; test_case_3():= BEGIN // slack varaibles // s.t.  x,  y ≥ 0 //     -2x -4y ≤-22 //     -3x -2y ≤-20 //     -4x -5y ≤-40 // min  6x +5y = z LOCAL expected,result; expected:=[320/7,[20/7,40/7,46/7,0,0],[[0,0,1,6/7,-8/7,46/7],[1,0,0,-5/7,2/7,20/7],[0,1,0,4/7,-3/7,40/7],[0,0,0,10/7,3/7,-320/7]],{3,1,2}]; result:=simplex_core([[-2,-4,1,0,0,-22],[-3,-2,0,1,0,-20],[-4,-5,0,0,1,-40],[6,5,0,0,0,0]],{3,4,5},5,0,0); RETURN sum(abs(expected - result)) == 0; END; test_case_4():= BEGIN // s.t.  x,  y,  z ≥ 0 //      3x +2y + z = 10 //      2x +5y +3z = 15 // min -2x -3y -4z = Z LOCAL expected,result; expected:=[-130/7,[15/7,0,25/7,0,0],[[1,1/7,0,3/7,-1/7,15/7],[0,11/7,1,-2/7,3/7,25/7],[0,25/7,0,-2/7,10/7,130/7]],{1,3}]; result:=simplex_core([[3,2,1,1,0,10],[2,5,3,0,1,15],[-2,-3,-4,0,0,0]],{4,5},5,2,0); RETURN sum(abs(expected - result)) == 0; END; test_case_5():= BEGIN // artificial variables, non zero offset // s.t.  x1,  x2,  x3,  x4,  x5,  x6,  x7 ≥ 0 //      4x1 - x2       +x4 - x5 + x6          = 6 //     -7x1 +8x2 + x3                - x7     = 7 //       x1 + x2      +4x4 -4x5               = 12 // min -3x1 +2x2 + x3 - x4 + x5           +87 = z LOCAL expected,result; expected:=[1441/17,[131/85,189/85,0,35/17,0,0,0,0,0,0],[[1,0,1/17,0,0,32/85,-1/17,32/85,1/17,-8/85,131/85],[0,1,3/17,0,0,28/85,-3/17,28/85,3/17,-7/85,189/85],[0,0,-1/17,1,-1,-3/17,1/17,-3/17,-1/17,5/17,35/17],[0,0,13/17,0,0,5/17,4/17,5/17,-4/17,3/17,-1441/17]],{1,2,4}]; result:=simplex_core([[4,-1,0,1,-1,1,0,1,0,0,6],[-7,8,1,0,0,0,-1,0,1,0,7],[1,1,0,4,-4,0,0,0,0,1,12],[-3,2,1,-1,1,0,0,0,0,0,-87]],{8,9,10},10,3,0); RETURN sum(abs(expected - result)) == 0; END; test_case_6():= BEGIN // slack variables 'dual problem' // s.t.  z1,z2,z3 ≥ 0 //      2z1 + 3z2 + 4z3 ≤ 6 //      4z1 + 2z2 + 5z3 ≤ 5 // max 22z1 +20z2 +40z3 = c LOCAL expected,result; expected:=[-320/7,[0,10/7,3/7,0,0],[[-6/7,1,0,5/7,-4/7,10/7],[8/7,0,1,-2/7,3/7,3/7],[46/7,0,0,20/7,40/7,320/7]],{2,3}]; result:=simplex_core([[2,3,4,1,0,6],[4,2,5,0,1,5],[-22,-20,-40,0,0,0]],{4,5},5,0,0); RETURN sum(abs(expected - result)) == 0; END; test_case_7():= BEGIN // lpsolve(2x+5y,[3x-y=1,x-y<=5],assume=nonnegint) LOCAL expected,result; expected:=[12,[1,2,1,0,0],[[1,0,0,-1,1,1],[0,0,1,-2,6,1],[0,1,0,-3,2,2],[0,0,0,17,-12,-12]],{1,3,2}]; result:=simplex_int([[3,-1,0,1,1],[1,-1,1,0,5],[2,5,0,0,0]],1,{1,2}); RETURN sum(abs(expected - result)) == 0; END; test_case_8():= BEGIN // lpsolve(x+3y+3z,[x+3y+2z<=7,2x+2y+z<=11], assume=lp_nonnegative,lp_maximize, lp_integervariables=[x,z]) LOCAL expected,result; expected:=[-10,[1,0,3,0,6,0],[[0,1,1,0,0,1,3],[0,-1,0,-2,1,3,6],[1,1,0,1,0,-2,1],[0,1,0,1,0,1,10]],{3,5,1}]; result:=simplex_int([[1,3,2,1,0,7],[2,2,1,0,1,11],[-1,-3,-3,0,0,0]],0,{1,3}); RETURN sum(abs(expected - result)) == 0; END; test_case_9():= BEGIN // lpsolve(-5x-7y,[7x+y<=35,-x+3y<=6],assume=integer) LOCAL expected,result; expected:=[-41,[4,3,4,1,1],[[1,0,0,0,0,0,0,0,0,0,-5/12,1/4,4],[0,1,0,0,0,0,0,0,0,0,1/3,0,3],[0,0,0,0,0,0,1,0,0,0,-1/6,-1/2,1],[0,0,0,0,1,0,0,0,0,0,-1,0,1],[0,0,0,0,0,1,0,0,0,0,-7/12,-1/4,1],[0,0,1,0,0,0,0,0,0,0,31/12,-7/4,4],[0,0,0,0,0,0,0,1,0,0,9/4,-7/4,3],[0,0,0,0,0,0,0,0,1,0,2/3,-1,1],[0,0,0,1,0,0,0,0,0,0,-17/12,1/4,1],[0,0,0,0,0,0,0,0,0,1,5/12,-5/4,1],[0,0,0,0,0,0,0,0,0,0,1/4,5/4,41]],{1,2,7,5,6,3,8,9,4,10}]; result:=simplex_int([[7,1,1,0,35],[-1,3,0,1,6],[-5,-7,0,0,0]],0,{1,2}); RETURN sum(abs(expected - result)) == 0; END; test_case_10():= BEGIN // max 3x1+2x2 // s.t. x1+x2<=6 // x1,x2>=0 and integer LOCAL expected,result; expected:=[-18,[6,0,0],[[1,1,1,6],[0,1,3,18]],[1]]; result:=simplex_int([[1,1,1,6],[-3,-2,0,0]],0,[1,2]); RETURN sum(abs(expected - result)) == 0; END; test_case_11():= BEGIN // max 3x1+2x2 // s.t. x1+x2<=6 // x1,x2>=0 and integer LOCAL expected,result; expected:=[-58,[2,4,0,92,90],[[0,1,1,0,0,0,0,0,0,0,0,-1,4],[1,0,-4,0,0,0,0,0,0,0,0,9/2,2],[0,0,52,1,0,0,0,0,0,0,0,-115/2,92],[0,0,51,0,1,0,0,0,0,0,0,-113/2,90],[0,0,50,0,0,1,0,0,0,0,0,-111/2,88],[0,0,48,0,0,0,1,0,0,0,0,-107/2,84],[0,0,46,0,0,0,0,1,0,0,0,-103/2,80],[0,0,4,0,0,0,0,0,1,0,0,-9/2,6],[0,0,3,0,0,0,0,0,0,1,0,-7/2,4],[0,0,2,0,0,0,0,0,0,0,1,-5/2,2],[0,0,1,0,0,0,0,0,0,0,0,1/2,58]],[2,1,4,5,6,7,8,9,10,11]]; result:=simplex_int([[2,9,1,0,40],[11,-8,0,1,82],[-3,-13,0,0,0]],0,[1,2]); RETURN sum(abs(expected - result)) == 0; END; #end
 « Next Oldest | Next Newest »

 Messages In This Thread Simplex Algorithm - ftneek - 11-11-2023, 11:21 PM RE: Simplex Algorithm - Albert Chan - 11-11-2023, 11:51 PM RE: Simplex Algorithm - Albert Chan - 11-12-2023, 12:38 AM RE: Simplex Algorithm - Albert Chan - 11-12-2023, 05:11 PM RE: Simplex Algorithm - ftneek - 11-12-2023, 05:40 PM RE: Simplex Algorithm - Albert Chan - 11-13-2023, 02:48 PM RE: Simplex Algorithm - ftneek - 11-12-2023, 12:00 AM RE: Simplex Algorithm - ftneek - 11-12-2023, 01:44 AM RE: Simplex Algorithm - Albert Chan - 11-12-2023, 11:04 PM RE: Simplex Algorithm - ftneek - 11-13-2023, 01:51 AM RE: Simplex Algorithm - Albert Chan - 11-12-2023, 08:46 PM RE: Simplex Algorithm - ftneek - 11-12-2023, 10:09 PM RE: Simplex Algorithm - ftneek - 11-13-2023, 05:34 PM RE: Simplex Algorithm - Albert Chan - 11-13-2023, 10:46 PM RE: Simplex Algorithm - Albert Chan - 11-14-2023, 02:25 AM RE: Simplex Algorithm - ftneek - 11-14-2023, 06:24 AM RE: Simplex Algorithm - Albert Chan - 11-14-2023, 09:23 AM RE: Simplex Algorithm - ftneek - 11-14-2023, 10:39 AM RE: Simplex Algorithm - Albert Chan - 11-15-2023, 04:36 PM RE: Simplex Algorithm - ftneek - 11-15-2023, 05:58 PM RE: Simplex Algorithm - Albert Chan - 11-16-2023, 11:42 AM RE: Simplex Algorithm - Albert Chan - 11-16-2023, 07:15 PM RE: Simplex Algorithm - ftneek - 11-17-2023, 07:47 AM RE: Simplex Algorithm - Albert Chan - 11-17-2023, 11:04 AM RE: Simplex Algorithm - ftneek - 11-18-2023, 01:27 AM RE: Simplex Algorithm - ftneek - 11-18-2023, 10:31 PM RE: Simplex Algorithm - Albert Chan - 11-19-2023, 12:57 AM RE: Simplex Algorithm - ftneek - 11-19-2023, 07:05 AM RE: Simplex Algorithm - Albert Chan - 11-19-2023, 04:58 PM RE: Simplex Algorithm - Albert Chan - 11-20-2023 06:12 PM RE: Simplex Algorithm - ftneek - 11-21-2023, 08:36 AM RE: Simplex Algorithm - Albert Chan - 11-21-2023, 02:05 PM RE: Simplex Algorithm - ftneek - 11-19-2023, 08:02 PM RE: Simplex Algorithm - ftneek - 11-20-2023, 12:18 AM RE: Simplex Algorithm - Albert Chan - 11-20-2023, 02:14 AM RE: Simplex Algorithm - ftneek - 11-20-2023, 09:02 AM RE: Simplex Algorithm - Albert Chan - 11-20-2023, 11:42 AM RE: Simplex Algorithm - Albert Chan - 11-20-2023, 03:34 PM RE: Simplex Algorithm - ftneek - 11-20-2023, 07:52 PM RE: Simplex Algorithm - Albert Chan - 11-21-2023, 05:58 PM RE: Simplex Algorithm - Albert Chan - 11-21-2023, 11:20 PM RE: Simplex Algorithm - Albert Chan - 11-22-2023, 06:44 PM RE: Simplex Algorithm - Albert Chan - 11-22-2023, 10:10 PM RE: Simplex Algorithm - Albert Chan - 12-24-2023, 03:46 PM RE: Simplex Algorithm - ftneek - 12-24-2023, 07:32 PM RE: Simplex Algorithm - Albert Chan - 12-24-2023, 08:05 PM RE: Simplex Algorithm - ftneek - 11-23-2023, 01:23 AM RE: Simplex Algorithm - Albert Chan - 11-23-2023, 06:35 AM RE: Simplex Algorithm - ftneek - 12-22-2023, 07:38 AM RE: Simplex Algorithm - Albert Chan - 12-22-2023, 04:07 PM RE: Simplex Algorithm - ftneek - 12-22-2023, 07:28 PM RE: Simplex Algorithm - Albert Chan - 12-23-2023, 04:44 AM RE: Simplex Algorithm - ftneek - 12-23-2023, 07:46 AM RE: Simplex Algorithm - Albert Chan - 12-23-2023, 10:23 AM RE: Simplex Algorithm - ftneek - 12-23-2023, 09:30 PM RE: Simplex Algorithm - ftneek - 01-07-2024, 02:40 AM

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