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

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

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

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:

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

hp41programs

Linear Programming: the Simplex method for the HP-41