simplex

# Linear Programming: the Simplex method for the HP-41

Overview

1°)  The Big-M Method
2°)  Artificial-variable free Algorithm

1°)  The Big-M Method

-The purpose is to find  m  non-negative real numbers:   x1 , ..... , xm  satisfying:

bi³ ai;1x1 + ....... + ai;m xm            i = 1 , .... , n                                ( n inequations )                   all the bi  , bi'  and  bi"
bi' = ai';1x1 + ....... +ai';m xm            i' = 1 , ..... , n'                              ( n' equations )                    must be non-negative*
bi"£ ai";1x1+ ....... +ai";m xm          i" = 1 , ..... , n"                             (n" inequations )

and such that      F  =  c1x1 + ....... +  cmxm       is maximized.    ( to find the minimum of F, seek the maximum of  -F )

* if some bi are negative, multiply the corresponding lines by -1.

Data registers:

- R01 thru R(p+1)(m+n"+p+1) are used for the simplex array. ( p = n + n' + n" )
- The coefficients of the system must be stored column by column from R01 to Rm(p+1)   ( store 0 in Rp+1 )
- When the program stops,

R00 = F  maximum
R01 = x1 , .... , Rmm = xm
Rm+1 ...etc... = The slack and/or artificial variables.

Flags:    F01 and F02. If at least one of these flags is set when the program stops, there is no solution ( or no bounded solution ).

Subroutines:  /

A short analysis of the program:

-Lines 01 thru 123 store the other coefficients of the simplex matrix.
If there are only inequations of the first type ( ³ ) it's quite easy: the right half is the unit matrix.
But if there are equations and/or inequations of the second type ( £ ) it's much more complex:
all the last row has to be changed because of the artificial variables.

-The idea is to weight these variables by a coefficient -M where M is such a large positive number
that for a maximum, all the artificial variables will eventually turn out to be zero.
On a calculator, we have to make a choice. In this program:   M = 100   ( lines 41-42 and line 103 ).

-Lines 124 thru 214 are the heart of the simplex method itself.
-Lines 215 thru 287 move the results to registers R00 to Rmm ... and verify:
1- if a solution is found
2- if all the artificial variables are really 0 ( lines 274 to 282 )
-Otherwise F01 and/or F02 are set.

Warning:

-Because of the choice M = 100 , all the cj  must be significantly smaller than 100.
-For instance, if  F = 2400 x + 1200 y  it would be better to find the maximum of  2.4 x + 1.2 y and to multiply the result by 1000.

-This recommendation does not apply if there are only inequations of the first type ( ³ )  i-e if  n' = n" = 0.

-This program finds only 1 solution ( even if there are several solutions ).

Notes:

-Status registers  M N O and a  are used.
-Therefore, "SPLEX" must not be called as more than a first level subroutine.
-Lines 137 and 214 are synthetic three-byte GTOs.
-If you don't have an HP-41CX, replace line 23 by these 6 lines:

