Simplex Algorithm

[Albert Chan]

simplex_le now have maximize default (d=∞), matched simplex_reduce.

simplex_le(a) ≡ simplex_le(a, ∞)

Code:

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)));
 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 := simplex2(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(a) END
  a := k ? pivot((a[j]:=a[j]/a[j][k]),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;

(11-14-2023 11:48 AM)Albert Chan Wrote:  2x + 4y >= 22
3x + 2y >= 20
4x + 5y >= 40
6x + 5y >= cost

Transpose above, we have dual problem.

2 z1 + 3 z2 + 4 z3 <= 6
4 z1 + 2 z2 + 5 z3 <= 5
22 z1 + 20 z2 + 40 z3 <= cost

> simplex_reduce([[2,3,4],[4,2,5]], [6,5], [22,20,40])

\(\frac{320}{7},[0,\frac{10}{7},\frac{3}{7},0,0],\left(\begin{array}{cccccc}
\frac{-6}{7} & 1 & 0 & \frac{5}{7} & \frac{-4}{7} & \frac{10}{7} \\
\frac{8}{7} & 0 & 1 & \frac{-2}{7} & \frac{3}{7} & \frac{3}{7} \\
\frac{46}{7} & 0 & 0 & \frac{20}{7} & \frac{40}{7} & \frac{320}{7}
\end{array}\right) \)

Last row for original problem solution: x = 20/7, y = 40/7, min(cost) = 320/7

> simplex_le([[2,3,4,6],[4,2,5,5],[22,20,40,0]])

\([\frac{320}{7},[0,\frac{10}{7},\frac{3}{7},0,0],\left(\begin{array}{cccccc}
\frac{-6}{7} & 1 & 0 & \frac{5}{7} & \frac{-4}{7} & \frac{10}{7} \\
\frac{8}{7} & 0 & 1 & \frac{-2}{7} & \frac{3}{7} & \frac{3}{7} \\
\frac{46}{7} & 0 & 0 & \frac{20}{7} & \frac{40}{7} & \frac{320}{7}
\end{array}\right) ]\)

[Albert Chan]

ftneek Wrote:[minRatio] algorithm uses this pivoting rule.

Quote:Step 4: Determine the variables xin to ENTER and xout to EXIT BVs.
-> Select Most negative cj <0, corresponding column variable -> xin to ENTER.
->Compute Ratios of (RHS entries bj>=0)/(Entries of that aij>0) -> Lowest Ratio -> Corresponding Row BV: xout to EXIT.
-> Go to step 5 to continue.

Selection is too restrictive. Especially since aij, bj can be anything.

It is possible code never get a chance of making ratios --> POS return 0
In HP Prime, A[0] = A last element, outside valid pivot selection area.

This is the reason canonicalForm() may get into trouble of divide by 0.
Yes, it patched with check aij≠0, but still just a patch.
With a correct minRatio, this should never happen!

Why? Because of added slack/art variables, equation rows are independent.
It is impossible to get variables row of all zeros, even after row operations.
And, slack/art variables are initially setup with 1's, a positive number.
If some other variables are negative, any pivoting will turn them positive.
In other words, positive coefficient(s) are always available, for the picking.

So, what should be the right algorithm for minRatio?
I think maxRatio algorithm is perfect for this! (as in min(x) = -max(-x))
It really is doing maximum ratio, does not care for sign of ratios!
Note: maxRatio is a piece of code inside complex2, not itself a function.

Bonus: for valid setup, canonicalForm() will never get a divide-by-0.
Conversely, divide-by-0 implied matrix was setup wrong, and should be stopped.

Note: I just posted latest complex_le code, thus not added here.

Code:

#cas
simplex(a):=
BEGIN
LOCAL b0,d,v,ign,l,column,art,n,l1,l2,I,J;
b0:={};
d:=dim(a);
v:=d(2)-1;
art:=0;
ign:=0;
n:=1;
FOR I FROM 1 TO d(1)-1 DO
 l:=row(a,I);
 FOR J FROM v downto 1 DO
  IF len(b0)-d(1)+n==0 THEN
   BREAK 2;
  END;
  column:=col(a,J);
  l1:=abs(trn(column));
  IF NOT contains(b0,J) AND l(J)-1==0 THEN
// This should detect artificial variables |last row|= sqrt(numColumns-1 for b,-2 for x1 x2)
   IF column(d(1))-1==0 AND l1-sqrt(2)==0 AND abs(trn(col(a,d(1))))-sqrt(v-2)<=0 THEN
    b0:=CONCAT(b0,{J});
    art:=art+1;
    n:=2;
    BREAK;
   ELSE
// If theres only 1 or 2 1's in a column, we can use it as initial bv.
    IF l1-1==0 OR l1-sqrt(2)==0 THEN
     b0:=CONCAT(b0,{J});
     BREAK;
    END;
   END;
  END;
 END;
END;
RETURN simplex2(a,b0,v,art,ign);
END;

simplex2(a,b0,v,art,ign):=
BEGIN
LOCAL aij,b,bv,c,cj,cmin,bmin,d,n,m,m1,xin,xout,z,l1,l2,I;
bv:=b0;
d:=dim(a);
m:=canonicalForm(a,bv);
b:=mat2list(col(m,d(2)));
c:=SUPPRESS(row(m,d(1)),d(2)-ign,d(2));
cmin:=MIN(c);
n:=1+(art>0);
b:=SUPPRESS(b,d(1)-n+1,d(1));
bmin:=MIN(b);

// Step 2: If all ci are =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(m,I))<=0 THEN
// UNBOUNDED SOLUTION
   RETURN [-inf,[],m];
  END;
 END;

// Step 4: Determine the bv to enter and exit
 xin:=POS(c,cmin);
 b:=mat2list(col(m,d(2)));
 xout:=minRatio(col(m,xin),b,n);

// Step 5: Canonical reduction with new bv
 IF bv[xout]!=xin THEN
  bv[xout]:=xin;
  m:=canonicalForm(m,bv);
  c:=SUPPRESS(row(m,d(1)),d(2)-ign,d(2));
  IF (cmin:=MIN(c))<0 THEN continue END;
 END;

 IF art>0 THEN
  IF m(d(1),d(2))!=0 THEN
// w?0? NO SOLUTION / EMPTY FEASIBLE REGION
   RETURN ["No Solution: empty feasible region (w?0)",[],m];
  END;
  m:=delrows(m,d(1));
  n:=1;
  ign:=art;
  art:=0;
  d=dim(m);
  c:=SUPPRESS(row(m,d(1)),d(2)-ign,d(2));
  cmin:=MIN(c);
 ELSE
  b:=mat2list(col(m,d(2)));
  b:=SUPPRESS(b,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(m,I),d(2)-ign,d(2)))>=0 THEN
     RETURN ["No Solution",[],m];
    END;
   END;
   xout:=POS(b,bmin);    
   bv[xout]:=minRatio(-m[xout],c,0);   /* max ratio */
   m:=canonicalForm(m,bv);
   c:=SUPPRESS(row(m,d(1)),d(2)-ign,d(2));
   cmin:=MIN(c);
  END;
 END;
END;

// Step 6: ci,bi=0 optimum achieved
l1:=sol(m,v);

m1:=m;
z:=-m(d(1),d(2));
c:=SUPPRESS(row(m1,d(1)),d(2)-ign,d(2));
FOR I FROM 1 TO len(c) DO
 IF c(I)==0 AND NOT contains(bv,I) AND abs(trn(col(m1,I)))!=1 AND abs(trn(col(m1,I)))!=0 THEN
  xin:=I;
  b:=mat2list(col(m1,d(2)));
  xout:=minRatio(col(m1,xin),b,1);
  bv[xout]:=xin;
  m1:=canonicalForm(m1,bv);
  l2:=sol(m1,v);
// IFINITE SOLUTIONS
  IF l1==l2 THEN BREAK END
  RETURN [z,l2-(l2-l1)*t,m];
 END;
END;

// UNIQUE SOLUTION
RETURN [z,l1,m];
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;

canonicalForm(M, BV):=
BEGIN
 LOCAL I, J, P;
 FOR I FROM 1 TO len(BV) STEP 1 DO
  IF (P:=M[I][J:=BV[I]])==0 THEN error(M) END
  M:=pivot((M[I]:=M[I]/P),I,J)
 END;
 RETURN M;
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;

Update:
1. skip some canonicalForm() calculations, if we know it will just return itself.
2. unbound solution, return [-inf, [], m], to signal no minimum.

[ftneek]

The changes you propose to canonicalForm causes test_simplex() to fail on my end. Looks like case 3 fails.
> simplex_test()
["[[-2,-4,1,0,0,-22],[-3,-2,0,1,0,-20],[-4,-5,0,0,1,-40],[6,5,0,0,0,0]] Error: Bad Argument Value"]

Check by reverting to:

Code:

canonicalForm(M, BV):=
BEGIN
LOCAL aij;
FOR I FROM 1 TO SIZE(BV) DO
 aij:=M(I,BV(I));
 IF aij !=0 THEN
  M:=mRow(aij^-1,M,I);
  M:=pivot(M,I,BV(I));
 END;
END;
RETURN M;
END;

> test_simplex()
{1,1,1,1,1,1}
I am also using sum(abs(expected - result)) == 0 as you suggested.

[Albert Chan]

Hi, ftneek

Have you updated all changes? I had run your supplied tests before posting, all passed.

BTW, your test_simplex() had a hidden character.

> test_simplex()      → {1,1,1,1,1,test​_case_6}
> test_case_6()       → 1

Copy/paste 1st answer, I get {1,1,1,1,1,test?_case_6}

[ftneek]

I suppose the issue was because I only tested the change to canonicalForm at the time, which I think also needed your change to use max ratio pivoting rule to work. Applying all changes your changes for sol, canonicalForm, and minRatio and extra skips to canonicalForm that I saw seem to work. The tests all passed and time(test_simplex()) on my virtual calculator seems to be ~.005s faster.

I will upload a new source file with the changes. For anyone wondering, simplex2() has been renamed to simplex_core(), but the arguments remain the same. The reason is it may be the starting point or core of other simplex algorithms.

[ftneek]

I uploaded a new source file. I added simplex_int() command, for integer programming problems. It should work for many examples with standard form input, but let me know if you encounter any that don't work. Also for simplicity I had to modify simplex_core() to return a 4th item, the list of basic variables. Not sure if there was a better option.

Example:

min x1-x2
s.t 2x1+x2<=5
-4x1+4x2<=5
x1,x2>=0 and integral (integer)

> simplex_int([[2,1,1,0,5],[-4,4,0,1,5],[1,-2,0,0,0]],{3,4},4,0,0,{1,2})

[-3,[1,2,1,1,0],[[1,0,0,0,1,-1/3,1],[0,1,0,0,1,0,2],[0,0,1,0,-3,2/3,1],[0,0,0,1,0,-4/3,1],[0,0,0,0,1,1/3,3]],{1,2,3,4}]

Thus min=-3 at x1=1,x2=2

Here's one example that currently doesn't work. I think because some non cas commands are used so decimals might be introduced.

Xcas > lpsolve(x1+x2,[1867*x1+1913*x2=3618894],assume=nonnegint,lp_verbose=true)

> simplex_int([[1867,1913,1,3618894],[1,1,0,0],[0,0,1,0]],{3},3,1,0,{1,2})

["No Solution",[],[[0,1,1,3618893/1913,-1410/1913],[1,0,-1913/1867,1/1867,1411/1867],[0,0,46/1867,−6756475144/3571571,-66773/3571571]]]

Here's another example that doesn't work, maybe for the same reason.

> simplex_int([[2,9,1,0,40],[11,-8,0,1,82],[-3,-13,0,0,0]],{3,4},4,0,0,{1,2})
"Error: Bad Argument Value"

[Albert Chan]

(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

[ftneek]

I agree, with only 1 pivot changing each iteration we can save a lot of operations. But I think in the beginning it still needs to go over each pivot at least once, and then the rest of the time it can pivot as you suggest. The reason is because simplex_core accepts a system in canonical or standard form. If it is standard form, I think it needs the pivots. For example for ">=" constraints you can just add negative slack variables if you pivot on them and use dual simplex, which is equivalent to multiplying the constraint by -1 and adding slack variable. Artificial variables also require a pivot, as well as choosing non slack/artificial as initial bv.

I still need to look over the other changes to see how it could have an effect. I will just say I probably prefer artificial variables to be set up outside of simplex_core. It can't be done inside of simplex? Or if you make it another function, it would probably be useful in a program for converting from symbolic input to standard form.

[Albert Chan]

(11-19-2023 07:05 AM)ftneek Wrote:  I agree, with only 1 pivot changing each iteration we can save a lot of operations. But I think in the beginning it still needs to go over each pivot at least once, and then the rest of the time it can pivot as you suggest.

I am not convinced that -1 slack variable will cause problem.
This was previous code, pivot only the art variables (assumed always 1's)

Code:

IF art THEN
 a:=extend(a,matrix(1,v+1));
 FOR J FROM v-art+1 TO v DO
  a[0] -= a[POS(b0,J)];
  a[0][J] := 0;
 END;
END;

But, on the safe side, it now go over all pivots.

Code:

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

Quote:I will just say I probably prefer artificial variables to be set up outside of simplex_core.
It can't be done inside of simplex?

If you prefer art row added from simplex() side, just move "IF art THEN ... END" there.
I personally like core to take care of this details. It allow the same matrix to be used for both.

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

If user doesn't like simplex picked pivots, just use the core with custom pivots.

[ftneek]

For now let's be safe on the first pivots. Your changes had a noticeable speed improvement. Hard to tell when adding more tests, but timing only the first 6 seemed to be a lower number than what I remember. I uploaded a new source file with your changes, but I've moved the IF art THEN inside of simplex(). The tests errored when it was inside simplex_core. simplex_core assumes it doesn't need to set anything up. I think we should keep it that way, at least for now.

[ftneek]

I just posted a new file with a fix for the simplex_int command.

> simplex_int([[2,9,1,0,40],[11,-8,0,1,82],[-3,-13,0,0,0]],[3,4],4,0,0,[1,2])
[-58,[2,4,0,92,90],...]
max=58

Would you know why this returns no solution, but Xcas is able to solve this problem? Is it just not solvable by cutting plane algorithm?
> simplex_int([[1867,1913,1,3618894],[1,1,0,0],[0,0,1,0]],[3],3,1,0,[1,2])
["No Solution",[],...]

Xcas > lpsolve(x1+x2,[1867*x1+1913*x2=3618894],assume=nonnegint,lp_verbose=true)
-> 1916, x1=1009,x2=907

I don't think the issue is having only one constraint. For example test 10 has 1 constraint and was able to be solved by simplex_int(), but I haven't figured out the right command to get xcas to solve that problem yet.

[Albert Chan]

(11-20-2023 12:18 AM)ftneek Wrote:  I just posted a new file with a fix for the simplex_int command.

> simplex_int([[2,9,1,0,40],[11,-8,0,1,82],[-3,-13,0,0,0]],[3,4],4,0,0,[1,2])
[-58,[2,4,0,92,90],...]
max=58

You beat me to it!
I also have simplex_int() bug fixed (culprit is INSERT command), and made code more compact.

Code:

// first attempt of Gomory's plane cutting algorithm
simplex_int(a,b0,v,art,ign,integers):=
BEGIN
 LOCAL b,d,x,m,I;
// Step 1: ignore integer restriction, apply simplex
 a:=simplex_core(a,b0,v,art,ign);
 ign:=art;
 art:=0;
 
 while len(x:=a[2]) DO
 
// Find non integer solution first index
  FOR I FROM 1 to len(integers) DO
   IF type(x[integers[I]])!=DOM_INT THEN break END
  END;

// STOP IF desired variables are integer
  IF I > len(integers) THEN return a END;
  
  a,b0 := a[3..4];
  d:=dim(a);
  b:=col(a,d[2]);
  
// Step 2: Detetmine constraint with largest fractional part
  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));
  b0:=append(b0,d(2)-ign);
  
// Step 3: solve with new constraint
  a:=simplex_core(a,b0,v+1,art,ign);
 END;
 RETURN a;
END;

[ftneek]

There seems to be some issue with the code you suggested, at least on my end.

> test_simplex()
["Unable to eval test in loop : [0.333333333333,0.0,4.66666666667,0.0] Error: Bad Argument Value"]

[Albert Chan]

(11-20-2023 09:02 AM)ftneek Wrote:  > test_simplex()
["Unable to eval test in loop : [0.333333333333,0.0,4.66666666667,0.0] Error: Bad Argument Value"]

HP Prime exact randomly turned float is an old bug, not fixed yet.
It may be due to memory corruption, you may want to reset calculator.

The error message is coming from test_case_7() siimplex_int() first simplex_core() call.

> simplex_core([[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)

[2/3,[1/3,0,14/3,0],[[1,-1/3,0,1/3,1/3],[0,-2/3,1,-1/3,14/3],[0,17/3,0,-2/3,-2/3]],{1,3}]

[Albert Chan]

Question about simplex_int(a,bv,v,art,ign,integers)

Code:

// Step 3: solve with new constraint
  a:=simplex_core(a,bv,v+1,0,ign);

From documentation, v is the number of variables (true + slack + artificial)
Since we keep addding slack variables to a, should (v+1) be more?

However, v is only used for simplex_core sol(a,v), so I guess doc is wrong.
It may meant to solve v leftmost variables.

And, why the +1?

What about simplex_int last argument, what is integers? How to setup problem with it?

---

I thought variable v is redundant = matrix columns - 1
Now that I know better, I am replacing art_row(v,art) as artrow(col,art)
If v=col-1, both matched. If not, using columns is safer.
Below will be on my next update.

>artrow(col,art):=append(makelist(k->k>col-art,2,col),0);
>artrow(9,2)      → [0,0,0,0,0,0,1,1,0]

[Albert Chan]

(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

[ftneek]

(11-20-2023 03:34 PM)Albert Chan Wrote:  Question about simplex_int(a,bv,v,art,ign,integers)

Code:

// Step 3: solve with new constraint
  a:=simplex_core(a,bv,v+1,0,ign);

From documentation, v is the number of variables (true + slack + artificial)
Since we keep addding slack variables to a, should (v+1) be more?

However, v is only used for simplex_core sol(a,v), so I guess doc is wrong.
It may meant to solve v leftmost variables.

And, why the +1?

What about simplex_int last argument, what is integers? How to setup problem with it?

---

I thought variable v is redundant = matrix columns - 1
Now that I know better, I am replacing art_row(v,art) as artrow(col,art)
If v=col-1, both matched. If not, using columns is safer.
Below will be on my next update.

>artrow(col,art):=append(makelist(k->k>col-art,2,col),0);
>artrow(9,2)      → [0,0,0,0,0,0,1,1,0]

Yes my mistake. Technically 'v' should increase every iteration, not just once. But the user may not care about these extra variable solutions, so for now it's probably not a big deal, but can be fixed. I think I meant to do something like ++v. But I think how the 'v' argument is handled may change in the future, probably when I figure out how to implement "unrestricted" variables (and/or symbolic input). simplex_core works with unrestricted variables, but again for now you have to put it in standard form yourself.

Ex: x1 is unrestricted (<=0 or >=0) and integer
-> x1=x1'-x1'' where x1',x1''>=0 and integer

From user's perspective, the original system only has v variables. But once in standard form we've added columns for slack and artificial variables (and fixed unrestricted variables), so from simplex_core perspective we should have v=col-1. As you say, the only difference v has right now is the number of left most solutions to read. If we had unrestricted variables, the solution user expects != solve original v leftmost variables, we have to recover what the unrestricted variables are from the xi=xi'-xi'' relation.

I gave an example a few comments ago, integers is a list of variables with integer restriction. For example
x1,x2 integer -> integers:=[1,2]
x2,x3 integer -> integers:=[2,3]
x1,x3 integer -> integers:=[1,3]
x3 integer -> integers:=[3]

Last the point I wanted to make was that for >= constraint we do not actually need to add an artificial variable at all since we use dual simplex, we either multiply constraint by -1 and add slack or add negative slack, which is the same thing because you still have to multiply the row by -1 to make the column 'canonical'. Thus for ">=" constraint giving input of negative slack for initial bv is fine if we pivot on it or multiply the row by -1 as the 'pivot' operation.

I will look at the rest of the changes a bit later. My main concern is that it could limit the number of ways the user can input a system and get to the same solution. simplex_core is flexible in how the constraints are ordered.

[ftneek]

(11-20-2023 06:12 PM)Albert Chan Wrote:  

(11-19-2023 04:58 PM)Albert Chan Wrote:  I am not convinced that -1 slack variable will cause problem.

I am now convniced that it will not cause problem, even if scaled.
.
.
.
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)

I uploaded a new file with the infinite solution correction you mentioned. But I wasn't able to include this portion of code in simplex_core yet.

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;

I added a new test 4. It is an example of input that is currently allowed by simplex_core, but the command didn't work when I tried your version, because the -1's are not being pivoted on. This what I mean by limiting the number of ways you can input a system and arrive at the same solution.

Current simplex_core
> 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)
[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}]

Your version
> 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)
[0,[0,0,0,0,0],[[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}]

Also, appending an artificial row should not happen inside simplex_core, at least the current way you propose. Reminder, simplex_core is able to continue iteration on a tableau from any stage. For sensitivity analysis, we may consider what happens to the solution if we add a constraint. But just adding a row/bv to the original matrix and applying simplex algorithm again can be expensive for large systems. Instead, we insert the row into the final tableau and continue computation from there. This an example of where your current implementation is wrong, if you are going to set up the artificial row, it should be offset by ign.

[Albert Chan]

(11-21-2023 08:36 AM)ftneek Wrote:  Reminder, simplex_core is able to continue iteration on a tableau from any stage.

I don't know this requirement. Easy fix. '+' = added, '-' = subtracted

simplex() Wrote:IF c > art-(n-=1) THEN continue END
+ IF art THEN
+ a[c+2]:=artrow(v+1,art)-sum(a[POS(bv,J)],J,v+1-art,v);
+ END
RETURN simplex_core(a,bv,v,art,0)

simplex_core() Wrote:- 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;

For efficiency sake, simplex_core assumed pivots are all 1's
This is why my version 4 get 20% speedup running test_simplex().

We can add another wrapper to turn pivots to 1, say, call this simplex2

Code:

simplex2(a,bv,v,art,ign):=
BEGIN
 LOCAL I;
 FOR I FROM 1 to len(bv) DO a:=pivot1(a,I,bv[I]) END
 RETURN simplex_core(a,bv,v,art,ign)
END;

ftneek Wrote:Current simplex_core
> 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)
[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}]

> simplex2([[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)
[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}]

These make simplex_core() restartable at any stage. And, run very fast.

BTW, is above tableau in standard form? Should all constraints be negated?

[Albert Chan]

Minor cosmetic issue.

Instead of returning [min(c), vars, matrix, bv], it might be better do [min(c), matrix, bv, vars]

The reason is the limited display of calculator.
Yes, we can press SHOW, and scroll for answer, but it is too much trouble.

The default , almost always, user see [min(c) ... bv]      // literally, "..."
With vars moved to the end, user see [min(c) ... vars]

What do you think?