c10-8 (Numerical Recipes in C), страница 4

In the output, the first numerical column contains the solution vector,along with the maximum value of the objective function. Where a slack variable (yi )appears on the left, the corresponding value is the amount by which its inequalityis safely satisfied. Variables that are not left-hand variables in the output tableauhave zero values. Slack variables with zero values represent constraints that aresatisfied as equalities.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).z10.8 Linear Programming and the Simplex Method439Routine Implementing the Simplex Methoda[i][k],i = 1 . . . m+2, k = 1 .

. . n+1(10.8.20)You will suffer endless agonies if you fail to understand this simple point. Also donot neglect to order the rows of a in 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.On output, the tableau a is indexed by two returned arrays of integers.

iposv[j]contains, for j= 1 . . . M , the number i whose original variable xi is now representedby row j+1 of a. These are thus the left-hand variables in the solution. (The first rowof a is of course the z-row.) A value i > N indicates that the variable is a yi ratherthan an xi , xN+j ≡ yj . Likewise, izrov[j] contains, for j= 1 . . .

N , the number iwhose original variable xi is now a right-hand variable, represented by column j+1of a. These variables are all zero in the solution. The meaning of i > N is the sameas above, except that i > N + m1 + m2 denotes an artificial or slack variable whichwas used only internally and should now be entirely ignored.The flag icase is set to zero if a finite solution is found, +1 if the objectivefunction is unbounded, −1 if no solution satisfies the given constraints.The routine treats the case of degenerate feasible vectors, so don’t worry aboutthem.

You may also wish to admire the fact that the routine does not require storagefor the columns of the tableau (10.8.18) that are to the right of the double line; itkeeps track of slack variables by more efficient bookkeeping.Please note that, as given, the routine is only “semi-sophisticated” in its testsfor convergence. While the routine properly implements tests for inequality withzero as tests against some small parameter EPS, it does not adjust this parameter toreflect the scale of the input data.

This is adequate for many problems, where theinput data do not differ from unity by too many orders of magnitude. If, however,you encounter endless cycling, then you should modify EPS in the routines simplxand simp2. Permuting your variables can also help.

Finally, consult [5].#include "nrutil.h"#define EPS 1.0e-6Here EPS is the absolute precision, which should be adjusted to the scale of your variables.#define FREEALL free_ivector(l3,1,m);free_ivector(l1,1,n+1);void simplx(float **a, int m, int n, int m1, int m2, int m3, int *icase,int izrov[], int iposv[])Simplex method for linear programming. Input parameters a, m, n, mp, np, m1, m2, and m3,and output parameters a, icase, izrov, and iposv are described above.{void simp1(float **a, int mm, int ll[], int nll, int iabf, int *kp,Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The following routine is based algorithmically on the implementation of Kuenzi,Tzschach, and Zehnder [4].

Aside from input values of M , N , m1 , m2 , m3 , theprincipal input to the routine is a two-dimensional array a containing the portion ofthe tableau (10.8.18) that is contained between the double lines. This input occupiesthe M + 1 rows and N + 1 columns of a[1..m+1][1..n+1]. Note, however, thatreference is made internally to row M + 2 of a (used for the auxiliary objectivefunction, just as in 10.8.18). Therefore the variable declared as float **a, mustpoint to allocated memory allowing references in the subrange440Chapter 10.Minimization or Maximization of Functionsfloat *bmax);void simp2(float **a, int n, int l2[], int nl2, int *ip, int kp, float *q1);void simp3(float **a, int i1, int k1, int ip, int kp);int i,ip,is,k,kh,kp,nl1;int *l1,*l3;float q1,bmax;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).if (m != (m1+m2+m3)) nrerror("Bad input constraint counts in simplx");l1=ivector(1,n+1);l3=ivector(1,m);nl1=n;for (k=1;k<=n;k++) l1[k]=izrov[k]=k;Initialize index list of columns admissible for exchange, and make all variables initiallyright-hand.for (i=1;i<=m;i++) {if (a[i+1][1] < 0.0) nrerror("Bad input tableau in simplx");Constants bi must be nonnegative.iposv[i]=n+i;Initial left-hand variables. m1 type constraints are represented by having their slackvariable initially left-hand, with no artificial variable. m2 type constraints have theirslack variable initially left-hand, with a minus sign, and their artificial variable handledimplicitly during their first exchange.

m3 type constraints have their artificial variableinitially left-hand.}if (m2+m3) {Origin is not a feasible starting sofor (i=1;i<=m2;i++) l3[i]=1;lution: we must do phase one.Initialize list of m2 constraints whose slack variables have never been exchanged outof the initial basis.for (k=1;k<=(n+1);k++) {Compute the auxiliary objective funcq1=0.0;tion.for (i=m1+1;i<=m;i++) q1 += a[i+1][k];a[m+2][k] = -q1;}for (;;) {simp1(a,m+1,l1,nl1,0,&kp,&bmax);Find max.

coeff. of auxiliary objecif (bmax <= EPS && a[m+2][1] < -EPS) {tive fn.*icase = -1;Auxiliary objective function is still negative and can’t be improved, hence nofeasible solution exists.FREEALL return;} else if (bmax <= EPS && a[m+2][1] <= EPS) {Auxiliary objective function is zero and can’t be improved; we have a feasiblestarting vector. Clean out the artificial variables corresponding to any remainingequality constraints by goto one and then move on to phase two.for (ip=m1+m2+1;ip<=m;ip++) {if (iposv[ip] == (ip+n)) {Found an artificial variable for ansimp1(a,ip,l1,nl1,1,&kp,&bmax);equality constraint.if (bmax > EPS)Exchange with column correspondgoto one;ing to maximum pivot element}in row.}for (i=m1+1;i<=m1+m2;i++)Change sign of row for any m2 conif (l3[i-m1] == 1)straints still present from the inifor (k=1;k<=n+1;k++)tial basis.a[i+1][k] = -a[i+1][k];break;Go to phase two.}simp2(a,m,n,&ip,kp);Locate a pivot element (phase one).if (ip == 0) {Maximum of auxiliary objective func*icase = -1;tion is unbounded, so no feasiFREEALL return;ble solution exists.}one:simp3(a,m+1,n,ip,kp);Exchange a left- and a right-hand variable (phase one), then update lists.10.8 Linear Programming and the Simplex Method441}The preceding routine makes use of the following utility functions.#include <math.h>void simp1(float **a, int mm, int ll[], int nll, int iabf, int *kp,float *bmax)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.{int k;float test;if (nll <= 0)No eligible columns.*bmax=0.0;else {*kp=ll[1];*bmax=a[mm+1][*kp+1];for (k=2;k<=nll;k++) {if (iabf == 0)test=a[mm+1][ll[k]+1]-(*bmax);Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.