Simplex Algorithm

[ftneek]

Attached is an implementation of the simplex algorithm as developed in
An Introduction to Linear Programming and Game Theory 3rd Edition, Paul R. Thie and Gerard E. Keough (ISBN: 978-0470232866).

Thanks to Albert Chan for his contributions. I have also included their simplex_le() method.

I originally posted it in this thread, but there's been a few changes since then.

2023/11/11: initial upload
2023/11/17: some speed improvements. simplex2 was renamed to simplex_core(), shown examples will still work if you just change the name.
2023/11/18: added simplex_int() for integer programming problems. simplex_core() returns a list of basic variables as 4th item when there is a solution. Updated and added some integer tests. See post 26 for an example.
2023/11/19: speed improvements, new simplex() wrapper. See post 27 for example inputs.
2023/11/19: fixed exact-> floating point bug for integer problems, more tests
2023/11/21: fixed bug in infinite solution code, new test
2023/12/21:

• new files for simplex, tests, and Game Theory
• improved simplex() command can handle maximization problems, integer restrictions, binary restrictions, or unrestricted variables. No need to manually input slack/artificial variables or negate the objective function row.
• adds solveGame() to solve 2-person zero-sum games
• minor bug fixes, changes to documentation, returns x as a column
• simplex_core() and simplex_int() accept linear programs in canonical form; simplex_core2() and simplex_int2() accept programs in standard form (simplex2 was renamed to simplex_core2)
  

2023/12/22: new files for simplex and tests, bug fix
2023/12/23:

• new files for simplex and tests
• simplex_core returns ∅ for infeasible model with no solution/empty feasible region
• simplex_core and simplex_int accept an additional parameter
• incorporates changes to simplex_le
  

2024/01/06:

• new simplex, tests, and hashin files
• switches to Bland's pivoting rule when a repeated basis is detected to prevent cycling
• infeasible models return "∅" instead of ∅
  

2024/11/17:

• hashin file has been replaced with basisID, small change reflected in simplex file
• files redistributed via zip
  

[Albert Chan]

Hi, ftneek

simplex_le code had a bug, using some HOME functions instead of CAS.
This sometimes turn exact numbers to approximate.

Also, I added "=" constraints, by adding "≥" on top of "≤" constraints.
(≤ constraint) + (≥ constraint) = (= constraint)

It is not efficient, but it is a 1 liner, worthy of adding in, for convenience.

Here is my latest update

Code:

simplex_le(a, dir) :=
BEGIN
 LOCAL b := abs(dir);
 IF b<len(a) THEN a := extend(-a[1 .. b],a) END
 b := dim(a);
 a := transpose(extend(
    [col(a,1 .. b(2)-1)],
    [col(identity(b(1)),1 .. b(1)-1)],
    [col(a,b(2))]
 ));
 IF dir>0 THEN a[0] := -a[0] END /* maximize */
 a := simplex(a);
 IF dir>0 THEN a[1] := -a[1] END /* undo flip */
 return a;
END;

Example: minimize Z = -2x - 3y - 4z
Constraints (3x+2y+z=10) AND (2x+5y+3z=15), x,y,z ≥ 0,

> m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> simplex_le(m, -2)            /* '-' = minimize cost, top 2 are "=" constraints */

[-130/7, [15/7, 0, 25/7, 0, 0, 0, 0], ...]

We can't do ±0. For "≤" constraints only, dir = ±∞, (+ = max, - = min)

HP Prime Help screen simplex_reduce example. (it always maximize)

> simplex_reduce([[3,2,2],[1,1,1]], [3,4], [1,2,3])

[9/2, [0,0,3/2,0,5/2], [[3/2,1,1,1/2,0,3/2],[-1/2,0,0,-1/2,1,5/2],[7/2,1,0,3/2,0,9/2]]]

> simplex_le([[3,2,2,3], [1,1,1,4], [1,2,3,0]], +∞)

[9/2, [0,0,3/2,0,5/2], [[3/2,1,1,1/2,0,3/2],[-1/2,0,0,-1/2,1,5/2],[7/2,1,0,3/2,0,9/2]]]

[ftneek]

Thanks, I modified the original post with an updated hpprgm file.

[Albert Chan]

(11-11-2023 11:51 PM)Albert Chan Wrote:  Also, I added "=" constraints, by adding "≥" on top of "≤" constraints.
(≤ constraint) + (≥ constraint) = (= constraint)

It is not efficient, but it is a 1 liner, worthy of adding in, for convenience.
...
Example: minimize Z = -2x - 3y - 4z
Constraints (3x+2y+z=10) AND (2x+5y+3z=15), x,y,z ≥ 0,

