From: george@unixg.ubc.ca (George chow)
Subject: Here's my implementation of the Simplex Method
Followup-To: comp.sys.hp48
Originator: george@unixg.ubc.ca
Keywords: simplex method, linear programming
Organization: University of British Columbia, Vancouver, B.C., Canada
Date: Wed, 15 Apr 1992 23:09:46 GMT


The following directory contains a set of routines for the Primal Simplex
Algorithm with full tableau. The routines allow one to manually step through
the method by choosing each pivot or to have the automatic procedure execute
the entire method for you. (I haven't tried mixing the two...)

For the dual problem

  P = max { c^T x : A x <= b, x >= 0 }
  D = min { y^T b : y^T A >= c^T, y >= 0 }

this implemenation requires an initial tableau with the following layout:

	+---+---------+---------+
	|   |         |         |
	|   |         |         |
	| b |    A    |    I    |
	|   |         |         |
	|   |         |         |
	+---+---------+---------+
	| 0 |  -c^T   |    0    |
	+---+---------+---------+

After each step, a copy of the tableau can be returned to be stack. The
returned tableau has layout:

	+---+---------+---------+
	|   |         |         |
	| ^ |         |         |
	| b |    U    |    B    |
	|   |         |         |
	|   |         |         |
	+---+---------+---------+
	| a |    n    |    y    |
	+---+---------+---------+

For instance, the problem:

  max { 3 x1 + 2 x2 } subject to
  
    x1 + 2 x2 <= 3
  2 x1 -   x2 <= 2
   -x1 + 2 x2  = 1
    x1        >= 0
           x2 >= 0

has initial tableau:

       x1 x2  e1 e2 e3
 [[ 3   1  2   1  0  0 ]
  [ 2   2 -1   0  1  0 ]
  [ 1  -1  2   0  0  1 ]
  [ 0  -3 -2   0  0  0 ]]

where e3 is artificial.

The algorithm is split into 3 phases:
  Phase 0: drive all artificial variable out of the basis
  Phase 1: eliminate all negative components of b.
  Phase 2: eliminate all unwanted negative components of n and y.
  
Here is a description of each routine:

INIT: This must be run before any of the other routines will work correctly.
      INIT precomputes all of the dimensions and indices needed by the other
      routine. It takes as argument from the stack the initial tableau and the
      index of the artificial variable block. One must keeps
      all the artificial variables to the extreme right hand side. The index
      that is provided is the beginning of the artificial variable sections.
      In the above example, we would have an index of 3.
      
FIND.P0.COLUMN: takes a row index from the stack and returns the index of the
                first non-zero element of that row or 0 if none. This is used
	        for picking the pivotal column in Phase 0.
		
                                                                ^
\GB.STATE: returns the index of the first negative component of b or 0 if none.

FIND.TARGET.COLUMN: takes a row index from the stack and returns the index of
                    the first negative element of that row or 0 if none. This
		    is used for picking the pivotal column in Phase 1 and 2.
		    
GEN.RATIO: generates the list of ratios b[i]/A[i][j] for the argument row i.
           These ratios are used with MIN.CRITERION to pick the pivotal row.
	   
MIN.CRITERION: returns the index of the element in the argument list with the
               small value.

PRIMAL.STEP: takes as argument a row and a column index. Performs the required
             row operation on the tableau.
             
DOIT: This starts the automatic Simplex method. This will step thru all three
      phases of the Simplex method and return the final tableau.
  
      
Full Example:
  1. Download the directory to your hp48sx
  2. CRDIR a new directory inside SIMPLEX.EXPLICIT called TEST
  3. Change into the directory TEST
  4. Key in the above matrix and store as INIT.TABLEAU
  5. Put INIT.TABLEAU and 3 on the stack
  6. Set user-flag 1 (1 SF) and press [CST]
  7. Execute INIT
  8. You should get a copy of INIT.TABLEAU on the stack
  
  9. Since we have an artificial variable, we need to do Phase 0:
 10. Do: 
        3 2 STEP
     You should get back:
        1: [[4 0 4 1 0 1]
            [4 0 3 0 1 2]
	    [-1 1 -2 0 0 -1]
	    [-3 0 -8 0 0 -3]]
     You have just pivoted with row 3 and column 2 (which is a -1)
     Since we only have one artifical variable, we're done in Phase 0 and can
     proceed to Phase 1.
 11. We want to use the -2 in row 3 as our pivot. Do:
        3 3 STEP
     You should get back:
        1: [[2 2 0 1 0 -1]
            [2.5 1.5 0 0 1 .5]
	    [.5 -.5 1 0 0 .5]
	    [1 -4 0 0 0 1]]
	                         ^
 12. We now have all elements of b positive so Phase 1 is done. We proceed to
     Phase 2. We take the second 2 in row 1 as our pivot. Do:
        1 2 STEP
     You should get back:
        1: [[1 1 0 .5 0 -.5]
            [1 0 0 -.75 1 1.25]
	    [1 0 1 .25 0 .25]
	    [5 0 0 2 0 -1]]
	    
	               ^
 13. Since we now have b positive and all needed elements of n and y positive,
     we're done. Our optimal program is:
         x = [1 1]
	 y = [2 0 -1]
	 n = [0 0]
	 e = [0 1 0]
     with optimal value 5.
     
 14. Using the automatic routine gives us the same result. Do:
         INIT.TABLEAU INIT DOIT
     You should get three tableaux, the final of which gives you the same
     optimal program.
  

Notes: 
  1. There is no smarts in DOIT. Conceivably, it can get trapped in a 
     cycle very easily. I wrote this as an aid for my homework, not to solve
     `real' optimization problems.
  2. This implementation isn't very fast. A typical homework problem with a
     5x7 tableau takes about 70-90 seconds for DOIT to process.
  3. Generating the initial tableau is a hassle but I haven't found a fast
     way for the machine to do it. Where possible, start typing into your HP
     ASAP. ;)
  4. If you have any comments/bug-reports/problems, please email me. 
  5. I'm working on the implicit method.


====
Here's the uuencoded directory:

Checksum: #F06Eh
Size: 2241

begin 600 SIMPLEX.EXP
M2%!(4#0X+426*O!_$@\````"14<"Z"F0,P`S*2````4`H````````````)@`
M````````D0```````````````````)$``````````0````````"2````````
M``$`````````````````````````````````````````"```````````````
M``````$``````````@```````````````````),````````````````````!
M``````````````````````$``````!`!``````````$``````````P``````
M```!``````````(``````````0```````````````````````````````0``
M`````````````````)0`````````D0`````````!`````````)$`````````
M!`````````"1``````````````````````````````````````````$`````
M``````````````"5````````````````````DP`````````"`````````)<`
M`````````````````````````````````````````$8#,#`T135`IP)T*L"B
M`@T`,$55!(7D`@M04DE-04PN4U1%4"LQP*("!0#`H@(%`,"B`@4`@.0"!$1/
M251(+D"0Y)1$M1(#B0!`D.241$70V0(>-D)E(T@N<!`D1>5B%2259R/-##+L
M(LFB,C$<^BYRN!_5+W*X'U0V@N0"!DU!5%))6'DVTLP@N,FQ]QS8^T%E(T@N
M,##TQ)1G(\T,<K@?5#:"Y`(#4D]7>3;2S"#)HI+0&E0V@N0"!D-/3E-44GDV
MTLP@2"YP$"1%Y6(5)(7D`@-#3TQ(+F`P].0T126%Y`('05)4+E9!4@FMD=`:
MR:*2T!I(+C`P],3D/QI4-H+D`@=605(N14Y$>3;2S"`Y-K(2`Y4!L``EE=04
MQ.0R154$M=#9`AXV(KH?WJ(R>!Q4-H+D`@9-051225AY-M*['\;702LJP332
MY@("<')M+B``-];F`@5P:79O=&TN8&`6-D;W)N=A(\FB@N0"`U)/5Z`QTN8"
M`6G#+M+F`@%I;2X@`"?7J1Y4-H+D`@9-051225AY-D*G`FTN$)#6Y@("<&,K
M,6!\'6TN4`"79O=&5_`:A_M!92-M+F!@%C9&]R:79R/-#$(K*IWJ,7@>^B[2
MV0+)HH+D`@-#3TR@,=+F`@%J5#:"Y`(&34%44DE8>39"IP)M+A"0UN8"`6HK
M,2"Z'\;7T>8"!F9A8W1O<E0V@N0"!DU!5%))6'DV0J<";2X@`"?7Y@(!:BLQ
M8'P=[JV1T!H'U,$D(RLQ4/TB3#*2+"I(+C`P],0$&B-M+A"@1F4C2"Y@T!1$
M)96$E6<C="K0Y@("<')M+A"@MA(#HOMA?!UM+E``EV;W1E?P&@?4P20CR:*"
MY`(#4D]7H#'2Y@(!:50V@N0"!DU!5%))6'DV0J<";2X0D-;F`@)P8RLQT.8"
M`6EM+B``)Q=V'@?4P20CPRZ2+"H3PZ'O(D@N8-`41"65A%7](OXUDF,C*S'@
M.``-34E.+D-2251%4DE/3@V=+>!A(WO/@;T?M*+RK1K!--+F`@)I>&TN$&#G
M82-(+F`P].0T1265+"J@,=+F`@%IPRYRN!]M+A!@Y[L>^B[2V0)4-M+F`@%V
M>3;2S"!M+A"01F4C;2X@D(:79R/-#+(2`[4O@KT?U2]B."J`,]+F`@)I>/XU
MDF,C*S%0$0`)1T5.+E)!5$E/"9TMX&$CP332Y@(!8QXVDBPJ2"Y@,/3D-$4E
M!1HC;2X0D/8](YTM0&4C2"Y@T!1$)96$E6<C="K0Y@(!:<FBLA(#QM=!92-(
M+F#0%$0EE8259R-T*M#F`@%I;2X0,+82`\;74?`:*S'P<2.=+3"_']^JL1(#
MM2_2V0+#+G*X'[2BXKL>^B[2V0+8^_&M&BLQ4/TB*S%0_2),,H+D`@9#3TY3
M5%*#Q^%?(SDVLA(#<`$@893D1.1"%21U5$3E,O3$5-7D)-'9`AXVXBTJP332
MY@(#<F]W;2X0H.9A([2B,@PCPRY"92-(+F#0%$0EE8259R-T*M#F`@-R;W=M
M+A"@MA(#QM=!*RJ^ZZ'O(ITM@+T?;2X0H+82`]4OT@XCA_M!*RI=[$%E(VTN
M$*"69R/T"(+D`@=605(N14Y$7>R1@!ZY-N)?(SDVLA(#,P%P\.TR115$5730
MV0(>-D(K*LFB$DPC;2X0D.9A(\,P,NPB5#:"Y`(&34%44DE8>39"IP)M+A"0
MEBPJ*S%@?!VTHN*['OHNTMD"V/O1Y@(!:2LQ4/TB[3!RN!^THM+%'E0VTN8"
M`6EY-D*/($@N8##TY#1%)=7%'@GHD6LC_C628R,K,6`0``Y&24Y$+E`P+D-/
M3%5-3@Z=+>!A(]ZB$DPC;2XP(/=VU^8"`6H>-D(K*L,P,NPB5#:"Y`(&34%4
M4DE8>39"IP)M+C`@]W;7Y@(!:BLQ8'P=M*+2J1[Z+M+9`MC[T>8"`6HK,5#]
M(NTP<K@?M*+2Q1Y4-M+F`@%J>39"CR!(+G!@%23E4N1$U,4>">B1:R/^-9)C
M(RLQL!(`!$1/250$G2W@82/#+H+D`@=!4E0N5D%2^B[2V0)(+G`0)$7E8A4D
MA>0"!D-/3E-44J`QTN8"`6EM+A"0=K@?2"[@8)3D1.0"!>,R],14U>2$Y`(+
M4%))34%,+E-415!,,K(2`]4O,@,C2"YP\.TR115$572X'UTPTMD"2"X@893D
M1.1"%21U5$3E,O3$5-7D=+@?2"Z0<%3DY"(51)7TA.0"#4U)3BY#4DE415))
M3TZ]^X'D`@M04DE-04PN4U1%4"LQ0&DCV/LQ`R-(+C`@]72%Y`(21DE.1"Y4
M05)'150N0T],54U.A_O1!2.=+7"X'T@ND'!4Y.0B%425](3D`@U-24XN0U))
D5$5224].O?N!Y`(+4%))34%,+E-415`K,4!I(]C[D6,C*S$`
`
end

  
Newsgroups: comp.sys.hp48