Post Reply 
Simplex Algorithm
11-19-2023, 12:57 AM (This post was last modified: 11-20-2023 11:21 PM by Albert Chan.)
Post: #27
RE: Simplex Algorithm
(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
Find all posts by this user
Quote this message in a reply
Post Reply 


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: 4 Guest(s)