> m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> simplex_le(m, -2)            /* '-' = minimize cost, top 2 are "=" constraints */

[-130/7, [15/7, 0, 25/7, 0, 0, 0, 0], ...]

Some benchmarks.
On my old laptop, running HP Prime emulator (2018 10 16)

> time(simplex_le(m, -2))      → time = 0.0195 s

(11-11-2023 04:54 AM)ftneek Wrote:  > simplex2([[3,2,1,1,0,10],[2,5,3,0,1,15],[-2,-3,-4,0,0,0],[0,0,0,1,1,0]],{4,5},5,2,0)
[-130/7,[15/7,0,25/7,0,0],...]

I have no idea how this is setup, but it run faster , time = 0.0105 s
Use simplex wrapper is easier, but cost a bit more, time = 0.0125 s

> simplex([[3,2,1,1,0,10],[2,5,3,0,1,15],[-2,-3,-4,0,0,0],[0,0,0,1,1,0]])

(11-11-2023 01:43 AM)Albert Chan Wrote:  > m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> m[1] /= m[1][1]
> m := pivot(m,1,1)
> m[2] /= m[2][2]
> m := pivot(m,2,2)

\(\left(\begin{array}{cccc}
1 & 0 & \frac{-1}{11} & \frac{20}{11} \\
0 & 1 & \frac{7}{11} & \frac{25}{11} \\
0 & 0 & \frac{-25}{11} & \frac{115}{11}
\end{array}\right) \)

> simplex_le(m, -∞)      /* "=" constraints are gone! */

[-130/7, [15/7, 0, 25/7, 15/7, 0], ...]

Fastest is remove excess variables, turning it into slack variables.

Note that we scale pivot as 1 before clearing other cells.
This does not affect other constraints, not even cost functions.
In other words, row operations are safe.

Time (including pivot operations) = 0.0075 s
Production code need pivots selection, but should not cost much.

[ftneek]

Thank you for the benchmarks. Now that I thought about it a little, since the simplex_le() method assumes the n constraints are <=, we know we need to add n slack variables. Thus if we make a list of the last n column indices (not counting/using b column) to use as initial basic variables, we can call simplex2() directly as we know all the arguments. This should speed up simplex_le() a bit. This is the idea, but I have a little trouble actually making the list at the moment.

Code:

 a := simplex2(a,makelist(I,I,b(2),b(2)+b(1)),b(2)+b(1)-1,0,0);

We have 0 artificial variables and should ignore 0 columns.

A similar idea could be used for artificial variables. Since simplex_le() assumes the top dir constraints are "=", this makes things easy. After the added identity matrix, we should add dir more columns (identity matrix like). But we should also add a row of 0's with 1's in positions of the new dir columns (the last dir columns, not counting/using b column). Again we can make a list of indices of the added columns and call simplex2 directly, bypassing the simplex() wrapper. Except this time we call simplex2(a,list,v,dir,0) since we have dir artificial variables. This is an "initial tableau" so we should not ignore any columns yet.

(11-12-2023 12:38 AM)Albert Chan Wrote:  

(11-11-2023 04:54 AM)ftneek Wrote:  > simplex2([[3,2,1,1,0,10],[2,5,3,0,1,15],[-2,-3,-4,0,0,0],[0,0,0,1,1,0]],{4,5},5,2,0)
[-130/7,[15/7,0,25/7,0,0],...]

I have no idea how this is setup, but it run faster , time = 0.0105 s

It is just manually entering the arguments simplex() tries to find. 2 constraints -> 2 basic variables. since all constraints are "=", they are the positions of 1's in the last row, hence positions {4,5}. We have 3 variables and 2 artificial variables -> v=5. 2 "=" constraints -> 2 artificial variables. Because we have artificial variables and it's an initial tableau, the last argument is 0.

(11-12-2023 12:38 AM)Albert Chan Wrote:  

(11-11-2023 11:51 PM)Albert Chan Wrote:  > m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> simplex_le(m, -2)
[-130/7, [15/7, 0, 25/7, 0, 0, 0, 0], ...]

> time(simplex_le(m, -2))      → time = 0.0195 s

Compare the time when supplying the equivalent inputs for simplex2 using <= + >= method, to see the time saved by bypassing the simplex wrapper.

> simplex2([[-3,-2,-1,1,0,0,0,-10],[-2,-5,-3,0,1,0,0,-15],[3,2,1,0,0,1,0,10],[2,5,3,0,0,0,1,15],[-2,-3,-4,0,0,0,0,0]],{4,5,6,7},7,0,0)
[-130/7,[15/7,0,25/7,0,0,0,0],...]

