(11-11-2023 04:54 AM)ftneek Wrote: Note I think there is a bug with my wrapper simplex() method, sometimes the artificial variable counter is not updated. So far the bad effect is for some systems with artificial variables. But also note that simplex2() works as intended with the correct arguments. I need to determine a correct (and hopefully efficient) way to identify artificial variables without false positives, as that would have a bad effect for systems without artificial variables.
It is confusing that artificial variable row is also a valid cost row.
Ask simplex() wrapper to figure out from matrix alone is just asking for trouble.
From what I read from simplex_core(), art variables must be next to col b, in a bunch.
When we are done with art variables, ign=art; art=0, and it work with smaller matrix.
We don't really need art variable row, only number of art variables, which already provided.
Another improvement is a major speedup for problem with lots of constraints.
M := canonicalForm(M, BV) update M by processing a whole list of pivots.
But, each iteration only adjust 1 pivot. Instead of len(BV) pivots to process, all we need is one.
Assuming pivots were good to begin with (col of 0's, only one 1's), no error checking is needed.
It is just a 1 liner: pivot1(m,j,k) := pivot((m[j]:=m[j]/m[j][k]), j,k);
An Introduction to Linear Programming and Game Theory 3rd Edition, example 3.6.1
Minimize, with 2 '=' constraint:
> a:=[[1,-2,-3,-2,3],[1,-1,2,1,11],[2,-3,1,1,0]]
> simplex_le(a, -2)
\([14,[19,8,0,0],\left(\begin{array}{ccccc}
0 & 1 & 5 & 3 & 8 \\
1 & 0 & 7 & 4 & 19 \\
0 & 0 & 2 & 2 & -14
\end{array}\right) ]\)
For conveniece, simplex_le(m, 0) just return m with identity matrix squeezed inside.
> b := simplex_le(a,
0)
\(\left(\begin{array}{ccccccc}
1 & -2 & -3 & -2 & 1 & 0 & 3 \\
1 & -1 & 2 & 1 & 0 & 1 & 11 \\
2 & -3 & 1 & 1 & 0 & 0 & 0
\end{array}\right) \)
We use the identity matrix as artificial variables.
No art row needed. Less typing, less errors!
> simplex(b,
2)
\([14,[19,8,0,0,0,0],\left(\begin{array}{ccccccc}
1 & 0 & 7 & 4 & -1 & 2 & 19 \\
0 & 1 & 5 & 3 & -1 & 1 & 8 \\
0 & 0 & 2 & 2 & -1 & -1 & -14
\end{array}\right) ]\)
simplex(b,2) wrapper called simplex_core(b, [5,6], 6,2,0)
simplex(b) ≡ simplex(b, 0)
Without art variables, it acted the same as previous versions.
> m := [[-2,-4,1,0,0,-22],[-3,-2,0,1,0,-20],[-4,-5,0,0,1,-40],[6,5,0,0,0,0]]
> simplex(m)
[320/7, [20/7,40/7,46/7,0,0], ...]
Code:
#cas
simplex(a):=
BEGIN
LOCAL c,v,colJ,I,J,b0,n,(art:=0);
IF dim(dim(a))!=2 THEN a,art:=a END
n:=art;
c,v:=dim(a).-1; /* constraints, variables */
b0:=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 b0[I] or l1norm(colJ[I]:=0) THEN continue END
b0[I]:=J;
IF c > art-(n-=1) THEN continue END
RETURN simplex_core(a,b0,v,art,0)
END;
error(2,a) /* too few pivots */
END;
simplex_core(a,b0,v,art,ign):=
BEGIN
LOCAL a2,b,c,d,bmin,cmin,n,l1,l2,I,J;
IF art THEN a:=append(a,art_row(v,art)) END
FOR I FROM 1 to len(b0) DO a:=pivot1(a,I,b0[I]) 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 J FROM 1 TO (d(2)-1-ign) DO
IF c(J) < 0 AND MAX(col(a,J))<=0 THEN
// UNBOUNDED SOLUTION
RETURN [-inf,[],a];
END;
END;
// Step 4: Determine b0 to enter and exit
J:=POS(c,cmin);
I:=minRatio(col(a,J),col(a,0),n);
// Step 5: Canonical reduction with new b0
IF b0[I]!=J THEN
b0[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 continue END
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];
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));
bmin:=MIN(b);
// IF ∃bi<0, use Dual Simplex
IF bmin<0 THEN
FOR I FROM 1 TO len(b) DO
IF b(I)<0 AND MIN(SUPPRESS(row(a,I),d(2)-ign,d(2)))>=0 THEN
RETURN ["No Solution",[],a];
END;
END;
I:=POS(b,bmin);
J:=minRatio(-a[I],c,0); /* max ratio */
b0[I]:=J;
a:=pivot1(a,I,J);
c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));
cmin:=MIN(c);
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(b0,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];
END;
// UNIQUE SOLUTION
RETURN [-a[0][0],l1,a];
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;
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 j,k, c := len(a[1]);
LOCAL p := makelist(n:=abs(n)); /* max pivots */
FOR j FROM n DOWNTO 1 DO
p[j] := k := pivot_index(suppress(a[j],c));
IF k==0 and a[j][c] THEN error(3,a) END
a := k ? pivot1(a,j,k) : delrows(a,j);
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);
art_row(v,art) := append(makelist(k->k>v-art,1,v),0);
#end
Update: to detect matrix is has options, I had used dim(dim(a))<2
This assumed dim(plain number)=1, which might be change in the future.
New test for detecting options:
dim(dim(a)!=2
Search
Member List
Calendar
Help