0
LBL 14
STO IND Y
ISG Y
GTO 14
RDN                        ( another possibility is to execute CLRG before storing the coefficients of the matrix )

 01  LBL "SPLEX"   02  STO 00   03  RDN   04  STO O   05  +   06  X<>Y   07  STO N   08  +   09  1   10  ST+ 00   11  +   12  STO M   13  +   14  RCL O   15  +   16   E3   17  /   18  RCL 00   19  +   20  RCL M   21  *   22  ISG X   23  CLRGX   24  FRC   25  RCL M   26   E5   27  /   28  +   29  RCL 00   30  RCL M   31  ST+ Y   32  *   33  +   34  STO Y   35   E-5   36  +   37  RCL O   38  X=0?   39  GTO 00   40  -   41  2   42  10^X   43  1   44  LBL 01   45  ST- IND Z   46  X<>Y   47  ST- IND T   48  ISG Z 49  X<>Y   50  ISG T   51  GTO 01   52  LBL 00   53  RCL 00   54  .1   55  %   56  +   57  RCL M   58  STO Z   59  DSE X   60  +   61  *   62  DSE X   63  RCL M   64  1   65  +   66   E5   67  /   68  +   69  SIGN   70  LBL 02   71  STO IND L    72  DSE L   73  GTO 02   74  RCL M   75  RCL 00   76  DSE X   77  ST+ Y   78  RCL N   79  +   80   E3   81  /   82  +   83  STO a   84  RCL N   85  RCL M   86  DSE X   87  X=Y?   88  GTO 00   89   E-3   90  STO Z   91  *   92  RCL 00   93  X<> N   94  +   95  ISG X   96  ST+ Y 97  0   98  LBL 03   99  RCL IND Y  100  + 101  ISG Y 102  GTO 03 103   E2 104  * 105  ST+ IND Y 106  RDN 107  + 108  0 109  DSE N 110  GTO 03 111  LBL 00 112  RCL 00 113  RCL O 114  + 115  DSE X 116  RCL M 117  ST+ Y 118  ST* Y 119   E2 120  / 121  + 122   E3 123  / 124  ISG X 125  STO 00 126  CF 01 127  CF 02 128  LBL 04 129  RCL 00 130  STO O 131  LBL 05 132  ISG O 133  GTO 05 134  RCL 00 135  ISG X 136  STO M 137  GTO 10 138  LBL 05 139  RCL IND O 140  X<=0? 141  GTO 05 142  RCL O 143  1 144  - 145  STO M 146  RCL 00 147  LAST X 148  - 149  INT 150  STO N 151  CLST 152   E99 153  LBL 06 154  RCL IND M  155  X<=0? 156  GTO 06 157  RCL IND N 158  X<>Y 159  / 160  X>Y? 161  GTO 06 162  RCL M 163  STO Z 164  LBL 06 165  CLX 166  SIGN 167  ST- M 168  RDN 169  DSE N 170  GTO 06 171  X<>Y 172  ENTER^ 173  X=0? 174  GTO 05 175  RCL M 176  INT 177  - 178  STO O 179  RCL IND Y 180  LBL 07 181  ST/ IND Y 182  ISG Y 183  GTO 07 184  RCL 00 185  INT 186  STO N 187  LBL 08 188  RCL 00 189  FRC 190  RCL N 191  + 192  RCL O 193  X=Y? 194  GTO 00 195  RCL M 196  RCL Z 197  + 198  RDN 199  RCL IND T 200  SIGN 201  LBL 09 202  CLX 203  RCL IND Y  204  LAST X 205  * 206  ST- IND Z 207  ISG Y 208  CLX 209  ISG Z 210  GTO 09 211  LBL 00 212  DSE N 213  GTO 08 214  GTO 04 215  LBL 10 216  RCL IND M 217  X<0? 218  GTO 00 219  X>0? 220  SF 01 221  RCL M 222  INT 223  LAST X 224  DSE Y 225  DSE X 226  CLX 227  INT 228   E3 229  / 230  + 231  1 232  STO O 233  LBL 11 234  RCL IND Y 235  1 236  X#Y? 237  GTO 11 238  ST- O 239  R^ 240  STO N 241  RDN 242  LBL 11 243  RDN 244  ABS 245  - 246  DSE Y 247  GTO 11 248  RCL O 249  X=0? 250  X#Y? 251  GTO 00 252  RCL N 253  RCL 00 254  INT 255  MOD 256    0 257  X<> IND Y 258  GTO 10 259  LBL 00 260  CLX 261  LBL 10 262  STO IND M 263  ISG M 264  GTO 10 265  RCL 00 266  0 267  LBL 12 268  RCL IND Y 269  STO IND Y 270  CLX 271  SIGN 272  + 273  ISG Y 274  GTO 12 275  CLX 276  LBL 13 277  X#0? 278  RCL IND a 279  X#0? 280  SF 02 281  PI 282  DSE a 283  GTO 13 284  RCL 00 285  CHS 286  STO 00 287  CLA 288  END