[Albert Chan]

(11-12-2023 12:38 AM)Albert Chan Wrote:  Fastest is remove excess variables, turning it into slack variables.

Note that we scale pivot as 1 before clearing other cells.
This does not affect other constraints, not even cost functions.
In other words, row operations are safe.

Time (including pivot operations) = 0.0075 s
Production code need pivots selection, but should not cost much.

I don't know a good name for this algorithm, so I just use remove_eq_c(a, n)
It remove top n "=" constraints. (n variables, if constraints are independent)
In other words, only "≤" constraints remain.

simplex_le(a, ±n) ≡ simplex_le(remove_eq_c(a,n), ±∞)

simplex_le(a, dir) only had 1 change. Instead of more constraints, it now do less variables.

< IF b<len(a) THEN a := extend(-a[1 .. b],a) END
> IF b<len(a) THEN a := remove_eq_c(a, b) END

Code:

pivot_index(a) :=
BEGIN    
    LOCAL b := contains((a:=abs(a)), 1.0);
    IF b THEN RETURN b END
    IF (b:=max(a)) THEN b := contains(a,b) END
    RETURN b;
END;

remove_eq_c(a, n) :=
BEGIN
  LOCAL p, r, c := len(a[1]);
  FOR r FROM n DOWNTO 1 DO
    p := pivot_index(suppress(a[r],c));
    IF p==0 and a[r][c] THEN error(a) END
    a := p ? pivot((a[r]/=a[r][p]),r,p) : delrows(a,r);
  END
  return a;
END;

simplex_le(a, dir) :=
BEGIN
 LOCAL b := abs(dir);
 IF b<len(a) THEN a := remove_eq_c(a, b) END
 b := dim(a);
 a := transpose(extend(
    [col(a,1 .. b(2)-1)],
    [col(identity(b(1)),1 .. b(1)-1)],
    [col(a,b(2))]
 ));
 IF dir>0 THEN a[0] := -a[0] END /* maximize */
 a := simplex2(a,range(b(2),sum(b)-1),sum(b)-2,0,0);
 IF dir>0 THEN a[1] := -a[1] END /* undo flip */
 return a;
END;

Benchmark example:

> m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> time(r := simplex_le(m, -2))                     → 0.0077
> r

\([\displaystyle
\frac{-130}{7},[\frac{15}{7},0,\frac{25}{7},\frac{15}{7},0],\left(\begin{array}{cccccc}
1 & \frac{1}{7} & 0 & 1 & \frac{1}{7} & \frac{15}{7} \\
0 & \frac{11}{7} & 1 & 0 & \frac{11}{7} & \frac{25}{7} \\
0 & \frac{25}{7} & 0 & 0 & \frac{25}{7} & \frac{130}{7}
\end{array}\right) ]\)

Update Nov 12, 2023: simplex_le directly call simplex2, instead of wrapper simplex.
see https://www.hpmuseum.org/forum/thread-20...#pid180152

[ftneek]

(11-12-2023 05:11 PM)Albert Chan Wrote:  > m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> time(r := simplex_le(m, -2))                     → 0.0077
> r

\([\displaystyle
\frac{-130}{7},[\frac{15}{7},0,\frac{25}{7},\frac{15}{7},0],... ]\)

(11-11-2023 04:54 AM)ftneek Wrote:  > simplex2([[3,2,1,1,0,10],[2,5,3,0,1,15],[-2,-3,-4,0,0,0],[0,0,0,1,1,0]],{4,5},5,2,0)
[-130/7,[15/7,0,25/7,0,0],...]

Reminder: to avoid incorrect answers for the Dual solution, the 15/7 should really be a 0. This likely just requires either a small modification or new sol() command. The sol command reads a solution for variable I if there there's a 1 and rest 0's in a column. If the column is not in canonical form, the I-th variable is 0.

Specifically, your output has 3 columns in "canonical" form, columns 1,3,4. But column 4 has the 1 in the same position as column 1, thus the answer read for x4 should probably be 0. This issue does not arise when the system is in standard form and the algorithm is run traditionally.

[Albert Chan]

(11-12-2023 05:40 PM)ftneek Wrote:  Reminder: to avoid incorrect answers for the Dual solution, the 15/7 should really be a 0. This likely just requires either a small modification or new sol() command.

Code:

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.0 THEN continue END
 IF (k:=contains(c,1.0))==0 THEN continue END
 IF b[k] THEN x[j]:=b[k]; b[k]:=0 END
END;
RETURN x;
END;

