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
|
Search
Member List
Calendar
Help
