f10-8

10.8LinearProgrammingandtheSimplexMethod 423Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).rather a triangular decomposition of A, itsCholesky decomposition (cf.§2.9). The updating formula used for the Cholesky decomposition of Ais of order N2and can be arranged to guarantee that the matrix remains positive definite and nonsingular, even in the presence offinite roundoff. This method is due to Gill and Murray [1,2]. CITED REFERENCES AND FURTHER READING: Dennis,J.E., andSchnabel,R.B. 1983, NumericalMethods forUnconstrained Optimizationand Nonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall). [1] Jacobs,D.A.H.(ed.)1977, TheStateoftheArtinNumericalAnalysis (London:AcademicPress), Chapter III.1, §§3–6 (by K. W. Brodlie). [2] Polak,E.1971, ComputationalMethodsinOptimization (NewYork:AcademicPress),pp.56ff.[3] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 467–468. 10.8 Linear Programming and the Simplex Method The subject of linear programming , sometimes called linear optimization , concernsitselfwiththefollowingproblem: For Nindependentvariables x1,...,x N, maximize the function z=a01x1+a02x2+···+a0NxN (10.8.1 ) subject to the primary constraints x1≥0,x 2≥0, ... x N≥0( 10.8.2 ) and simultaneously subject to M =m1+m2+m3additional constraints, m1of them of the form ai1x1+ai2x2+···+aiNxN≤bi (bi≥0) i=1 ,...,m 1 (10.8.3 ) m2of them of the form aj1x1+aj2x2+···+ajNxN≥bj≥0 j=m1+1 ,...,m 1+m2(10.8.4 ) and m3of them of the form ak1x1+ak2x2+···+akNxN=bk≥0 k=m1+m2+1 ,...,m 1+m2+m3(10.8.5 ) The various aij’s can have either sign, or be zero. The fact that the b’s must all be nonnegative (as indicated by the final inequality in the above three equations) is a matter of convention only, since you can multiply any contrary inequality by −1. There is no particular significance in the number of constraints Mbeing less than, equal to, or greater than the number of unknowns N. 424 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Asetofvalues x1...x Nthatsatisfiestheconstraints(10.8.2)–(10.8.5)iscalled afeasiblevector . Thefunctionthatwearetryingtomaximizeis calledthe objective function. The feasible vector that maximizes the objective function is called the optimal feasible vector . An optimal feasible vector can fail to exist for two distinct reasons: (i)thereare nofeasiblevectors,i.e.,thegivenconstraintsareincompatible, or (ii) there is no maximum, i.e., there is a direction in Nspace where one or more of the variables can be taken to infinity while still satisfying the constraints, giving an unbounded value for the objective function. As you see, the subject of linear programmingis surroundedby notational and terminological thickets. Both of these thorny defenses are lovingly cultivated by acoterie of stern acolytes who have devoted themselves to the field. Actually, the basic ideas of linear programming are quite simple. Avoiding the shrubbery, we want to teach you the basics by means of a couple of specific examples; it shouldthen be quite obvious how to generalize. Why is linear programming so important? (i) Because “nonnegativity” is the usual constraint on any variable x ithat represents the tangible amount of some physical commodity, like guns, butter, dollars, units of vitamin E, food calories, kilowatt hours, mass, etc. Hence equation (10.8.2). (ii) Because one is often interested in additive (linear) limitations or bounds imposed by man or nature: minimumnutritionalrequirement,maximumaffordablecost,maximumonavailable labor or capital, minimum tolerable level of voter approval, etc. Hence equations(10.8.3)–(10.8.5). (iii) Because the function that one wants to optimize may be linear, or else may at least be approximatedby a linear function — since that is the problem that linear programming cansolve. Hence equation (10.8.1). For a short, semipopular survey of linear programming applications, see Bland [1]. Here is a specific example of a problem in linear programming, which has N=4,m1=2,m2=m3=1, hence M =4: Maximize z=x1+x2+3 x3−1 2x4 (10.8.6 ) with all the x’s nonnegative and also with x1+2 x3≤740 2x2−7x4≤0 x2−x3+2 x4≥1 2 x1+x2+x3+x4=9(10.8.7 ) Theanswerturnsouttobe(to2decimals) x1=0,x2=3 .33,x3=4 .73,x4=0 .95. In the rest of this section we will learn how this answer is obtained. Figure 10.8.1 summarizes some of the terminology thus far. 10.8LinearProgrammingandtheSimplexMethod 425Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).additional constraint (inequality) additional constraint (inequality)the optimal feasible vectorsome feasible vectorsx1 primary constraint x2a feasible basic vector (not optimal)primary constraintadditional constraint (equality) z = 3.1 z = 2.9z = 2.8z = 2.7z = 2.6z = 2.5z = 2.4z = 3.0 Figure 10.8.1. Basic concepts of linear programming. The case of only two independent variables, x1,x2, is shown. The linear function z, to be maximized, is represented by its contour lines. Primary constraints require x1and x2to be positive. Additional constraints may restrict the solution to regions (inequality constraints) or to surfaces of lower dimensionality (equality constraints). Feasible vectorssatisfy all constraints. Feasible basic vectors also lie on the boundary of the allowed region. The simplex method steps among feasible basic vectors until the optimal feasible vector is found. FundamentalTheorem of LinearOptimization Imaginethatwestartwithafull N-dimensionalspaceofcandidatevectors. Then (inmind’seye,atleast)wecarveawaytheregionsthatareeliminatedinturnbyeach imposed constraint. Since the constraints are linear, every boundary introduced by thisprocessisaplane,orratherhyperplane. Equalityconstraintsoftheform(10.8.5)force the feasible region onto hyperplanes of smaller dimension, while inequalities simply divide the then-feasible region into allowed and nonallowedpieces. When all the constraints are imposed, either we are left with some feasible region or else there are no feasible vectors. Since the feasible region is boundedby hyperplanes,it is geometricallya kind of convexpolyhedronor simplex (cf. §10.4). If there is a feasible region, can the optimal feasible vector be somewhere in its interior, away from the boundaries? No, because the objective function is linear. This means that it always has a nonzero vector gradient. This, in turn, means thatwe could always increase the objective function by running up the gradient until we hit a boundary wall. Theboundaryofanygeometricalregionhasonelessdimensionthanitsinterior. Therefore,wecannowrunupthegradientprojectedintotheboundarywalluntilwe 426 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).reach an edge of that wall. We can then run up that edge, and so on, down through whatever number of dimensions, until we finally arrive at a point, a vertexof the original simplex. Since this point has all Nof its coordinates de fined, it must be the solution of Nsimultaneous equalities drawn from the original set of equalities and inequalities (10.8.2) –(10.8.5). Points that are feasible vectors and that satisfy Nof the original constraints as equalities, are termed feasible basic vectors .I f N> M, then a feasible basic vector has at least N−Mof its componentsequal to zero, since at least that many of the constraints (10.8.2) will be needed to make up the total of N. Put the other way,at most Mcomponentsof a feasible basic vector are nonzero. In the example (10.8.6)–(10.8.7),you can check that the solution as given satis fies as equalities the lastthreeconstraintsof(10.8.7)andtheconstraint x1≥0,fortherequiredtotalof4. Put together the two preceding paragraphs and you have the Fundamental Theorem of LinearOptimization : If an optimalfeasible vectorexists, then thereis a feasible basic vector that is optimal. (Didn ’t we warn you about the terminological thicket?) The importance of the fundamental theorem is that it reduces the optimization problem to a “combinatorial ”problem, that of determining which Nconstraints (out of the M+Nconstraints in 10.8.2 –10.8.5)should be satis fied by the optimal feasible vector. We haveonly to keep tryingdifferentcombinations,andcomputing the objective function for each trial, until we find the best. Doing this blindly would take halfway to forever. The simplex method ,first published by Dantzig in 1948 (see [2]), is a way of organizingthe procedureso that (i)aseriesofcombinationsistriedforwhichtheobjectivefunctionincreasesateachstep, and (ii) the optimal feasible vector is reached after a number of iterations that isalmostalwaysnolargerthanoforder MorN,whicheveris larger. Aninteresting mathematicalsidelightisthatthissecondproperty,althoughknownempiricallyever sincethesimplexmethodwasdevised,wasnotprovedtobetrueuntilthe1982work of Stephen Smale. (For a contemporary account, see [3].) SimplexMethod fora Restricted NormalForm A linear programming problem is said to be in normal form if it has no constraintsintheform(10.8.3)or(10.8.4),butratheronlyequalityconstraintsofthe form (10.8.5) and nonnegativity constraints of the form (10.8.2). Forourpurposesitwillbeusefultoconsideranevenmorerestrictedsetofcases, with this additional property: Each equality constraint of the form (10.8.5) must haveat least onevariablethathasa positivecoef ficientand thatappearsuniquelyin that one constraint only . We can then choose one such variable in each constraint equation, and solve that constraint equation for it. The variables thus chosen arecalledleft-hand variables orbasic variables , and there are exactly M(=m 3)o f them. The remaining N−Mvariables are called right-handvariables ornonbasic variables. Obviously this restricted normal form can be achieved only in the case M≤N, so that is the case that we will consider. You may be thinking that our restricted normal form is so specialized that it is unlikely to include the linear programming problem that you wish to solve. Not at all! We will presently show how anylinear programming problem can be 10.8LinearProgrammingandtheSimplexMethod 427Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).transformed into restricted normal form. Therefore bear with us and learn how to apply the simplex method to a restricted normal form. Here is an example of a problem in restricted normal form: Maximize z=2 x2−4x3 (10.8.8 ) with x1,x2,x3, and x4all nonnegative and also with x1=2−6x2+x3 x4=8+3 x2−4x3(10.8.9 ) This example has N=4,M =2; the left-hand variables are x1and x4; the right-hand variables are x2and x3. The objective function (10.8.8) is written so as to depend only on right-hand variables; note, however, that this is not an actual restriction on objective functions in restricted normal form, since any left-handvariables appearing in the objective function could be eliminated algebraically by use of (10.8.9) or its analogs. For any problemin restricted normal form, we can instantly read off a feasible basic vector(althoughnot necessarilythe optimalfeasiblebasic vector). Simplyset all right-handvariables equal to zero, and equation(10.8.9)then gives the values oftheleft-handvariablesforwhichtheconstraintsaresatis fied. Theideaofthesimplex method is to proceed by a series of exchanges. In each exchange, a right-hand variableandaleft-handvariablechangeplaces. Ateachstagewemaintainaproblemin restricted normal form that is equivalent to the original problem. It is notationally convenient to record the information content of equations (10.8.8) and (10.8.9) in a so-called tableau, as follows: x2 x3 z 0 2 −4 x1 2 −6 1 x4 8 3 −4(10.8.10 ) You should study (10.8.10) to be sure that you understand where each entry comes from,and how to translate backand forth betweenthe tableauand equationformats of a problem in restricted normal form. Thefirst step in the simplex method is to examine the top row of the tableau, whichwewillcallthe “z-row.”Lookattheentriesincolumnslabeledbyright-hand variables(wewill callthese “right-columns ”). We wanttoimagineinturntheeffect of increasing each right-hand variable from its present value of zero, while leaving all the other right-hand variables at zero. Will the objective function increase ordecrease? The answer is given by the sign of the entry in the z-row. Since we want to increase the objective function, only right columns having positive z-row entries are ofinterest. In (10.8.10)thereis onlyonesuchcolumn,whose z-rowentryis 2. The second step is to examine the column entries below each z-row entry that was selected by step one. We want to ask how muchwe can increase the right-hand variablebeforeoneoftheleft-handvariablesisdrivennegative,whichisnotallowed. If the tableau element at the intersection of the right-handcolumn and the left-hand 428 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).variable’s row is positive, then it poses no restriction: the corresponding left-hand variablewilljustbedrivenmoreandmorepositive. If alltheentriesinanyright-hand column are positive, then there is no bound on the objective function and (having said so) we are done with the problem. If one or more entries below a positive z-row entry are negative, then we have tofigure out which such entry first limits the increase of that column ’s right-hand variable. Evidentlythelimitingincreaseisgivenbydividingtheelementintheright- hand column (which is called the pivot element ) into the element in the “constant column”(leftmost column) of the pivot element ’s row. A value that is small in magnitude is most restrictive. The increase in the objective function for this choiceofpivot elementis thenthat valuemultipliedby thez-rowentryof that column. We repeat this procedure on all possible right-hand columns to find the pivot element with the largest such increase. That completes our “choiceof a pivot element. ” In the above example, the only positive z-row entry is 2. There is only one negativeentrybelowit,namely −6,sothisisthepivotelement. Itsconstant-column entryis 2. Thispivotwillthereforeallow x 2tobeincreasedby 2÷|6|,whichresults in an increase of the objective function by an amount (2×2)÷|6|. The third step is to dothe increase of the selected right-hand variable, thus makingit aleft-handvariable;andsimultaneouslytomodifytheleft-handvariables, reducingthepivot-rowelementtozeroandthusmakingitaright-handvariable. For our above example let ’s do this first by hand: We begin by solving the pivot-row equation for the new left-hand variable x2in favor of the old one x1, namely x1=2−6x2+x3 → x2=1 3−1 6x1+1 6x3 (10.8.11 ) We then substitute this into the old z-row, z=2 x2−4x3=2/bracketleftbig1 3−1 6x1+1 6x3/bracketrightbig −4x3=2 3−1 3x1−11 3x3 (10.8.12 ) and into all other left-variable rows, in this case only x4, x4=8+3/bracketleftbig1 3−1 6x1+1 6x3/bracketrightbig −4x3=9−1 2x1−7 2x3 (10.8.13 ) Equations (10.8.11) –(10.8.13) form the new tableau x1 x3 z2 3−1 3−11 3 x21 3−1 61 6 x4 9 −1 2−7 2(10.8.14 ) Thefourthstepistogobackandrepeatthe firststep,lookingforanotherpossible increaseoftheobjectivefunction. Wedothisasmanytimesaspossible,thatis,until all theright-handentriesin thez-rowarenegative,signalingthat nofurtherincreaseis possible. Inthepresentexample,thisalreadyoccursin(10.8.14),so weare done. The answer can now be read from the constant column of the final tableau. In (10.8.14) we see that the objective function is maximized to a value of 2/3for the solution vector x 2=1 /3,x4=9,x1=x3=0. Nowlookbackovertheprocedurethatledfrom(10.8.10)to(10.8.14). Youwill find that it could be summarized entirely in tableau format as a series of prescribed elementary matrix operations: 10.8LinearProgrammingandtheSimplexMethod 429Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).•Locate the pivot element and save it. •Save the whole pivot column. •Replaceeachrow,exceptthepivotrow,bythatlinearcombinationofitself and the pivot row which makes its pivot-column entry zero. •Divide the pivot row by the negative of the pivot. •Replace the pivot element by the reciprocal of its saved value. •Replace the rest of the pivot column by its saved values divided by the saved pivot element. This is the sequence of operations actually performed by a linear programming routine, such as the one that we will presently give. You should now be able to solve almost any linear programmingproblem that starts in restricted normal form. The only special case that might stump you is if an entry in the constant column turns out to be zero at some stage, so that aleft-hand variable is zero at the same time as all the right-hand variables are zero. This is called a degenerate feasible vector . To proceed, you may need to exchange the degenerate left-hand variable for one of the right-hand variables, perhaps evenmaking several such exchanges. Writingthe GeneralProblem inRestricted NormalForm Here is a pleasant surprise. There exist a couple of clever tricks that render trivial the task of translating a general linear programming problem into restricted normal form! First, we need to get rid of the inequalities of the form (10.8.3)or (10.8.4),for example,the first three constraints in (10.8.7). We do this by addingto the problem so-called slack variables which, when their nonnegativity is required, convert the inequalities to equalities. We will denote slack variables as yi. There will be m1+m2of them. Once they are introduced, you treat them on an equal footing with the original variables xi; then, at the very end, you simply ignore them. For example, introducing slack variables leaves (10.8.6) unchanged but turns (10.8.7) into x1+2 x3+y1= 740 2x2−7x4+y2=0 x2−x3+2 x4−y3=1 2 x1+x2+x3+x4=9(10.8.15 ) (Notice how the sign of the coef ficient of the slack variable is determined by which sense of inequality it is replacing.) Second,we needto insure that there is a set of Mleft-handvectors, so that we can set up a starting tableau in restricted normal form. (In other words, we need to find a“feasible basic starting vector. ”) The trick is again to invent new variables! Thereare Mofthese,andtheyarecalled artificialvariables ;we denotethemby zi. Youputexactlyonearti ficial variableintoeachconstraintequationonthefollowing 430 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).model for the example (10.8.15): z1= 740 −x1−2x3−y1 z2=−2x2+7 x4−y2 z3=1 2−x2+x3−2x4+y3 z4=9−x1−x2−x3−x4(10.8.16 ) Our example is now in restricted normal form. Now you may object that (10.8.16) is not the same problem as (10.8.15) or (10.8.7)unless all the zi’s are zero . Right you are! There is some subtlety here! We must proceed to solve our problem in two phases. First phase: We replace our objective function (10.8.6) by a so-called auxiliary objective function z/prime≡− z1−z2−z3−z4=−(7491 2−2x1−4x2−2x3+4 x4−y1−y2+y3) (10.8.17 ) (where the last equality follows from using 10.8.16). We now perform the simplex method on the auxiliary objective function (10.8.17) with the constraints (10.8.16).Obviouslytheauxiliaryobjectivefunctionwillbemaximizedfornonnegative z i’si f all the zi’s are zero. We therefore expect the simplex method in this first phase to produce a set of left-hand variables drawn from the xi’s and yi’s only, with all the zi’s being right-handvariables. Aha! We then cross out the zi’s, leaving a problem involvingonly xi’sand yi’sinrestrictednormalform. Inotherwords,the firstphase producesaninitialfeasiblebasicvector. Secondphase: Solvetheproblemproduced by thefirst phase, using the original objective function, not the auxiliary. And what if the first phase doesn’tproduce zero values for all the zi’s? That signals that there is noinitial feasible basic vector, i.e., that the constraints given to us are inconsistent among themselves. Report that fact, and you are done. Hereishowtotranslateintotableauformattheinformationneededforboththe first and second phases of the overall method. As before, the underlying problem to be solved is as posed in equations (10.8.6) –(10.8.7). x1 x2 x3 x4 y1 y2 y3 z 0 1 1 3−1 20 0 0 z1 740 −1 0−2 0 −1 0 0 z2 0 0−2 0 7 0−1 0 z31 20−1 1 −2 0 0 1 z4 9−1−1−1 −1 0 0 0 z/prime−7491 22 4 2 −4 1 1−1 (10.8.18 ) This is not as daunting as it may, at first sight, appear. The table entries inside the box of double lines are no more than the coef ficients of the original problem (10.8.6)–(10.8.7) organized into a tabular form. In fact, these entries, along with 10.8LinearProgrammingandtheSimplexMethod 431Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).the values of N,M,m1,m2, and m3, are the only input that is needed by the simplex method routine below. The columns under the slack variables yisimply recordwhethereachofthe Mconstraintsisoftheform ≤,≥,or =;thisisredundant informationwith the values m1,m 2,m 3, as longas we are sure to enter the rows of thetableauinthecorrectrespectiveorder. Thecoef ficientsoftheauxiliaryobjective function (bottom row) are just the negatives of the column sums of the rows above, so these are easily calculated automatically. The output from a simplex routine will be (i) a flag telling whether a finite solution,nosolution,oranunboundedsolutionwasfound,and(ii)anupdatedtableau. The outputtableau that derivesfrom(10.8.18),givento two signi ficantfigures,is x1 y2 y3··· z 17.03 −.95 −.05 −1.05 ··· x2 3.33 −.35 −.15 .35 ··· x3 4.73 −.55 .05 −.45 ··· x4 .95 −.10 .10 .10 ··· y1 730 .55 .10 −.10 .90 ··· (10.8.19 ) A little counting of the xi’s and yi’s will convince you that there are M+1 rows (including the z-row) in both the input and the output tableaux, but that only N+1−m3columnsof the outputtableau(includingthe constantcolumn)contain any useful information, the other columns belonging to now-discarded arti ficial variables. In the output, the first numerical column contains the solution vector, alongwiththemaximumvalueoftheobjectivefunction. Whereaslackvariable( yi) appears on the left, the corresponding value is the amount by which its inequality is safely satis fied. Variables that are not left-hand variables in the output tableau have zero values. Slack variables with zero values represent constraints that are satisfied as equalities. RoutineImplementingthe SimplexMethod ThefollowingroutineisbasedalgorithmicallyontheimplementationofKuenzi, Tzschach, and Zehnder [4]. Aside from input values of M,N,m1,m2,m3, the principal input to the routine is a two-dimensionalarray acontainingthe portionof thetableau(10.8.18)thatis containedbetweenthedoublelines. Thisinputoccupies thefirstM+1rowsand N+1columnsof a. Note,however,thatreferenceismade internally to row M+2ofa(used for the auxiliary objective function, just as in 10.8.18). Therefore the physical dimensions of a, REAL a(MP,NP) (10.8.20 ) musthave NP≥N+1andMP≥M+2.You will suffer endless agonies if you fail to understand this simple point. Also do not neglect to order the rows of ain the same order as equations (10.8.1), (10.8.3), (10.8.4), and (10.8.5), that is, objective function, ≤-constraints, ≥-constraints, =-constraints. 432 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Onoutput,thetableau aisindexedbytworeturnedarraysofintegers. iposv(j) contains,for j=1 ...M,thenumber iwhoseoriginalvariable xiisnowrepresented byrow j+1ofa. Thesearethustheleft-handvariablesinthesolution. (The firstrow ofais of course the z-row.) A value i>Nindicates that the variable is a yirather thanan xi,xN+j≡yj. Likewise, izrov(j) contains,for j=1 ...N,thenumber i whose original variable xiis now a right-handvariable, representedby column j+1 ofa. These variablesareall zerointhesolution. Themeaningof i>Nis thesame asabove,exceptthat i>N +m1+m2denotesanarti ficialorslackvariablewhich was used only internally and should now be entirely ignored. Theflagicaseisreturnedaszeroifa finitesolutionisfound, +1iftheobjective function is unbounded, −1if no solution satis fies the given constraints. The routinetreats the case of degeneratefeasible vectors,so don ’tworryabout them. Youmayalso wishto admirethefact thattheroutinedoesnotrequirestoragefor the columns of the tableau (10.8.18) that are to the right of the double line; it keeps track of slack variables by more ef ficient bookkeeping. Please note that, as given, the routine is only “semi-sophisticated ”in its tests for convergence. While the routine properly implements tests for inequality with zero as tests against some small parameter EPS, it does not adjust this parameter to reflect the scale of the input data. This is adequate for many problems, where the input data do not differ from unity by too many orders of magnitude. If, however, you encounter endless cycling, then you should modify EPSin the routines simplx andsimp2. Permuting your variables can also help. Finally, consult [5]. SUBROUTINE simplx(a,m,n,mp,np,m1,m2,m3,icase,izrov,iposv) INTEGER icase,m,m1,m2,m3,mp,n,np,iposv(m),izrov(n),MMAX,NMAX REAL a(mp,np),EPSPARAMETER (MMAX=100,NMAX=100,EPS=1.e-6) C USES simp1,simp2,simp3 Simplex method for linear programming. Input parameters a,m,n,mp,np,m1,m2,a n d m3, and output parameters a,icase,izrov,a n d iposvare described above. Parameters: MMAXis the maximum number of constraints expected; NMAXis the maximum number of variables expected; EPSis the absolute precision, which should be adjusted to the scale of your variables. INTEGER i,ip,is,k,kh,kp,nl1,l1(NMAX),l3(MMAX)REAL bmax,q1 if(m.ne.m1+m2+m3)pause ’bad input constraint counts in simplx’ nl1=ndo 11k=1,n l1(k)=k Initialize index list of columns admissible for exchange. izrov(k)=k Initially make all variables right-hand. enddo 11 do12i=1,m if(a(i+1,1).lt.0.)pause ’bad input tableau in simplx’ Constants bimustbenon- negative. iposv(i)=n+i Initial left-hand variables. m1type constraints are represented by having their slackvari- able initially left-hand, with no artificial variable. m2type constraints have their slack variable initially left-hand, with a minus sign, and their artificial variable handled implic- itly during their first exchange. m3type constraints have their artificial variable initially left-hand. enddo 12 if(m2+m3.eq.0)goto 30 The origin is a feasible starting solution. Go to phase two. do13i=1,m2 Initializelistof m2constraints whoseslackvariableshavenever been exchanged out of the initial basis. l3(i)=1 enddo 13 do15k=1,n+1 Compute the auxiliary objective function. q1=0. do14i=m1+1,m 10.8LinearProgrammingandtheSimplexMethod 433Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).q1=q1+a(i+1,k) enddo 14 a(m+2,k)=-q1 enddo 15 10 call simp1(a,mp,np,m+1,l1,nl1,0,kp,bmax) Find max. coeff. of auxiliary objec- tive fn. if(bmax.le.EPS.and.a(m+2,1).lt.-EPS)then icase=-1 Auxiliary objective function is still negative and can’t be im- proved, hence no feasible solution exists. return else if(bmax.le.EPS.and.a(m+2,1).le.EPS)then Auxiliary objective function is zero and can’t be improved; we have a feasible starting vec- tor. Clean out the artificial variables corresponding to any remaining equality constraints bygoto 1’s and then move on to phase two by goto 30. do 16ip=m1+m2+1,m if(iposv(ip).eq.ip+n)then Found an artificial variable for an equality constraint. call simp1(a,mp,np,ip,l1,nl1,1,kp,bmax) if(bmax.gt.EPS)goto 1 Exchangewithcolumncorresponding tomax- imum pivot element in row. endif enddo 16 do18i=m1+1,m1+m2 Change sign of row for any m2constraints still present from the initial basis. if(l3(i-m1).eq.1)then do17k=1,n+1 a(i+1,k)=-a(i+1,k) enddo 17 endif enddo 18 goto 30 Go to phase two. endif call simp2(a,m,n,mp,np,ip,kp) Locate a pivot element (phase one). if(ip.eq.0)then Maximum of auxiliary objective function is unbounded, so no feasible solution ex- ists.icase=-1 return endif 1 call simp3(a,mp,np,m+1,n,ip,kp) Exchange a left- and a right-hand variable (phase one), then update lists. if(iposv(ip).ge.n+m1+m2+1)then Exchanged out an artificial variable for an equality constraint. Make sure it staysout by removing it from the l1list.do 19k=1,nl1 if(l1(k).eq.kp)goto 2 enddo 19 2 nl1=nl1-1 do21is=k,nl1 l1(is)=l1(is+1) enddo 21 else kh=iposv(ip)-m1-nif(kh.ge.1)then Exchanged out an m2type constraint. if(l3(kh).ne.0)then If it’s the first time, correct the pivot col- umn for the minus sign and the implicitartificial variable.l3(kh)=0 a(m+2,kp+1)=a(m+2,kp+1)+1. do 22i=1,m+2 a(i,kp+1)=-a(i,kp+1) enddo 22 endif endif endifis=izrov(kp) Update lists of left- and right-hand variables. izrov(kp)=iposv(ip) iposv(ip)=isgoto 10 Still in phase one, go backto 10. End of phase one code for finding an initial feasible solution. Now, in phase two, optimize it. 30 call simp1(a,mp,np,0,l1,nl1,0,kp,bmax) Test the z-row for doneness. if(bmax.le.EPS)then Done. Solution found. Return with the good news. icase=0return endif 434 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).call simp2(a,m,n,mp,np,ip,kp) Locate a pivot element (phase two). if(ip.eq.0)then Objective function is unbounded. Report and return. icase=1 return endif call simp3(a,mp,np,m,n,ip,kp) Exchange a left- and a right-hand variable (phase two), is=izrov(kp) update lists of left- and right-hand variables, izrov(kp)=iposv(ip)iposv(ip)=is goto 30 and return for another iteration. END The precedingroutine makes use of the following utility subroutines. SUBROUTINE simp1(a,mp,np,mm,ll,nll,iabf,kp,bmax) INTEGER iabf,kp,mm,mp,nll,np,ll(np) REAL bmax,a(mp,np) Determines the maximum of those elements whose index is contained in the supplied list ll, either with or without taking the absolute value, as flagged by iabf. INTEGER k REAL testif(nll.le.0)then No eligible columns. bmax=0. else kp=ll(1)bmax=a(mm+1,kp+1) do 11k=2,nll if(iabf.eq.0)then test=a(mm+1,ll(k)+1)-bmax else test=abs(a(mm+1,ll(k)+1))-abs(bmax) endifif(test.gt.0.)then bmax=a(mm+1,ll(k)+1) kp=ll(k) endif enddo 11 endifreturnEND SUBROUTINE simp2(a,m,n,mp,np,ip,kp) INTEGER ip,kp,m,mp,n,npREAL a(mp,np),EPS PARAMETER (EPS=1.e-6) Locate a pivot element, taking degeneracy into account. INTEGER i,k REAL q,q0,q1,qp ip=0do 11i=1,m if(a(i+1,kp+1).lt.-EPS)goto 1 enddo 11 return No possible pivots. Return with message. 1 q1=-a(i+1,1)/a(i+1,kp+1) ip=i do13i=ip+1,m if(a(i+1,kp+1).lt.-EPS)then q=-a(i+1,1)/a(i+1,kp+1) if(q.lt.q1)then ip=iq1=q else if (q.eq.q1) then We have a degeneracy. 10.8LinearProgrammingandtheSimplexMethod 435Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do12k=1,n qp=-a(ip+1,k+1)/a(ip+1,kp+1) q0=-a(i+1,k+1)/a(i+1,kp+1) if(q0.ne.qp)goto 2 enddo 12 2 if(q0.lt.qp)ip=i endif endif enddo 13 returnEND SUBROUTINE simp3(a,mp,np,i1,k1,ip,kp) INTEGER i1,ip,k1,kp,mp,np REAL a(mp,np) Matrix operations to exchange a left-hand and right-hand variable (see text). INTEGER ii,kk REAL piv piv=1./a(ip+1,kp+1)do 12ii=1,i1+1 if(ii-1.ne.ip)then a(ii,kp+1)=a(ii,kp+1)*pivdo 11kk=1,k1+1 if(kk-1.ne.kp)then a(ii,kk)=a(ii,kk)-a(ip+1,kk)*a(ii,kp+1) endif enddo 11 endif enddo 12 do13kk=1,k1+1 if(kk-1.ne.kp)a(ip+1,kk)=-a(ip+1,kk)*piv enddo 13 a(ip+1,kp+1)=piv return END OtherTopics Briefly Mentioned Every linear programming problem in normal form with Nvariables and M constraints has a corresponding dualproblem with Mvariables and Nconstraints. The tableau of the dual problem is, in essence, the transpose of the tableau of the original (sometimes called primal) problem. It is possible to go from a solution of the dual to a solution of the primal. This can occasionally be computationally useful, but generally it is no big deal. Therevised simplex method is exactly equivalent to the simplex method in its choiceofwhichleft-handandright-handvariablesareexchanged. Itscomputational effort is not signi ficantly less than that of the simplex method. It does differ in the organization of its storage, requiring only a matrix of size M×M, rather than M×N, in its intermediate stages. If you have a lot of constraints, and memory size is one of them, then you should look into it. Theprimal-dual algorithm and thecomposite simplex algorithm are two dif- ferent methods for avoiding the two phases of the usual simplex method: Progress is made simultaneously towards finding a feasible solution and finding an optimal solution. There seems to be no clearcut evidence that these methods are superior 436 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).to the usual method by any factor substantially larger than the “tender-loving-care factor”(which re flects the programming effort of the proponents). Problemswheretheobjectivefunctionand/oroneormoreoftheconstraintsare replacedbyexpressionsnonlinearinthevariablesarecalled nonlinearprogramming problems. Theliteratureonsuchproblemsisvast,butoutsideourscope. Thespecial case of quadraticexpressions is called quadraticprogramming . Optimization prob- lemswherethevariablestakeononlyintegervaluesarecalled integerprogramming problems, a special case of discrete optimization generally. The next section looks at a particular kind of discrete optimization problem. CITED REFERENCES AND FURTHER READING: Bland, R.G. 1981, Scientific American , vol. 244 (June), pp. 126–144. [1] Dantzig, G.B. 1963, Linear Programming and Extensions (Princeton, NJ: Princeton University Press). [2] Kolata, G. 1982, Science, vol. 217, p. 39. [3] Gill,P.E., Murray, W.,andWright, M.H. 1991, NumericalLinearAlgebraandOptimization , vol. 1 (Redwood City, CA: Addison-Wesley), Chapters 7–8. Cooper,L.,andSteinberg,D.1970, IntroductiontoMethodsofOptimization (Philadelphia:Saun- ders). Gass, S.T. 1969, Linear Programming , 3rd ed. (New York: McGraw-Hill). Murty, K.G. 1976, Linear and Combinatorial Programming (New York: Wiley). Land,A.H., andPowell,S.1973, FortranCodesforMathematicalProgramming (London:Wiley- Interscience). Kuenzi, H.P., Tzschach, H.G., and Zehnder, C.A. 1971, Numerical Methods of Mathematical Optimization (New York: Academic Press). [4] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §4.10. Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [5] 10.9 Simulated Annealing Methods Themethodofsimulatedannealing [1,2]isa techniquethathasattractedsignif- icant attention as suitable for optimization problems of large scale, especially ones wherea desiredglobalextremumis hiddenamongmany,poorer,localextrema. For practicalpurposes,simulatedannealinghaseffectively “solved”thefamous traveling salesman problem offinding the shortest cyclical itinerary for a traveling salesman who must visit each of Ncities in turn. (Other practical methods have also been found.) Themethodhasalsobeenusedsuccessfullyfordesigningcomplexintegrated circuits: The arrangement of several hundred thousand circuit elements on a tiny siliconsubstrateisoptimizedsoas tominimizeinterferenceamongtheirconnectingwires [3,4]. Surprisingly,the implementationof thealgorithmis relativelysimple. Notice that the two applications cited are both examples of combinatorial minimization . Thereisanobjectivefunctiontobeminimized,asusual;butthespace over which that function is de fined is not simply the N-dimensional space of N