[ftneek]

Ok, your change to sol() don't seem to change the output of simplex2() (on the examples from the thread that I tried), so I think I will accept it.

Also, I misspoke earlier. Before, the 15/7 effected the primal solution (though it does not really matter as it is an artificial variable, they usually end up being 0). What matters for the Dual solution is the values of the final c row.

In short, for this example:
Min -> Max
2 "=" constraints -> y1,y2 unrestricted
x1,x2,x3>=0 -> 3 "<=" constraints
Constant terms of constraints -> Objective function coefficients
Objective function coefficients -> Constraint terms of constraints
Constraint coefficients matrix A^t -> A

Quote:Solving a minimization problem containing equality constraints, a sign change adjustment is necessary when determining an optimal solution point to the dual. The value of each unrestricted variable in the optimal solution point to the dual is the negative of the entry in the bottom row of the associated artificial variable column. (Because for a maximization problem, the corresponding minimization problem is considered, which necessitates an initial sign change in the objective function coefficients when entered in to the initial tableau.)

Thus y1=2/7,y2=-10/7 is the solution to the associated dual maximization problem. Negative solutions are valid, since y2 is unrestricted. This is generally relevant and important information for price analysis and business related questions.

When I get around to it (busy today), I think it would be helpful to compile some test cases of the examples shared on the forum so far. Any proposed changes to the code base should not affect the output of simplex2() with correct arguments.

[Albert Chan]

(11-12-2023 01:44 AM)ftneek Wrote:  Thank you for the benchmarks. Now that I thought about it a little, since the simplex_le() method assumes the n constraints are <=, we know we need to add n slack variables. Thus if we make a list of the last n column indices (not counting/using b column) to use as initial basic variables, we can call simplex2() directly as we know all the arguments. This should speed up simplex_le() a bit.

Nice idea. Updated to simplex_le latest version.

< a := simplex(a);
> a := simplex2(a,range(b(2),sum(b)-1),sum(b)-2,0,0);

range(...) are column location of 1's of identity matrix.
variables = rows + cols - 2 = sum(b) - 2
No artificial variables, thus the doubled 0's

(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.

Call simplex2 directly, simplex_le is more robust. Run faster too (was 0.0077 s)

> m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> time(simplex_le(m, -2))                        → 0.0058

Just curious, if simplex2 can assume no artificial variables, can it speed up more?

[ftneek]

(11-12-2023 11:04 PM)Albert Chan Wrote:  Just curious, if simplex2 can assume no artificial variables, can it speed up more?

simplex2 makes no assumption. Well, it assumes that it was provided a matrix in standard or canonical form and enough basic variables to use. It is able to take a tableau from any stage and continue iteration, but it needs inputs to know where to start from. But I suppose if you know in advance that you have no "=" constraints, then you know you will have 0 artificial variables. I guess there could be a separate version of simplex2, that assumes no artificial variables, and the parts/logic involving artificial variables could be altered or removed from the method. But that would require maintaining 2 versions of essentially the same code.

Artificial variables only really affect the indices to use in checks of b column and c row, and there is an IF art>0 inside the while loop. It is only used once, but the if statement is reached when the algorithm decides to exit the while loop or use Dual Simplex. Once artificial variables have served their purpose, ign:=art (so we will not consider the artificial entries in the c row, they are allowed to be negative) and art:=0, the IF art>0 statement has been used and won't be needed again. Also n:=1, now we only ignore the c row/z entry in the b column, as the artificial row was deleted. I'm not sure the time saved by removing those parts of the code and maintaining a separate version would be worth it. But maybe it could be possible to simplify the arguments of simplex2.

There are potentially optimizations that could be made as well. For example sometimes the matrix is already given in canonical form, or some columns already are. Then it makes sense to omit row operations from the basic variable columns already in canonical form, but I'm not sure if detecting that would have a smaller computation cost than just keeping it as is. The IF aij!=0 in canonicalForm is to prevent division by zero, and may already handle it.

Again, I would like to set up some tests to ensure that any changes to the code do not effect the output of simplex2.

[Albert Chan]

Another optimization, remove_eq_c(a,n) now also return pivots information.
simplex_le(a, ±n) use pivots to make less slack variables (smaller identity matrix)

Code:

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;

remove_eq_c(a, n) :=
BEGIN
  LOCAL j,k, c := len(a[1]);
  LOCAL p := makelist(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;

simplex_le(a, dir) :=
BEGIN
 LOCAL b, p := [];
 LOCAL n := abs(dir);
 IF n<len(a) THEN a,p := remove_eq_c(a,n) 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)));
 n := sum(b)-2 - len(p);         /* variables */
 p := extend(p,range(b[2],n+1)); /* pivots */
 IF dir>0 THEN a[0] := -a[0] END /* maximize */
 a := simplex2(a,p,n,0,0);
 IF dir>0 THEN a[1] := -a[1] END /* undo flip */
 return a;
