Simplex Algorithm
11-19-2023, 12:57 AM (This post was last modified: 11-20-2023 11:21 PM by Albert Chan.)
Post: #27
 Albert Chan Senior Member Posts: 2,637 Joined: Jul 2018
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
 « 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: 4 Guest(s)