( 440 bytes / SIZE = 1+(p+1)(m+n"+p+1 )     with p = n+n'+n"

 STACK INPUTS OUTPUTS T n  = the number of inequations of the 1st kind ( ³ ) / Z n' = the number of equations ( = ) / Y n" = the number of inequations of the 2nd kind ( £ ) / X m = the number of unknowns max (F)

Example:

-Find  4 non-negative numbers  x, y, z and t satisfying:

56 ³  x + 2y + 3z +4t
41 ³ 5x +  y - 3z - t
61 ³ 4x + 7y +2z - 3t
25 = 2x + 3y - z + 2t                          ( always write the ">=" inequations first, then the "=" equations , and finally the "<=" inequations )
7  £   x + y + z + t
17 £ 2x - y + z + 3t

such that F = 2x + 4y +3z + 5t is maximized.

1- SIZE 092 ( or greater )

2- We store these 35 numbers in registers R01 thru R35 like this:

56   1   2   3   4                  R01  R08  R15  R22  R29
41   5   1  -3  -1                 R02  R09  R16  R23  R30
61   4   7   2  -3                 R03  R10  R17  R24  R31
25   2   3  -1   2      in        R04  R11  R18  R25  R32           respectively.    we must store 0  at the place of  F = R07 here
7    1   1   1   1                 R05  R12  R19   R26  R33
17   2  -1   1   3                 R06  R13  R20   R27  R34
0    2   4   3   5                 R07  R14  R21   R28  R35

3- There are 3 inequations of the first type, 1 equation , 2 inequations of the second type and 4 unknowns, therefore:

3 ENTER^
1 ENTER^
2 ENTER^
4 XEQ SPLEX  and ( 3mn49s later ) we obtain 76 in the X-register and in R00  ( F01 and F02 are cleared )

RCL 01 = 2
RCL 02 = 7
RCL 03 = 8
RCL 04 = 4

-Thus the maximum of F is 76 and it is achieved for x = 2, y = 7, z = 8, t = 4.

-We can also verify that:

R05 = 0 ; R06 = 52 ; R07 = 0    ( the n slack variables of the three first inequations )
R08 = R09 = R10 = 0               ( the n'+n" artificial variables which are always 0 when a solution exists )
R11 = 14 ; R12 = 0                  ( the n" slack variables of the two last inequations )

2°)  An Artificial-variable free Algorithm

The purpose is to find  n  non-negative real numbers:   x1 , ..... , xn   satisfying:

bi ³ ai;1x1 + ....... + ai;n xn            i = 1 , ....... , m                              ( m inequations )           the bi  , bi'  may be positive or not !
bi' = ai';1x1 + ....... +ai';n xn            i' = 1 , ..... , m'                              ( m' equations )

and such that      F  =  c1x1 + ....... +  cnxn       is maximized.    ( to find the minimum of F, seek the maximum of  -F )

-Here, we first treat each equation by pivoting on the largest coefficient to develop a complete basic variable set ( SF 00 ).

-Then, the dual simplex pivoting rule is used to transform all the negative bi into non-negative numbers ( if possible ) ( SF 01 )

•  pivot = ai,j < 0  with the largest  cj / ai,j

-Then, the primal simplex pivoting rule is applied to transform all the positive ci into non-positive numbers ( CF 01 )

•  pivot = ai,j > 0  with the smallest  bi / ai,j

-Finally, the tableau is solved.

Data Registers:                                               ( Registers R01 thru R(n+1)(m+m'+1) are to be initialized before executing "SPLX" )

•  Store all the coefficients in column order ( with 0 in Rm+m'+1 at the place of  F )

>>>>  When the program stops,  R00 = F maximum

R01 = x1 , .... , Rnn = xn    ;     Rnn+1 ...etc... = The slack variables.

Flags:  F00 & F01

>>>>   If  F01 is set at the end, the problem has no solution.

Subroutines: /

-Lines 140-218-220-221  are synthetic three-byte GTO's

 01  LBL "SPLX"   02  1   03  ST+ Z   04  +   05  ST+ Z   06  X<>Y   07  R^   08  STO M   09  +   10  ST* Y   11  ST* Z   12   E-5   13  *   14  STO 00   15  LASTX   16  +   17  X<> Z   18   E3   19  /   20  ST+ 00   21  ISG 00   22  +   23  ISG X   24  X=0?   25  GTO 00   26  CLRGX   27  +   28  SIGN   29  LBL 16   30  STO IND L   31  ISG L   32  GTO 16   33  LBL 00   34  RCL 00   35  INT   36  RCL M   37   E3   38  /   39  +   40  STO Q   41  SF 00   42  SF 01   43  LBL 13   44  DSE Q   45  GTO 13   46  CF 00   47  GTO 01   48  LBL 13 49  RCL Q   50  INT   51  STO O   52  RCL 00   53  FRC   54  +   55  ISG X   56  STO M   57  CLST   58  LBL 14   59  RCL IND M   60  ABS   61  X<=Y?   62  GTO 14   63  RCL M   64  STO Z   65  LBL 14   66  RDN   67  ISG M   68  GTO 14   69  1/X   70  X<>Y   71  GTO 00   72  LBL 01   73  RCL 00   74  INT   75  STO O   76  LBL 02   77  DSE O   78  GTO 02   79  CF 01   80  GTO 04   81  LBL 02   82  RCL IND O   83  CHS   84  X<=0?   85  GTO 02   86  RCL O   87  RCL 00   88  FRC   89  +   90  ISG X   91  STO M   92  RCL 00   93  ISG X   94  STO N   95  CLST   96   E99 97  LBL 03   98  RCL IND M   99  CHS 100  X<=0? 101  GTO 03 102  RCL IND N 103  X<>Y 104  / 105  X>Y? 106  GTO 03 107  RCL M 108  STO Z 109  LBL 03 110  ISG M 111  RDN 112  ISG N 113  GTO 03 114  X<> Z 115  LBL 00 116  ENTER^ 117  ENTER^ 118  X=0? 119  GTO 02 120  RCL 00 121  INT 122  MOD 123  - 124  STO M 125  CLX 126  RCL O 127  RCL 00 128  FRC 129  + 130  GTO 00 131  LBL 04 132  RCL 00 133  STO O 134  LBL 05 135  ISG O 136  GTO 05 137  RCL 00 138  ISG X 139  STO M 140  GTO 10 141  LBL 05 142  RCL IND O 143  X<=0? 144  GTO 05 145  RCL O 146  1 147  - 148  STO M 149  RCL 00 150  LASTX 151  - 152  INT 153  STO N 154  CLST 155   E99 156  LBL 06 157  RCL IND M 158  X<=0? 159  GTO 06 160  RCL IND N 161  X<>Y 162  / 163  X>Y? 164  GTO 06 165  RCL M 166  STO Z 167  LBL 06 168  CLX 169  SIGN 170  ST- M 171  RDN 172  DSE N 173  GTO 06 174  X<>Y 175  ENTER^ 176  X=0? 177  GTO 05 178  RCL M 179  INT 180  - 181  LBL 00 182  STO P 183  RCL IND Y 184  LBL 07 185  ST/ IND Y 186  ISG Y 187  GTO 07 188  RCL 00 189  INT 190  STO N 191  LBL 08 192  RCL 00 193  FRC 194  RCL N 195  + 196  RCL P 197  X=Y? 198  GTO 08 199  RCL M 200  ST+ L 201  CLX 202  RCL IND L 203  SIGN 204  LBL 09 205  CLX 206  RCL IND Y 207  LASTX 208  * 209  ST- IND Z 210  ISG Y 211  CLX 212  ISG Z 213  GTO 09 214  LBL 08 215  DSE N 216  GTO 08 217  FS? 00 218  GTO 13 219  FS? 01 220  GTO 01 221  GTO 04 222  LBL 10 223  RCL IND M 224  X<0? 225  GTO 00 226  X>0? 227  SF 01 228  RCL M 229  INT 230  LASTX 231  DSE Y 232  DSE X 233  CLX 234  INT 235   E3 236  / 237  + 238  1 239  STO O 240  LBL 11 241  RCL IND Y 242  X#Y? 243  X=0? 244  FS? 30 245  GTO 00 246  X#Y? 247  GTO 11 248  ST- O 249  ENTER^ 250  + 251  LBL 11 252  RDN 253  DSE Y 254  GTO 11 255  RCL O 256  X#0? 257  GTO 00 258  R^ 259  RCL 00 260  INT 261  MOD 262  0 263  X<> IND Y 264  GTO 10 265  LBL 00 266  CLX 267  LBL 10 268  STO IND M 269  ISG M 270  GTO 10 271  RCL IND 00 272  CHS 273  X<> 00 274  1 275  ISG Y 276  LBL 12 277  RCL IND Y 278  STO IND Y 279  X<0? 280  SF 01 281  CLX 282  SIGN 283  + 284  ISG Y 285  GTO 12 286  RCL 00 287  CLA 288  END

( 446 bytes / SIZE 1+(n+m+1)(m+m'+1) )

 STACK INPUTS OUTPUTS Z m  = the number of inequations  ( ³ ) / Y m' = the number of equations ( = ) / X n = the number of unknowns max (F)

Example1:   The same one.

Find  4 non-negative numbers  x, y, z and t satisfying:

56 ³  x + 2y + 3z +4t
41 ³ 5x +  y - 3z - t
61 ³ 4x + 7y +2z - 3t
25 = 2x + 3y - z + 2t
7 £   x + y + z + t
17 £ 2x - y + z + 3t

such that F = 2x + 4y +3z + 5t is maximized.

1-SIZE 071  or more.  Then re-write the LP as follows:  the "³" inequations first ( if any ), then the "=" equations ( if any )

56 ³  x + 2y + 3z +4t
41 ³ 5x +  y - 3z - t
61 ³ 4x + 7y +2z - 3t
-7 ³  -x  -  y  -  z  -  t                          the "£" inequations must be multiplied by (-1) to become  "³" inequations
-17 ³ -2x + y  -  z - 3t
25 =  2x + 3y - z + 2t

2- We store these 35 numbers in registers R01 thru R35 like this:

56   1   2   3   4                  R01  R08  R15  R22  R29
41   5   1  -3  -1                 R02  R09  R16  R23  R30
61   4   7   2  -3                 R03  R10  R17  R24  R31
-7  -1  -1  -1  -1     in        R04  R11  R18  R25  R32           respectively.   ( we must store 0  at the place of  F that is in R07 here )
-17 -2   1  -1  -3                R05  R12  R19   R26  R33
25   2   3  -1   2                 R06  R13  R20   R27  R34
0   2   4   3   5                 R07  R14  R21   R28  R35

3- There are m = 5 inequations, m' = 1 equation and n = 4 unknowns, therefore:

5  ENTER^
1  ENTER^
4  XEQ "SPLX"   >>>>   Max = 76 = R00                  ---Execution time = 2mn45s---                ( instead of 3mn49s with the big-M method )

R01 = x = 2
R02 = y = 7
R03 = z = 8
R04 = t = 4       thus the maximum of F is 76 and it is achieved for x = 2, y = 7, z = 8, t = 4.

We can also verify that the 5 slack variables are

R05 = 0 , R06 = 52 , R07 = 0 ,  R08 = 14 ,  R09 = 0

Example2:   Almost the same one: same constraints but find the minimum of  F = 2x + 4y +3z + 5t

-Simply change all the signs in the corresponding raw i-e the last one.  Store

56   1   2   3   4                  R01  R08  R15  R22  R29
41   5   1  -3  -1                 R02  R09  R16  R23  R30
61   4   7   2  -3                 R03  R10  R17  R24  R31
-7  -1  -1  -1  -1     in        R04  R11  R18  R25  R32           respectively.
-17 -2   1  -1  -3                R05  R12  R19   R26  R33
25   2   3  -1   2                 R06  R13  R20   R27  R34
0  -2  -4  -3  -5                 R07  R14  R21   R28  R35

5  ENTER^
1  ENTER^
4  XEQ "SPLX"   >>>>   Max = -30.62295079 = R00   ( F01 is clear )          ---Execution time = 3mn16s---

-So  Fmin = 30.62295079   and it is achieved for

R01 = x = 7.967213120                 and the 5 slack variables are:
R02 = y = 2.278688521
R03 = z = 0                                    R05 = 39.01639346      R07 = 16.52459018
R04 = t = 1.114754093                  R06 = 0      R08 = 4.360655733     R09 = 0

Notes:

-We assume that all the equations are linearly independant.
-Especially, don't key in more equations ( = ) than unknowns.
-Line 69 ( 1/X ) is useful to produce a DATA ERROR message in this case, but roundoff-errors could fool this trick.

-If there is no equation and if all bi's are positive, both programs are equivalent.
-Otherwise, this method obviously uses less data registers, the program often runs faster and roundoff-errors are usually smaller !
-Moreover, synthetic register a is unused.

-If, however, there are not enough data registers, we can use the equations to sustitute one or several variables
into the inequations and thus reduce the size of the problem.

Simplification:

-If you never solve problems that involve equations:

delete lines 218-217-115 and replace lines 02 to 71 by

RCL Y    *                ISG X         FRC            SIGN
1             ST+ Y        CLRGX      ISG X         LBL 13
ST+ Z     X<>Y        X<>Y         STO 00       STO IND L
+             E3              E5              LASTX        ISG L
STO T    /                  /                  E-5             GTO 13
ST* Z     +                 +                +                 SF 01

-Then the inputs are:

 STACK INPUTS OUTPUTS Y m  = the number of inequations  ( ³ ) / X n = the number of unknowns max (F)

-Actually, this simplified version may also be used if there are equations:
-Replace them by 2 inequations, for example, instead of  25 =  2x + 3y - z + 2t  write

25  ³  2x + 3y - z + 2t
-25  ³ -2x  - 3y + z - 2t

-However, it might then require more data registers than the big-M method !

References:

  Francis Scheid - "Numerical analysis" - McGraw-Hill    ISBN:  07-055197-9
  Hossein Arsham - "Artificial-variable Free Solution Algorithms for LP Models"