END;

Benchmark example, time reduced to 0.0050 s. For this example, no slack variables!

> m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> simplex_le(m, -2)

\(\displaystyle [
\frac{-130}{7},[\frac{15}{7},0,\frac{25}{7}],\left(\begin{array}{cccc}
1 & \frac{1}{7} & 0 & \frac{15}{7} \\
0 & \frac{11}{7} & 1 & \frac{25}{7} \\
0 & \frac{25}{7} & 0 & \frac{130}{7}
\end{array}\right) ]\)

An Introduction to Linear Programming and Game Theory 3rd Edition, Paul R. Thie and Gerard E. Keough
Chapter 3.7 Redundant system, example 2. Again, no slack variables!

> m := [[1,2,0,1,20], [2,1,1,0,10], [-1,4,-2,3,40], [1,4,3,2,0]]
> simplex_le(m, -3)

\([35,[5,0,0,15],\left(\begin{array}{ccccc}
0 & \frac{3}{2} & \frac{-1}{2} & 1 & 15 \\
1 & \frac{1}{2} & \frac{1}{2} & 0 & 5 \\
0 & \frac{1}{2} & \frac{7}{2} & 0 & -35
\end{array}\right) ]\)

We have redundant "=" constraint, net of only 2 "=" constraints.

> remove_eq_c(m, 3)

\([\left(\begin{array}{ccccc}
0 & 1 & \frac{-1}{3} & \frac{2}{3} & 10 \\
1 & 0 & \frac{2}{3} & \frac{-1}{3} & 0 \\
0 & 0 & \frac{11}{3} & \frac{-1}{3} & -40
\end{array}\right) ,[2,1] ]\)

This simulate my previous version, not using pivot information.

> simplex_le(Ans[1], -inf)

\([35,[5,0,0,15,0,0],\left(\begin{array}{ccccccc}
0 & \frac{3}{2} & \frac{-1}{2} & 1 & \frac{3}{2} & 0 & 15 \\
1 & \frac{1}{2} & \frac{1}{2} & 0 & \frac{1}{2} & 1 & 5 \\
0 & \frac{1}{2} & \frac{7}{2} & 0 & \frac{1}{2} & 0 & -35
\end{array}\right) ]\)

Update: HP Prime had a weird behavior list /= k interpeted as list := list * (1/k)

> m := [3.0]
> m /= 3.0
> m[1] - 1      → −3.5527136788e−15

Thus, in remove_eq_c, to turn pivot as 1. (exactly), we do a[j] := a[j]/a[j][k]
I have made changes to ftneek code as well. Now pivot of 1. is exactly 1
sol(a,v) (and others) are reverted back to original, comparing pivot 1 instead of fuzzy 1.0

[ftneek]

Albert, I do think it is interesting that it seems your "shortcut" method is able to return the same primal solution in the examples you've shown. If the shortcut method correctly provides solutions to larger systems with mixed constraints, and you are not interested in the dual solution, I suppose it could be worth keeping as an option if speed is truly a concern on larger systems.

But the algorithm was designed to simultaneously solve the dual problem as well. If you do not use artificial variables, "half" the information the algorithm was designed to provide is missing. I guess it makes sense then that your shortcut takes roughly half as much time... But I think for this shortcut to be absolutely useful, there should be some way to recover the final tableau obtained by calling simplex2 on the system in proper standard form.

Additionally, slack variables generally represent “unused” resources. For example if the constraint is x1 + x2<=30, and x1=5 and x2=10, then the slack variable should equal 15, because that is the remaining quantity of whatever the constraint represents (x1+x2+s1=30 -> s1=15). It is nice to provide this information as well with the solutions, eg provide x3=15.

Quote:Since an artificial variable cannot be removed from the basis, the original system of constraints contains one redundant equation.

Also, I think this may be a nice piece of information to include in the result. If an artificial variable remains a basic variable at the final tableau, then the system contains a redundant equation. This should be easy enough to check. But do you have any suggestions about how to present this information to the user? I think it is probably best to keep the return type as a list of 3 items if possible.

Proper standard form arguments for Example 3.7.2:
> simplex2([[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],[0,0,0,0,1,1,1,0]],{5,6,7},7,3,0)

[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]]]

[Albert Chan]

