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:   x₁ , ..... , x_m  satisfying:

          b_i³ a_i;1x₁ + ....... + a_i;m x_m            i = 1 , .... , n                                ( n inequations )                   all the b_i  , b_i'  and  b_i"
          b_i' = a_i';1x₁ + ....... +a_i';m x_m            i' = 1 , ..... , n'                              ( n' equations )                    must be non-negative*
          b_i"£ a_i";1x₁+ ....... +a_i";m x_m          i" = 1 , ..... , n"                             (n" inequations )

and such that      F  =  c₁x₁ + ....... +  c_mx_m       is maximized.    ( to find the minimum of F, seek the maximum of  -F )
 

* if some b_i 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 = x₁ , .... , Rmm = x_m
        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 c_j  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 )
 
 

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

 
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:   x₁ , ..... , x_n   satisfying:

                          b_i ³ a_i;1x₁ + ....... + a_i;n x_n            i = 1 , ....... , m                              ( m inequations )           the b_i  , b_i'  may be positive or not !
                          b_i' = a_i';1x₁ + ....... +a_i';n x_n            i' = 1 , ..... , m'                              ( m' equations )

and such that      F  =  c₁x₁ + ....... +  c_nx_n       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 b_i into non-negative numbers ( if possible ) ( SF 01 )

       •  pivot = a_i,j < 0  with the largest  c_j / a_i,j

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

       •  pivot = a_i,j > 0  with the smallest  b_i / a_i,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 = x₁ , .... , Rnn = x_n    ;     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
 
 

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

 
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  F_min = 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 b_i'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:
 
 

 
-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:

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