Hi, ftneek

It is not just a shortcut (as in hack), there is no lost of information.
Row operations represent *exactly* the original problem, but with less variables to handle.
We are adding equal amount to both side. All constraints, even cost function is unaffected.

> m := [[1,2,0,1,20], [2,1,1,0,10], [-1,4,-2,3,40], [1,4,3,2,0]]
> m2 := remove_eq_c(m, 3)[1]

\(\left(\begin{array}{ccccc}
0 & 1 & \frac{-1}{3} & \frac{2}{3} & 10 \\
1 & 0 & \frac{2}{3} & \frac{-1}{3} & 0 \\
0 & 0 & \frac{11}{3} & \frac{-1}{3} & -40
\end{array}\right) \)

Once variable are unrelated to cost function, it automatically becomes "slack".
You can adjust them whatever you want, to fit your need (max or min)

Yes, identity matrix on the "wrong" side, but it does not matter.
(if you prefer, swap columns to make identity matrix look right)

We don't need any more slack variables.

> simplex2(m2, [2,1], 4, 0, 0)

\([35,[5,0,0,15],\left(\begin{array}{ccccc}
0 & \frac{3}{2} & \frac{-1}{2} & 1 & 15 \\
1 & \frac{1}{2} & \frac{1}{2} & 0 & 5 \\
0 & \frac{1}{2} & \frac{7}{2} & 0 & -35
\end{array}\right) ]\)

ftneek Wrote:If an artificial variable remains a basic variable at the final tableau, then the system contains a redundant equation. This should be easy enough to check. But do you have any suggestions about how to present this information to the user?

Problem matrix of 4 rows, Solution matrix of 3 rows --> 1 redundant "=" constraint.

[Albert Chan]

(11-13-2023 05:34 PM)ftneek Wrote:  Proper standard form arguments for Example 3.7.2:
> simplex2([[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],[0,0,0,0,1,1,1,0]],{5,6,7},7,3,0)

[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]]]

You can do this with simplex_le too, but not much typing is saved.
Constraints already are "=", artificial variables need to start with -1, thus result negated.

> simplex_le([[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]],-3)

[35,[5,0,0,15,0,0,0],[[0,3/2,-1/2,1,-1/4,0,-1/4,15],[0,0,0,0,-3/2,1,1/2,0],[1,1/2,1/2,0,-3/4,0,1/4,5],[0,1/2,7/2,0,5/4,0,1/4,-35]]]

Inside numbers not quite match, because simplex_le use pivots order {5,6,1}. But that shouldn't matter.

[ftneek]

(11-14-2023 02:25 AM)Albert Chan Wrote:  You can do this with simplex_le too, but not much typing is saved.
Constraints already are "=", artificial variables need to start with -1, thus result negated.

> simplex_le([[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]],-3)

[35,[5,0,0,15,0,0,0],[[0,3/2,-1/2,1,-1/4,0,-1/4,15],[0,0,0,0,-3/2,1,1/2,0],[1,1/2,1/2,0,-3/4,0,1/4,5],[0,1/2,7/2,0,5/4,0,1/4,-35]]]

Inside numbers not quite match, because simplex_le use pivots order {5,6,1}. But that shouldn't matter.

Does that input still save computing time? That is closer, potentially promising. Order of pivots does not matter. But the coefficients of the objective function/z row should match simplex2, including in sign, else the dual solution will not be correct. To me it seems those columns need to be multiplied by -1, except for the entries where the b coefficient is 0... But to me this just seems to be creating new, strange rules and complicating things. For example adding negative slack variables to "=" constraints, why? That does serve a purpose, but I don't think it's explained in this textbook. I know another book that refers to that has "soft and hard constraints" for when you sometimes want to allow the "=" constraints to be violated for a cost, my upcoming project actually uses that idea. If it does not save computing time, I suggest just using the intended standard form.

As I mentioned in the other thread, artificial variables are not always required to solve a system, but can always be used as starting place for the simplex algorithm, and their use cuts down on computation time (when compared to size of <= + >= matrix). Learning how to put systems in standard form is not very difficult, there are only a few rules to follow for when the following arise: ">=","=","<=", unrestricted variables. I think that's everything.. Once you learn the few rules you should be able to handle most systems. But I agree it is not convenient to have to enter this information manually. A program to convert to standard form for the user would be immensely convenient and useful.

That was what I was going to start with originally, but handling user input would have taken too much time for me then. So I just wrote the main simplex2 code (the original simplex()) with the assumption that it will be provided the correct inputs to use, ideally by a helper method that sets everything up for you, which was my goal for simplex() (current version).

The algorithm does make some assumptions, and in return makes some guarantees. Namely if it is provided a non degenerate matrix in standard form, it is guaranteed to return the optimal solution to the primal and dual system (or determine no solution) in a finite number of iterations. Even if the matrix is degenerate, there are a few things we can try that should guarantee finiteness. If the assumption is not met, I don't know what is guaranteed.

If you want an algorithm that does not use artificial variables, there are some out there. I encourage you to read the article/link in the source comment. If I were to implement another simplex algorithm, it would probably be that one unless I find something more promising. It uses no artificial variables, requires less multiplications and additions, and has the dual solution embedded in the matrix. But I follow this text this semester, as my course uses it heavily.

I don't want to make it seem like I'm squashing this idea. If you just want a faster simplex algorithm for large problems, that's fine. But (unless you can figure out how to consistently match the signs of simplex2) it should be made clear that the intention is to only provide solutions to the primal, and that finding the dual solution would require converting to dual problem AND rewriting in standard form (maybe not if using your method) AND another round of simplex iteration, probably expensive for a large system. Needing to run a second tableau, the computation time is probably at least back to running simplex2 as originally intended... My original intention was for extensions of the program to make it easier to use, extract the desired information, or add functionality without limiting its capabilities. In my opinion without the dual solution available, the capabilities are being limited. But it should useful if you want only the primal, fast.

[Albert Chan]

(11-14-2023 06:24 AM)ftneek Wrote:  

(11-14-2023 02:25 AM)Albert Chan Wrote:  You can do this with simplex_le too, but not much typing is saved.
Constraints already are "=", artificial variables need to start with -1, thus result negated.

> simplex_le([[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]],-3)

[35,[5,0,0,15,0,0,0],[[0,3/2,-1/2,1,-1/4,0,-1/4,15],[0,0,0,0,-3/2,1,1/2,0],[1,1/2,1/2,0,-3/4,0,1/4,5],[0,1/2,7/2,0,5/4,0,1/4,-35]]]

Inside numbers not quite match, because simplex_le use pivots order {5,6,1}. But that shouldn't matter.

Does that input still save computing time? That is closer, potentially promising. Order of pivots does not matter. But the coefficients of the objective function/z row should match simplex2, including in sign, else the dual solution will not be correct. To me it seems those columns need to be multiplied by -1, except for the entries where the b coefficient is 0...

Constraint (a*x ≤ b) is setup as matrix [a|b], with a 1 column b
To turn "≤" into "=", we do a*x + (slack var) = b + (artificial var)

But we can't setup this way, thus a*x + (slack var) + (−artficial var) = b

(−artificial var) = [5/4,0,1/4]
(artificial var) = −[5/4,0,1/4]      // answer do match simplex2

Purpose of wrapper function is not to save computing time, but to save typing time.
But yes, to match result of standard form, you might as well just use simplex2.

Again, simplex_le() is just a wrapper for simplex2().
I have nothing against the standard form, except for typing it all in.





simplex_le can do dual problem, with identity matrix squeeze inside.
Video 1st dual probelm (time = 6 .. 17 minutes)

> simplex_le([[2,4,40],[1,1,12],[5,1,40],[20,30,0]],+∞)

\([320,[4,8,0,0,12],\left(\begin{array}{cccccc}
0 & 1 & \frac{1}{2} & -1 & 0 & 8 \\
1 & 0 & \frac{-1}{2} & 2 & 0 & 4 \\
0 & 0 & 2 & -9 & 1 & 12 \\
0 & 0 & 5 & 10 & 0 & 320
\end{array}\right) ]\)

Note that P column is always [0,0 .. 0,1] thus not needed to carry along.

[ftneek]

(11-14-2023 09:23 AM)Albert Chan Wrote:  simplex_le can do dual problem, with identity matrix squeeze inside.
Video 1st dual probelm (time = 6 .. 17 minutes)

> simplex_le([[2,4,40],[1,1,12],[5,1,40],[20,30,0]],+∞)

That is impressive, much less typing, and matches the standard form output of simplex2. Can simplex_le() already solve dual problem with "=" constraints? I think the only thing it really needs is a way to set up artificial variables.

Odd, for whatever reason the video explained it the opposite way from my understanding. According to our source, primal is given by b column while dual is found in c row. Anyhow, the final tableau matches simplex2 exactly!

If you would like to make another optimization, there is potentially one regarding duality. Sometimes the dual or the primal is easier to solve (it might have positive variables or only needs slack variables, or a certain quality that makes it "easier"). Choosing the "easier" problem to solve may require less computation, and the solution to the harder problem is provided by duality.

Finally, I had a little time to start making some "tests". But for some reason test 3 returns 0, even after directly copying the result from the command history... Would you happen to know why? I'm not sure if this is the best way to maintain tests, and it would be nice if they don't all show up in the catalog, only a single wrapper that runs them. Since we only are checking the equivalency of 2 items, does it really need to be a cas program? I would appreciate tips on how to organize it, I think it's important to ensure refactoring simplex2 does not cause correct solutions to change.

Code:

#cas
test_simplex():=
BEGIN
RETURN {test_case_1(),test_case_2(),test_case_3(),test_case_4()};
END;

// Example 3.7.2
test_case_1():=
BEGIN
RETURN simplex2([[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],[0,0,0,0,1,1,1,0]],{5,6,7},7,3,0)==[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]]];
END;

test_case_2():=
BEGIN
RETURN simplex2([[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],[0,0,0,0,0,1,1,1,0]],{6,7,8},8,3,0)==[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]]];
END;

test_case_3():=
BEGIN
RETURN 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]]];
END;

test_case_4():=
BEGIN
RETURN 4;
END;
#end

[Albert Chan]

Hi, ftneek

Why simplex2() able to handle negative b column, but HP Prime default simplex_reduce() can't ?

Is the b≥0 requirement arbitrarily added, or does it serve a purpose?

[ftneek]

simplex2 incorporates Dual Simplex algorithm. Using (non dual) simplex, in the end if c row is all positive we assume we have reached the optimum. But if there is a negative in b column, the answer is considered infeasible as the assumption is all variables are non negative. Using dual simplex, if the "final tableau" would have a negative b iteration continues. If a c becomes negative again, regular simplex is used until all c are positive and there is a negative b.

So I think b>=0 requirement means only regular simplex will be used.

For whatever reason test 3 kept returning 0 on my virtual and physical calculator despite both sides of "==" being visually the same. I believe there to be some kind of issue, at least on my end. Why would "==" return 0 if the objects are the same? It seemed to work fine on the other 2 tests. Even if I were to reset my virtual and physical calcs, I hope to figure out the cause so that this issue doesn't randomly pop up again sometime in the future.

Using size(difference())==0 seems to work.. But I guess using this, it would be possible for a false pass if the values are all the same but ordered differently.

Code:

#cas

test_simplex():=
BEGIN
RETURN {test_case_1(),test_case_2(),test_case_3(),test_case_4(),test_case_5(),test​_case_6()};
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]]];
result:=simplex2([[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],[0,0,0,0,1,1,1,0]],{5,6,7},7,3,0);
RETURN size(difference(result,expected))==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]]];
result:=simplex2([[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],[0,0,0,0,0,1,1,1,0]],{6,7,8},8,3,0);
RETURN size(difference(result,expected))==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]]];
result:=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);
RETURN size(difference(result,expected))==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]]];
result:=simplex2([[3,2,1,1,0,10],[2,5,3,0,1,15],[-2,-3,-4,0,0,0],[0,0,0,1,1,0]],{4,5},5,2,0);
RETURN size(difference(result,expected))==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],[[1,0,1/17,0,0,32/85,-1/17,1/17,-8/85,131/85],[0,1,3/17,0,0,28/85,-3/17,3/17,-7/85,189/85],[0,0,-1/17,1,-1,-3/17,1/17,-1/17,5/17,35/17],[0,0,13/17,0,0,5/17,4/17,-4/17,3/17,-1441/17]]];
result:=simplex2([[4,-1,0,1,-1,1,0,0,0,6],[-7,8,1,0,0,0,-1,1,0,7],[1,1,0,4,-4,0,0,0,1,12],[-3,2,1,-1,1,0,0,0,0,-87],[0,0,0,0,0,0,0,1,1,0]],{6,8,9},9,2,0);
RETURN size(difference(result,expected))==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]]];
result:=simplex2([[2,3,4,1,0,6],[4,2,5,0,1,5],[-22,-20,-40,0,0,0]],{4,5},5,0,0);
RETURN size(difference(result,expected))==0;
END;

#end

edit: I added the rest of the examples shared in the threads as tests and added the problem statements to the comments. For any change to simplex2, test_simplex() should still return all 1's (at minimum). Feel free to add more examples as you encounter them. If you want to check an answer before adding it, lpsolve on xcas can be used as a reference. It might be nice to have a method that can read problems stored in a file, once we have a method to convert a system to standard form. That would probably make it easy to transfer problems found online to the calculator if needed. If anyone has a hint for not exposing all the test cases in the catalog, please let me know.