A C CODE FOR SOLVING THE GENERALIZED ASSIGNMENT PROBLEM Nejat Karabakal Department of Industrial & Operations Engineering University of Michigan Ann Arbor, MI 48109-2117 Technical Report 92-32 May 1992

kblgap A C Code for Solving the Generalized Assignment Problem* Nejat Karabakal Department of Industrial and Operations Engineering The University of Michigan, Ann Arbor, Michigan 48109 Version 1.0 of May 1992 Contents 1 Introduction 2 User Manual 2.1 Compiling and Running 2.2 Setting Parameters........ 2.3 An Example...................... 1 1 1 2 4 3 Code Documentation 3.1 Basic Definitions and Data Structure 3.2 Input Functions...... readin()...... setparams (). 3.3 Upper Bounding Functions..... MAM()........ fbknap()...... subgradient()......... knapdp()............. 3.4 Lower Bounding Function (Heuristic) lowerbound()........ 3.5 Fathoming Functions......... fathomroot().......... fathomnode()......... 3.6 Branch-and-Bound Function..... BAB()....... 3.7 Priority Queue Handling Functions. putdangling().... getdangling().... removeinferiors...... 3.8 Main Function and Utilities.... printsol()..... checkincumbent ()...... main()....... 7.............................. 7 7 8.............................. 8.............................. 8............................... 1..... 10.............................. 1 1............................. 1 2................................................... 15.............................................. 2 1 ~............................................. 21..............................~................. 23................................................. 2 3.............................................. 24..................................................26............................................... 2 6................................................28............................................ 28................................................28.............................. 2 8...................................... 2 8................................................29................................................29................................................30................................................30 *This research was supported by National Science Foundation grants PYI-G-DMC-8352346 and DMM-9018515 to The University of Michigan.

1 Introduction This report is a documentation and user manual for the C code, kblgap, developed for solving the following formulation of the generalized assignment problem (GAP): m n maximize cij xij i=1 j=l n subject to aijx <i bi i = 1,.., m j=1 m = 1 = 1,...,n i=1 xij {0,1} all i,j The code implements a branch-and-bound algorithm featuring an efficient Lagrangian relaxation methodology for bounding. Computational tests show that kblgap is, on the average, 100 times faster than the best branch-and-bound code (Martello-Toth code) available for the GAP to date. In addition, it is capable of solving larger problems (up to 500 variables) optimally in reasonable times. For more information, see Karabakal, Bean, and Lohmann (1992). Also, see Martello and Toth (1990) for computational comparisions of the Martello and Toth algorithm with other existing branch-and-bound algorithms, making it possible to compare kblgap with many other codes indirectly. 2 User Manual 2.1 Compiling and Running The code is written with portability considerations in mind. The author has successfully compiled it under MS-DOS using Turbo C compiler, tcc kblgap.c and under Unix, cc kblgap.c -okblgap -0 to produce the executable program. The executable kblgap is invoked by the command kblgap (input-file-name) If (input-file-name) is omitted, kblgap prompts the user to enter an input file name. The output goes to the standard output (default is the screen). If a file output is desired, output redirection can be used (under Unix and MS-DOS), e.g., kblgap (input-file-name) > output.txt The input file format is given below. All numbers are integers and seperated by at least one space. 1

m Cll all a1 C21 a21 Cml aml bl n C12 a12 C22 a22 Cm2 am2 b2... Cln...aln... C2n...a2n... Cmn... amn... bm 2.2 Setting Parameters The maximum problem size is given by three parameters that are used for dimensioning arrays. MAXM MAXN MAXB maximum number of knapsack constraints, maximum number of multiple-choice constraints, largest knapsack RHS capacity. They are specified at the beginning with #define statements (see section 3.1). If these limits are violated, kblgap terminates with a descriptive error message. These three parameters are "static" in the sense that whenever any of them are modified, the code must be recompiled. In contrast, the following parameters, used for fine-tuning the algorithmic performance and for controlling the detail of output, are "dynamic" in the sense that they can be specified in a seperate configuration file called kblgap. cfg. INFINITY ROOTMAMITLIM MAMITLIM ROOTSUBITLIM SUBITLIM SUBPATIENCE MAXBRANCH PRINTSTEPS PRINTINPUT PRINTOPTIM ZEROTOL a big number, maximum multiplier adjustment method (MAM) iterations made at the root node of the branch-and-bound tree, maximum MAM iterations made at any other node, maximum subgradient iterations made at the root node, maximum subgradient iterations allowed at any other node, number of iterations the subgradient algorithm permits the upper bound to diverge before halving the scale parameter u (see subgradient() function description for more information), maximum number of nodes that can be generated in the branch-and-bound tree, set to 1 if branching steps are to be printed, 0 otherwise, set to 1 if input data are to be printed, 0 otherwise, set to 1 if the final optimum solution is to be printed, 0 otherwise, zero tolerance. Each time kblgap is run, it first checks the existence of kblgap. cfg. If it exists, values of the above parameters are read in from this file; otherwise, the default values assigned at the beginning of the code take effect (see section 3.1). Each line of the configuration file consists of two quantities, the value of the parameter followed by a descriptive string. Parameters must be supplied in the above order. A sample kblgap. cfg file is provided below. 2

32000 INFINITY 10 ROOTMAMITLIM 5 MAMITLIM 25 ROOTSUBITLIM 10 SUBITLIM 4 SUBPATIENCE 5000 MAXBRANCH 0 PRINTSTEPS 1 PRINTINPUT 1 PRINTOPTIM 0.001 ZEROTOL Iteration limits for upper bounding procedures, MAM() and subgradient(), can be specified zero to prevent invoking a particular procedure. The parameter, MAXBRANCH, controls the accuracy of the solution. If an optimum solution is sought, a big integer is assigned to MAXBRANCH. On the other hand, a fairly good feasible solution with an error bound can be obtained by keeping this parameter low. In the latter case, the tightest upper bound available at the time of termination is used for computing the maximum error of the best feasible solution. A line of output with the following format is printed: Best solution found = (1), Upper bound = (2) (Error bound = (3)%) where 1. Objective value of the best feasible solution, 2. Tightest upper bound on optimum objective, 3. (2)1) x 100. (1) See section 2.3 for an example. If the parameter PRINTSTEPS is turned on (set to one), a line with the following format is printed for every node generated in the branch-and-bound tree. Node (1) From (2): (3) [(4)] ==> [(5)] (6) where 1. node's serial number, 2. parent node, 3. variable fixed at this node, 4. [current incumbent, upper bound of the subproblem at the parent node] pair, 5. [potentially better incumbent, upper bound of the subproblem at the current node] pair, 6. whether the node is fathomed or dangles after bounding. 3

Node 0 corresponds to the root node and its output line has a slightly different interpretation. Node 0: [(1)] ==> [(2)] (3) where 1. the interval [abslb, absub], where abslb = Ej min{cij} and absub = Ej max{cij} These quantities represent absolute lower and upper bounds, respectively, on the optimum objective of the GAP. 2. [best feasible solution, upper bound] pair for the original problem, 3. whether the root is fathomed or dangles after bounding (if it is fathomed, the best feasible solution is optimal; if it dangles, there is a duality gap). If the number of nodes generated is expected to be high, PRINTSTEPS switch should be turned off to prevent excessive output. See section 2.3 for examples. 2.3 An Example Consider a GAP with 5 knapsack and 10 multiple-choice constraints. Given below are an input data file 5bylO. cl set up for this problem and the code's output. 19 14 23 22 (cij) = 18 19 21 23 12 23 21 19 10 24 13 11 22 16 18 21 19 12 13 10 18 19 21 25 17 10 22 14 25 22 23 25 10 22 16 13 19 12 21 17 22 22 18 11 24 24 24 18 5 18 (aij) = 15 11 22 24 19 5 5 8 24 11 13 22 15 14 24 19 12 13 20 9 11 23 7 20 14 18 12 22 15 22 14 24 15 18 / 26 19 21 9 24 12 8 21 (bi) = 24 16. 12 23 28 18 5 16 / 22/ 5 10 19 14 21 24 18 5 23 22 18 5 18 11 18 19 21 15 11 15 21 23 23 22 24 19 12 23 21 19 5 20 26 24 24 28 19 8 21 13 25 14 25 12 17 9 22 10 24 13 11 22 16 24 23 7 24 15 18 19 12 13 10 18 19 22 20 14 19 21 9 17 10 22 14 25 22 24 18 12 12 8 21 10 22 16 13 19 12 13 22 15 16 12 23 22 22 18 11 24 24 11 22 14 18 5 16 4

-— PARAMETERS SET Root iter limit Node iter limit Patience MAM SUB 5 20 5 10 4 PRINT FLAGS Input YES Branches YES OptimSol YES GENERAL Infinity Zero Tol Max Branch 32000 0.0010 10000 Program Limits: Max Rows 5, Max Columns 20, Max Knapsack RHS 200. PROBLEM DATA: S rows by 10 columns (Input file: Sby10.cl) Objective 19 14 21 23 22 18 18 19 21 21 23 23 12 23 21 coefficients: 19 10 24 13 11 22 16 21 19 12 13 10 18 19 25 17 10 22 14 25 22 25 10 22 16 13 19 12 17 22 22 18 11 24 24 Knapsack constraints: 24 18 5 8 24 23 7 24 15 18 <= 26 5 18 11 13 22 20 14 19 21 9 <= 24 15 11 15 14 24 18 12 12 8 21 <= 24 22 24 19 12 13 22 15 16 12 23 <= 28 19 5 20 9 11 22 14 18 5 16 <= 22 Node 0: [134, 225] ==> [201,215] Node Node Node Node Node Node Node Node Node Node 1 From 0: 2 From 0: 3 From 0: 4 From 0: 5 From 0: 6 From 5: 7 From 5: 8 From 5: 9 From 5: 10 From 5 x(1,9) x(2,9) x(3,9) x(4,9) x(5,9) x(1,2) x(2,2) x(3,2) x(4,2) [201,215] [201,215] [201,215] [201,215] [201,215] [209,212] [209,212] [209,212] [209,212] DANGLING [201,204] [201,200] [201,208] [201,205] [209,212] [209,204] [209,209] > [209,208] [209,206] DANGLING FATHOM DANGLING DANGLING DANGLING FATHOM FATHOM FATHOM FATHOM FATHOM: x(5,2) [209,212] ==> [209,209] Optimum = 209 No of branches = 10 cx = 209 1234567890 1)..11...... 2)1.......1 3)......11. 4)......... 5).1..1...1. USED 13 14 24 22 21 RHS 26 24 24 28 22 5

To illustrate the pieces of information printed for each branch-and-bound node, consider node 9. It is generated by branching from node 5 by fixing x42. Hence, the subproblem at node nine has two fixed variables, x42 and x59. Upper bounding procedures applied to this subproblem conclude that its optimum cannot exceed 206. Therefore, node 9 is fathomed. If the parameter MAXBRANCH is set to 0, the output would look like the following: -— PARAMETERS SET --- —------------------------------- MAM SUB PRINT FLAGS GENERAL Root iter limit 10 25 Input NO Infinity 32000 Node iter limit 5 10 Branches YES Zero Tol 0.0010 Patience 4 OptimSol YES Max Branch 0 Program Limits: Max Rows 5, Max Columns 20, Max Knapsack RHS 200. PROBLEM DATA: 5 rows by 10 columns (Input file: SbylO.cl) Node 0: [134, 225] ==> [201,215] DANGLING Branch limit exceeded. Best solution found = 201, Upper bound = 215 (Error bound = 6.97%) No of branches = 0 cx = 201 1234567890 USED RHS 1)..11..1.. 20 26 2)1........1 14 24 3).......11. 20 24 4).....1.... 22 28 5).1..1..... 16 22 6

3 Code Documentation 3.1 Basic Definitions and Data Structure #include <stdio.h> /* Maximum problem size, modify to solve bigger problems */ #define MAXM 5 /* max no. of knapsack constraints */ #define MAXN 20 /* max no. of multiple-choice constraints */ #define MAXB 200 /* biggest knapsack capacity */ /* Mnemonic definitions */ #define YES 1 #define NO 0 typedef int bool; /* Structure of a branch-and-bound tree node */ typedef struct NODE { int nodeno, /* node's serial no */ bound, /* upper bound of the node's partial problem */ ifix, jfix, /* indices of the fixed variable */ lambda[MAXN+1]; /* multipliers used in upper bounding */ struct NODE *next, *parent; /* pointers (see below) */ } NODE, *NODEPTR; The branch-and-bound tree consists of nodes as defined above. Each outstanding node is a member of either or both of the two linked lists maintained for efficient operations: 1. Path from node to the root: Each node points to its parent through the *parent link, except for the root node. The root's *parent is set to NULL indicating the end of a list traversal. By following the path from a particular node to the root, the algorithm determines which variables are fixed in the partial problem associated with this node. Another use of the *parent pointer is to initialize the multiplier vector with the best multipliers computed at the parent node. 2. Dangling nodes: This is a priority queue consisting only of dangling nodes, i.e., those without children. It is ordered by decreasing bound values. Nodes are connected via *next links. The start of the queue is pointed by firstdangling, so that firstdangling->bound always gives the highest upper bound among all yet-to-be-solved subproblems. When firstdangling = NULL, the branch-and-bound algorithm is terminated because this condition implies there exists no upper bound greater than the incumbent (see BAB() function for more information). /* Set defaults for parameters */ /* (They can be modified by the use of configuration file kblgap.cfg) */ 7

int INFINITY = 32000, ROOTMAMITLIM = 5, MAMITLIM = 5, ROOTSUBITLIM = 20, SUBITLIM = 10, SUBPATIENCE = 4, MAXBRANCH = 10000, PRINTSTEPS = YES, PRINTINPUT = YES, PRINTOPTIM = YES; float ZEROTOL = 0.001; /* Global variables */ int m, n, /* actual no of knapsack and MC constraints */ c[MAXM+1] [MAXN+1], /* objective values */ a[MAXM+l] [MAXN+1], /* knapsack weights */ b[MAXM+1], /* knapsack capacities, or RHS values */ abslb, /* sum of min objective values over every MC set */ nbranch = 0, /* branch counter */ incumbent, xinc[MAXN+l1]; /* incumbent value and solution */ Feasible solutions handled in various routines are kept in n-dimensional vectors, called positions. Let x is a feasible solution to the GAP and pos is the corresponding position vector. Then, for each j, pos [j] = i iff xij = 1. Incumbent solution xinc is a position vector. NODEPTR firstdangling = NULL; /* initially, there is no dangling nodes */ 3.2 Input Functions void readin(inpname) char inpname[]; { -FILE *inp; int i, j; if (inpname[O] == '*') { printf("Enter input file: "); } scanf("%s", inpname); if ((inp = fopen(inpname, "r")) == NULL) { printf( "ERROR: Cannot open input file: %s\n", inpname); exit(1); } 8

fscanf inp, "%d~d\n"I, &m, &n); if (m > MAXM) { printf("ERROR: Row limit exceeded.") printf (" (Data = %d, Limit %d) \n", m, MAXM); exit (1) } if (n > MAXN){ printf("IERROR: Column limit exceeded.") printf (" (Data = %d, Limit %d)\n", n, MAXN); exit (1) } for (i1l; i <= m; i++){ for (j1; j<=n; j++){ fscanf(inp, "/7.d"I, &c[i][j]); if (c[i][j] < 0) { printf("ERROR: Nonnegative objective value expected."); printf(" (Data: cD'.d,%d] = %d)\n", i, j, c[iJ[jJ); exit (1) fscanf(inp, "\n"); for (j1; j<=n; j++){ fscanf(inp, "U.", &a[i][j]); if (a[i][j] <= 0) { printf("ERROR: Positive knapsack weight expected."); printf(" (Data: a[Yd,Yd] = '/d)\n", i, j, a[i][j]); exit(1); fscanf~inp, "\n"); for (i1l; i <= m; i++) { fscanf(inp, "7.", &b[i]); if (b~il > MAXB) { printf("ERROR: Knapsack RHS too big."); printf(" (Data: b['hdl = '/d, Limit 'hd)\n", i, b~i], MAXB); exit (1) fscanf(inp, "\n"); fclose(inp); 9

printf("PROBLEM DATA: Xd rows by %d columns", m, n); printf("1 (Input file: /.s)\n\n", inpname); if (PRINTINPUT) { printf("Objective coefficients:\n"); f or (i1; i <= m; i++) { for (j1; j<=n; j++) printf("1%3d", c~ilj]j); printf("\n"); printf("\nKnapsack constraints:\n"); for (i1; i <= m; i++){ for (j1; j<n; j++) printf("%3d", a~iJ[j]); printf (" <= Xd\n", b [i]); } printf("\n"); } } /* end of readin * void setparams() char iflag[3], bflag[3], of lag[3]; FILE *cfg; mnt i, j; char dummy[20]; if ((cfg = fopen("1kblgap.cfg", "r")) = NULL){ printf( "Cannot open kblgap-cfg, defaults are in effect.\n"); } else { fscanf (cfg, fscanf (cfg, fscanf (cfg, fscanf (cfg, fscanf (cfg, fscanf (cfg, fscanf (cfg, fscanf (cfg, fscanf (cfg, fscanf(cfg, fscanf (cfg, 1 1 44 d "Yed "Yed "Yed 117.d 117ild 11 7od 11711d "Yod "Yod Ilyof Y.s\nll I Y.s\nll I Y.'s \ n I II Y.s\nll I Y.s\nllI Y.s\nll.9 Y.s\nllI Y.s\nll Y.s\nll ) Y.s\nll Y.s\nll j &INFINITY, dummy); &ROOTMAMITLIM, d3ummy);_ &MAMITLIM, dummy); &ROOTSUBITLIM, dummy); &SUBITLIM, dummy); &SUBPATIENCE, dummy); &MAXBRANCH, dummy); &PRINTSTEPS, dummy); &PRINTINPIJT, dummy); &PRINTOPTIM, dummy); &ZEROTOL, dummy); f close (cf g); } 10

if ((ROOTMAMITLIM <= 0 && ROOTSUBITLIM <= 0) II (MAMITLIM <= 0 && SUBITLIM <= 0)){ printf("IERROR: One of iter limits for upper bounding must be positive.\n"); printf("I (Data: ROOTMAMITLIM = Xd, ROOTSUBITLIM = dn IROOTMAMITLIM, ROOTSUBITLIM); printf(" MAMITLIM = %d, SUBITLIM = %d)\n"., MAMITLIM, SUBITLIM); exit (1) if (PRINTINPUT) strcpy(iflag,"YES"); else strcpy(iflag,"NO") if (PRINTSTEPS) strcpy(bflag,"YES"l); else strcpy(bflag,"NO") if (PRINTOPTIM) strcpy(oflag,"YES"); else strcpy(oflag,"NO") printf("\n --- PARAMETERS SET"); printf(" --- —----------------------— \n"); printf ( I I MAM SUB PRINT FLAGS printf ( 11Root iter limit 7%3d 7/.3d Input. ROOTMAMITLIM, ROOTSUBITLIM, if lag, INFINITY); printf ( "Node iter limit 7*3d /.3d Branches %s MAMITLIM, SUBITLIM, bf lag, ZEROTOL); printf ( Patience 7%3d Optimsol %s SUBPATIENCE, of lag, MAXBRANCH); printf("~ --- - - -- - - -- - -- - - -- - - printf("-\n"); --- —----------— W printf ( "Program Limits: Max Rows Yd, Max Columns 7.d, Max MAXM, MAXN, MAXB); printf("~ --- - - -- - - -- - -- - - -- - - printf("-\n"); --- —----------— W GENERAL\n"); Infinity %d\n", Zero Tol %6.4f\n", Max Branch %d\n", Knapsack RHS Yd.dn", }/* end of setparams * 3.3 Upper Bounding ]Functions Two procedures are run in serial to obtain ani upper bound on the optimum objective of the GAP or any of its subproblems obtained by fixing some of the variables. First, a multiplier adjustment method (MAM) is called, followed by a, subgradient algorithm. Both procedures solve a Lagrangian dual of the GAP for upper bounding. Associated Lagrangian relaxation is obtained by dualizing multiple-choice constraints. MAMO and subgradiento share a commoni input list: *Forj'= 1,...,n fI if x - =1I and xe,- = 0,f4 $i jfixed~j]= ~0 if xi 1i ---,,.,x,- are all free Let S ={Il'Ijf ixed[j]#:~0}. 11I

* Sum of objective values of fixed variables obj fixed = Ck 3, k = if ixedlij] 3JES * For i = 1...,Im bbar[i] = b- S ak3, k =jfixedlj] {a ESlk=i'} * For j=1.., Jlambda[j] =initial multiplier associated with multiple-choice constraint j. The outputs are the vector of final multipliers and the function value itself that gives the tightest upper bound obtained. The following code for the MAM, described in detail in Karabakal, Bean, and Lohmann (1992), is based on the idea of steepest descent improvements per iteration in solving the Lagrangian dual. Computationally, it outperforms its predecessors suggested by Fisher, Jaikumar, and Van Wassenhove (1986) and Guignard and Rosenwein (1989). jint MAM(jf ixed, objf ixed, bbar, lambda) int j fixed[] obj fixed, bbar[]l, lamnbda[]l int iter, i, j, k, rhs, istar, jstar, lastjstar, firsttime, row, p[MAXN+1], x[MAXM+1] [MAXN+11, pos [MAXN+1], aknap[MAXN+1], xknap[MAXN+1], seq[MAXN+1], dual, cknap[MAXN+1], vknap, slack[MAXN+1], delta[MAXM+1] [MAXN+1], deltal, delta2, imp, maximp, steplength, maxstep, bound; bool optcomp; void f bknapO checkincumbent; 1* initialize'* last'star = 0; firsttime = YES; iter = 0; f* initialize dual * for (j1, dual=0; j<=n; j++) if (jfixed~j]==0) dual += lambda~j]; while (1){ /* initialize column totals * for (j1; j<=n; j++) p[jl = 0; 12

/* solve m forward-backward knapsacks * for (i1; i<=m; i++) { k = 0; /* k counts no of elements in knapsack * for (j1; j<=n; j++) if (jfixed[jJ = 0) if (a[iII[jJ <= bbar[i]) seq[++k] = cknap [k] = c Ei] [jIi - lambda [jJ Iaknap [kJ a a[i] rjI] for (j1; j<=n; j++){ x[i][j]= 0; delta[i][j] = INFINITY; } if (k > 0){ rhs = bbar[i]; fbknap(k, rhs, cknap, aknap, xknap, &vknap, slack); if (firsttime) dual += vknap;P for (j=1; j<=k; j++) { if (xknap~j]) { x~i][seqlj]] 1; ++p[seq[j]]; } delta[i] [seq[j]] = slack [ji; } /* end of m forward-backward knapsacks * bound = objfixed + dual; if (f irsttime) { firsttime = NO; if (bound <= incumbent) return(bound); } opt comp = YES; /*assume optimum completion * for (j1; j<=n; j++) if (jfixed[j]=0) if ((pEji 1) optcomp = NO; break;} if (optcomp){ /* optimum completion, compute the position vector * f or (j=1; j <n; j ++) Pos[j I = j fixed EjJI for (j1; -j<=n; j++) if (jfixed[jl == 0) for (i1l; i<=m; i++) if (x Ei][j]I == 1) { posEj]I = i; break; } checkincumbent(pos); return (bound); } else {/* not optimum completion * 13

/* compute the steepest descent * maximp = -INFINITY; for (j1l; j<=n; j++) if (jfixed[j] == 0) if ((P[j]!= 1) && (j = lastjstar)){ deltal = delta2 = INFINITY; for (i1;- i<=m; i++) if ((p[j]==0) II ((p[j] > 1) && (x[i][jJ == 1)) if (delta[i] [ji < deltal) { delta2 =deltal; deltal = delta[i][jJ; row= else{ if (delta[i] [j] < delta2) delta2= delta[i] [j]; } if ((delta2!= INFINITY) && (delta2 > 0)){ imp = (p[j]==0)? deltal (p[jJ-2) * delta2 + deltal; if (imp > maximp){ maximp =imp; maxstep =delta2; jstar =j; istar = row; steplength = delta2;else if ((imp == maximp) && (delta2 > maxstep)){ maxstep =delta2; j star =;istar = row; steplength = delta2; } /* end of for jloop/ }/* end of else not opt comp * if (maximp == -INFINITY) return(bound); /* zero step length * lastjstar = jstar; if (p EjStar] == 0){ dual -= delta[istar] [jstar]; lambda~jstar] -= steplength; } else{ dual - (p[jstarl-2) * steplength + delta~istar][jistar]l) lambda~jstar] += steplength; bound = objf ixed + dual; if (bound <= incumbent) return(bound); else if (++iter > ((firstdangling == NULL)? ROOTMAMITLIM MAMITLIM)) return (bound) }/* end of while (1) loop * }/* end of MAM */ 14

MAM() calls fbknap() function to solve knapsacks with a post-optimality analysis, i.e., with the calculation of the slack vector below. Input: * nknap = number of items in the knapsack problem, * rhs = RHS capacity of the knapsack, * cknap[j] = objective value of item j, * aknap[j] = knapsack weight of item j, Output * Optimum solution 1 if item j is selected in the optimum xknap[j] 0 otherwise * vknap = optimum objective value, * Slacks ("Aj"s) minimum 6 such that cknap[j]+5 would cause if xknap[j] = 0 xknap[j] = 1 in an alternative optimum slack[j] = minimum 6 such that cknap[j]-b would cause if xknap[j] = 1 xknap[j] = 0 in an alternative optimum See Karabakal, Bean, and Lohmann (1992) for the formulation of forward and backward dynamic programming recursions and for the derivation of the slack vector. void fbknap(nknap, rhs, cknap, aknap, xknap, vknap, slack) int nknap, rhs, aknap[], xknap[], cknap[], *vknap, slack[]; { int f[MAXN+1][MAXB], g[MAXN+] [MAXB]; int soln[MAXN+l] [MAXB], j, beta; /* forward recursions */ for (beta=O; beta < aknap[1]; beta++) { f[1] [beta] = 0; soln[1] [beta] = 0; } for (beta=aknap[1]; beta <= rhs; beta++) if (cknap[1] > 0) { f[1][beta] = cknap[1]; soln[1] [beta] = 1; } else { f[1] [beta] = 0; soln[l] [beta] = 0; } 15

for (j=2; j <= nknap; j++) { for (beta=O; beta < aknap[j]; beta++){ f [j] [beta] = f Ej-1] [beta]; soln[j] [beta] = 0; } for (beta=aknap [j]; beta <= rhs; beta++) if (fU-1] [beta] >= f[j-l1][beta-aknap[j]] + cknap[j]){ f[j] [beta] = f[j-1][beta]; soln[j][beta] = 0; else{ f [j] [beta] = f [j-1] [beta-aknap[j]] + cknap[j]; soln[j] [beta] = 1; } /* determine optimum solution * *vknap = f[nknap][rhs]; beta = rhs; nknap; while (j){ xknap[j] =soln[j] [beta]; if (soln [j]I [beta] ) beta -=aknap Ej]I /* backward recursion * f or (beta=0; beta<=rhs; beta++) g[nknap] [beta] = *vknap; for (j~nknap-1; j >= 1; j —){ for (beta~rhs; beta >= rhs-aknap[j+1]+1; beta —) g~j] [beta] = g[j+1] [beta]; for (beta~rhs-aknap~j+1]; beta >= 0; beta —) if (g~j+1] [beta] < g[j+l][beta+aknap[j+1]] cknap[j+l]) g[j] [beta] g[j+1] [beta]; else g[j] [beta]l g[j+1][be'ta+aknap~j+l]] cknap[j+1]; } /* compute slacks (,A3 quantities) * slack[1] = INFINITY; if (xknap[l] == 0) { for (beta~aknap1l]; beta <= rhs; beta++) if (g[l] [beta] cknap[1] < slack[1]) slack[l] = g[1] [beta] - cknap[1]; } else for (beta=o; beta<=rhs; beta++) if (g [1] [beta] <slack[l] ) slack[El] = E[l] [beta]; 16

for (j=2; j<=nknap; j++){ slack[j] = INFINITY; if (xknap[j] == 0) { for (beta=aknap [j]; beta <= rhs; beta++) if (g [j]I[beta] f f[j - 1][beta-aknap [j] I cknap [j] I( slack [j]) slack[j] = g[j] [beta] - f [j-1] [beta-aknap[j]]- cknap[j]; } else for (beta=0;- beta <= rhs; beta++) if (g [j]I [beta]- f [j- 1] [beta] < slack [j]) slack[j] = g~ji [beta] - f[j-1] [beta]; }/* end of fbknap/ The subgradiento( function implements the subgradient algorithm described in Fisher (1981) for solving Lagrangian duals of integer programming problems. The following step size is chosen for the GAP: steps ize = u (L(A) - incumbent). I- zmx.) int subgradient(jfixed, objfixed, bbar, lambda) int jfixed[], objfixed, bbar[]; float lambda[]; int iter = 0, badcount = 0, pos[MAXN+1], p[MAXN+1], minp[MAXN+1], x[MAXM+1] [MAXN+1], aknap[MAXN+1], xknap[MAXN+1], seq[MAXN+1], i, j, k, rhs, bound; float u, dual, mindual = INFINITY, minlambda[MAXN+1], subgrad[MAXN+1], norm, stepsize, cknap[MAXN+1], vknap; bool optcomp; void knapdpo, lowerboundo, checkincumbento; if (firstdangling = NULL) u = 2.0; else u = 1.0; while (1) { /* initialize dual * for (j=1, dual=0; j<=n; j++) if (jfixed[j]==0) dual += lambda[j]; /* initialize column totals * for (j1; j<=n; j++) p~ji = 0; I 7

/* solve m. knapsacks */ for (=I; i<=m; i++){ k = 0; /* k counts no of elements in knapsack * for (j1; j<=n; j++) if (jfixed[j] == 0) if ( (a [i] [jJ <= bbar [i] ) && ( (c E[i][j]-lambda [j] I> 0. 0)) seq[++k] =; cknap[k]= c[i] [j] - lambda[j]; aknap[k] a= ] j } for (j1; j<=n; j++) x[i] [j] = 0; if (k > 0) { rhs = bbar[i]; knapdp(k, rhs, cknap, aknap, xknap, &vknap); dual += vknap; for (j=I. j<=k; j++) if (xknap[j]) { x[i][Eseq[j]] = 1; ++p~seq[j]]; } }/* end of knapsack solutions * bound = objfixed + (int) dual; optcomp = YES; for (j1l; j<=n; if ((p[j]!= /*assume optimum completion */ j++) if (jfixed[j]0=) 1)) { optcomp = NO; break; } if (optcomp) { /* optimum completion, compute the position vector */ f or (j=1; j <=n; j ++) PosEj I = j fixed [j]I; for (j1; j<=n; j++) if (jfixed[jl == 0) for (i1; i<=m; i++) if (x[i]E[ii == 1) { pos[jl = i; break; } checkincumbent(pos); return (bound); } else {- /* not optimui lowerbound(jfixed, x, bbar, p); if (bound <= incumbent) return (bound); else { /* neither opti rn completion */ imum completion nor inferior */ if (dual < mindual) { /* solution improved */ badcount = 0; mindual = dual; for (j1; j<=n; j++) { minlambda[j] = lambda[j]; minp[j] =pj] } } else { /* solution did not improve */ if (++badcount == SUBPATIENCE) { /* restore last best solution found and continue from there */ /* by halving the step size */ 18

badcount = O; dual = mindual; for (j=1; j<=n; j++) { lambda[j] = minlambda[j]; p[j] = minp[j]; } u = u / 2; } /* find a subgradient */ for (j=1; j<=n; j++) if (jfixed[j] == 0) subgrad[j] = 1 - p[j]; /* compute step size */ norm = 0; for (j=1; j<=n; j++) if (jfixed[j] == 0) norm += subgrad[j] * subgrad[j]; stepsize = u * (objfixed + dual - incumbent) / norm; if (stepsize > ZEROTOL) { /* update multipliers */ for (j=1; j<=n; j++) if (jfixed[j] == O) { lambda[j] -= stepsize * subgrad[j]; } if (++iter < ((firstdangling == NULL)? ROOTSUBITLIM: SUBITLIM)) { continue; /* either zero step size or iteration limit exceeded */ bound = objfixed + (int) mindual; for (j=1; j<=n; j++) lambda[j] = minlambda[j]; return(bound); } /* else neither opt comp nor inferior */ 3 /* else not opt comp */ /* while (1) ) /* end of subgradient */ The knapdp() function, which the subgradient algorithm calls for optimum knapsack solutions, is a linked-list implementation of the "Procedure P2" given by Toth (1980). It is a forward dynamic programming algorithm and uses the "further but cheaper"-type pruning technique to improve performance. Note that, under this pruning, the slack vector cannot be calculated and thus, this procedure is not suitable for solving the knapsack subproblems of the MAM(). Incidentally, the "Procedure P2" of Toth (1980) has a bug in it; it is corrected in the implementation here. The input-output parameters defined for function fbknap() are also valid for knapdp(). 19

void knapdp(nknap, rhs, cknap, aknap, xknap, vknap) int nknap, rhs, aknap[], xknap[]; float cknap[], *vknap; { int L[MAXN+1][MAXB], val[MAXN+l] [MAXB], prev[MAXN+l][MAXB], s[MAXN+l], h, k, j, y, xtemp, ptemp, mm; float F[MAXN+1] [MAXB], ftemp; /* initialize for mm = 1 */ L[l][O] = F[l][0] = val[l][0] = prev[i][0] = s[ll = 0; if ( (aknap[l] <= rhs) && (cknap[l] > 0).) { L[l][1] = aknap[1]; F[] [1] = cknap[l]; val[l][1] = 1; prev[l][1] = O; s[ll = 1; } for (mm=2; mm <= nknap; mm++) { L[mm][O] = F[mm][O] = val[mm ][0] = prev[mm][0] = 0; h = ( (s[mm-] == 0)? 0: 1); k = j = 0; y = aknap[mm]; while (1) { if ( L[mm-l][h] < y ) { if ( F[mm-l][h] > F[mm] [k ) { k++; L[mm] [k] = L[mm-l][ h]; F[mm] [k] = F[mm-l] [h; val[mm][k] = 0; prev[mm][k] = h; } if ( h == s[mm-1] ) break; else h++; } else { if ( L[mm-l][h] > y ) { ftemp = F[mm-l] [j + cknap[mm]; if ( ftemp > F[mm][k] ) { k++; L[mm][k] = y; F[mm][k] = ftemp; val[mm][k] = 1; prev[mm][k] = j; } y = L[mm-1][++j] + aknap[mm]; } else { /* L[mm-l][h] = y*/ ftemp = F[mm-l][j] + cknap[mm]; xtemp = 1; ptemp = j; if (F[mm-l][h] > ftemp) { ftemp = F[mm-l][h]; xtemp = 0; ptemp = h; } if (ftemp > F[mm][k]) { k++; L[mm][k] = y; F[mm] [k = ftemp; val[mm][k] = xtemp; prev[mm][k] = ptemp; } y = L[mm-l][++j] + aknap[mm]; if (h == s[mm-1]) break; else h++; } } } /* while (1) */ 20

/* step 3 */ while (1) { if Cy > rhs) break; else { f temp = Fllimm- 1][j]I + cknapli[mm]; if (ftemp > F[nm]l[k]){ k++; L[mmllk] = y; F[mm][k] = ftemp; val[rnm][k] = 1; prevlmm][k] = j if (j == S[mm-1]) break; else y = L[mm-l][++j] + aknap[rnm]; } *while (1) * /* step 4 */ s[mm] = k; } /* for (mm=2... * 1* trace links back to recover solution * *vknap = F [nknap] [s[nknap]]; k = s [nknap]; for (mm=nknap; mm >= 1; mm —){ xknap [mm] = val [mm] [k]; k = preyv[mm] [k]; } }/* end of knapdp/ 3.4 Lower Bounding Function (Heuristic) Given a Lagrangian solution, x, that satisfies knapsack and integrality constraints, but violates multiple-choice constraints, the following heuristic attempts to obtain overall feasibility- by modifying current assignments. void lowerbound(jfixed, x, bbar, p) mnt j fixed [], x [][MAXN+1], bb ar[]l, p ]; mnt pcopy[MAXN+1], xcopyEMAXM+1] [MAXN+1], capacity[MAXM+1], pos [MAXN+1], aknap[MAXN+1], xknap[MAXN+1], seq[MAXN+1], i, j, k, rhs, istar, feasible; float minratio, temp. cknap[MAXN+1], vknap; void knapdpQ), checkincumbent(); /* copy p to pcopy, x to xcopy * for (j1; j<=n; j++){ pcopylj] = pEji; for (i=1; i<=m; i++) xcopy [ili [ = xli]E[ii } 21

/* compute available knapsack capacities */ for (i1l; i<=m; i++) { capacity Ei] = bbarEi]; for (j1; j<=n; j++) if (jfixed[j] == 0) if (xcopy[i] [ji == 1) capacity[i] -= a[i] [j]; /* Order the variables x, such that xi I and p3- > 1 in the current solution * in ascending ci -/ai -. Following this order, set xij = 0 and increase the'slack space of knapsack i by a l,. Repeat until p, ~ 1 for all j. for (j1; j<=n; j++) if (jfixed[j] == 0) while (pcopy[j] > 1) { minratio = INFINITY; for (i1; i<=m; i++) if (xcopyllij[j] == 1){ temp = (float) c~il [ji / (float) a[i] [j]; if (temp < minratio){ minratio = temp; istar =; xcopy[istarjlj] = 0; — pcopy[jl; capacity[istar] += a[istar] [jJ; /* Using only variables XI~ with p3 0. solve all knapsack problems optimally * for (i1; i<=m; i++){ k = 0; /* k counts no of elements in knapsack * for (j1; j<=n; j++) if (jfixed[j] == 0) if ( (pcopy[j] = 0) && (a[i] [j] <= capacity[i])){ seq[++k] = cknap [kl = c Ei][l [Iaknap [k], = a [i][Ej]I if (k > 0){ rhs = capacity[i]; knapdp(k, rhs, cknap, aknap, xknap, &vknap); for (j1; j<=k; j++) if (xknap Ej] == 1) { xcopy[i] [seqEji] = 1; ++pcopy~seq~j]]; } /* end of knapsack solutions * 22

feasible = 1; for (j1; j<=n; j++) if (jfixed[jJ == 0) if (pcopy~jj!= 1) { feasible = 0; break; } if (feasible) { for (j=1; j<=n; j++) pos ji = jfixed[j]; for (j1; j<=n; j++) if (jfixed[j] == 0) for (i1; i<=m; i++) if (xcopyli][j] == 1) { pos[j] = i; break; } checkincumbent(pos); } /* end of lowerbound */ 3.5 Fathoming Functions Fathoming functions are designed to initialize parameters and organize calling sequences. They return one if the node is fathomed, zero if the node is dangling. Fathomed nodes are cleared from the memory immediately, whereas dangling nodes are inserted into a priority queue ordered by upper bounds. bool fathomroot() NODEPTR root; int jfixed[MAXN+1], i, j, maxi, max2, omin, omax, bound; bool fathom, fathomrodeo; void putdanglingo; /* determine absolute lower and upper bounds */ absib = bound = 0; for (j1; j<=n; j++) { cmin = INFINITY; cmax = -INFINITY; for (i1; i<=m; i++) { if (c~i][j] < cmin) cmin = c~i][i]; if (c[i] [ji > cmax) cmax = c~il [ii; absib += cmin; bound += cmax; incumbent = absib; /* create the root node */ if ((root = (NODEPTR) malloc(sizeof(NODE))) == NULL) { printf( "Cannot create root node."); exit(1); } root->nodeno = 0; root->parent = NULL; 23

/* initialize multiplers by selecting the second maximum objective value * /* of the variables from each multiple-choice set * for (j1; j<=n; j++){ maxl =max2 = -INFINITY; for (i1; i<=m; i++) if (a[i][jJ <= b[i]){ if (c[i][j] > mail) { max2 = maxi; maxi c[iJ[j]; } else if (c~i][j] > max2) max2 = c[iJ[j]; root->lambda[j] = (max2 == -INFINITY)? 0 max2; if (PRINTSTEPS) printf ("Node 0 D'.d, %.dJ => ", incumbent, bound); for (j=l; j<=n; j++) jfixed[j] = 0; /* no variables are-fixed at the root/ fathom = fathomnode(root, jfixed); if (!fathom) putdangling(root); return (f athom); }/* end of fathomroot * bool fathomnode(bbnode, jfixed) NODEPTR bbnode; mnt jfixed[]; int objfixed, bbar[MAXM+l], bound, i, jmamiter, subiter, ilambda[MAXN+l], 1UB; float dlambda [MAXN+l];.bool fathom; NODEPTR hold; mnt MAMO, subgradientO); void printsolO); if (bbnode->nodeno){ ++nbranch; mamiter = MAMITLIM; subiter =SUBITLIM; } else{ mamiter = BROOTMAMITLIM; subiter = ROOTSUBITLIM; } /* add node's fixed var to jfixed * if (bbnode->nodeno > 0) jfixed~bbnode->jfix] = bbnode->if ix; 24

for (P1; i<=m; i++) bbar[i] = b~il; objfixed = 0; for (j=1; j<=n; j++) if (jfixedlljj){ bbar[ jfixed~j] J =a[ jfixed~j] I[j]; objfixed += cE jfixedlij] ][jI; } for (j1; j<=n; j++) ilambda[j] = bbnode->l~ambda[j]; if (mamiter > 0) bound = MAM(jfixed, objfixed, bbar, ilambda); if (bound > incumbent) { for (j1; j<=n; j++) dlambda[j] = (float) ilambda[j]; if (subiter > 0) bound = subgradient(jfixed, objfixed, bbar, dlambda); if (bound > incumbent) { for (j1; j<=n; j++) bbnode->lambda[j] = (int) dlambda[j]; bbnode->bound = bound; fathom = (bound <= incumbent)? YES:NO; if (PRINTSTEPS) { printf("[%.d,%d] " incumbent, bound); if (fathom) printf("1FATHOM\n"); else printf("'DANGLING\n"); } if(nbrarich >= MAXBRANCH){ printf ("Branch limit exceeded.\n"); if (incumbent == abslb) { printf("No feasible solution found."); } else{ printf("\nBest-solution found = 7.d, ",incumbent); if (bbnode->nodeno == 0) UB = bound; else{ UB= (firstdangling == NULL)? -INFINITY firstdangling->bound; if (bbnode->ifix!= m){ hold = bbnode->parent; if (hold->bound > UB) UB = hold->bound; printf("Upper bound = %d", UB); printf("l (Error bound =42M n ( 100 * (float) (UB- incumbent)/ (float) incumbent)); printf("No of branches =%d\n", nbranch); if (PRINTOPTIM) printsol(xinc); exit (0) 25

/* correct jfixed */ if (bbnode->nodeno > 0) jfixed[bbnode->jfix] = 0; return(fathom); } /* end of fathomnode */ 3.6 Branch-and-Bound Function A node selected from the list of dangling nodes must have a solution violating one or more multiple-choice sets. At the beginning, the only dangling node is the root. Select a violated multiple-choice set with largest multiplier. Since any feasible solution includes exactly one out of m variables from this constraint, create m branches (Bean 1984). At each branch, one variable is fixed at one, and the rest are automatically set to zero. void BAB() { NODEPTR cnode, kid, search, getdangling(); int jfixed[MAXN+1], used[MAXM+1], i, j, jstar, maxlambda, serialno = 0; bool feasible, fathomnode(); void putdangling(, checkincumbent(); while (firstdangling!= NULL) { cnode = getdangling(); /* determine fixed and free vars */ for (j=l; j<=n; j++) jfixed[j] = 0; search = cnode; while (search->parent!= NULL) { jfixed[ search->jfix ] = search->ifix; search = search->parent; } /* select jstar = multiple choice set to branch (one having biggest multiplier) */ maxlambda = -INFINITY; jstar = 0; for (j=l; j<=n; j++) if (jfixed[j] == 0) if (cnode->lambda[j] >= maxlambda) { maxlambda = cnode->lambda[j]; jstar = j; } /* determine used knapsack capacities by fixed variables */ for (i=l; i<=m; i++) used[i] = 0; for (j=1; j<=n; j++) if (jfixed[j] > 0) used[jfixed[j]] += a[jfixed[j]][j]; 26

if (jstar == 0) { /* hit the bottom, jfixed is nonzero feasible = YES; for (i1; i<=m; i++) if (used~il > b[iJ) { feasible = NO; break; } if (feasible) checkincumbent~jfixed); else /* create m branches * for (i1; i<=m; i++) if ((used[i]+a[i][jstar] <= b[i]) && cnode->bound > incumbent){ if ((kid = (NODEPTR) malloc(sizeof(NODE))) == NULL){ printf( "Cannot create new branches."); exit (1) kid->nodeno = ++serialno; kid->Parent = cnode; kid->ifix = i; kid->jfix = jstar; for (j=1; j<=n; j++) kid->1ambda[j] = cnode->lambda[j]; if (PRINTSTEPS) printf ("Node Yd From %d: x(%d,%d) [%d,%d] ==> " kid->nodeno, cnode->nodeno, i, jstar, incumbent, cnode->bound); if (fathomriode(kid,jfixed)) f ree ((char *) kid); else /* node dangling * putdangling(kid); 1* insert into queue * }/* end of while * }/* end of BAB * 2'T

3.7 Priority Queue Handling Functions Three functions are needed for handling the priority queue of dangling nodes. 1. The function putdangling() puts a new node into the correct place so that the order of decreasing upper bounds is maintained. void putdangling(newnode) NODEPTR newnode; { NODEPTR search, last; newnode->next = NULL; search = firstdangling; if (search == NULL) firstdangling = newnode; else { while (search!= NULL) if (search->bound > newnode->bound) { last = search; search = search->next; } else break; newnode->next = search; if (search->nodeno == firstdangling->nodeno) firstdangling = newnode; else last->next = newnode; } } /* end of putdangling */ 2. The function getdangling() gets a dangling node with highest upper bound. NODEPTR getdangling() { NODEPTR hold = firstdangling; firstdangling = firstdangling->next; return(hold); } /* end of getdangling */ 3. The function removeinferiors() checks the list of dangling nodes for possible fathomes whenever a better incumbent is found. 28

void removeinferiorso0 NODEPTR last, prey; while (1) if (firstdangling == NULL) return; else { last = f irstdangling; preyv NULL; while (last->next!= NULL){ preyv last; last = last->next; if (last->bound <= incumbent) { /* delete last * free((char *)last); if (prey == NULL) {/* last was the only one firstdangling =NULL; return; } else Prev->next = NULL; else return; } /* end of removeinferiors * 3.8 Main ]Function and Utilities The function printsol() prints an input position vector (feasible solution,) in an easy-tovisualize format. void printsol(pos) mnt posE]; mnt i, j, used[MAXM+1], cx=0; for (i1l; i<=m; i++) usedEil = 0; for (j1; j<=n; j++){ i = pos~ji; used~il += a~i]Ej]; cx +=ciE; printf("\n cx = %d\n ",'cx); for (j1; j<=n; iji-i) printf("%ld", (j % 10)); Printf(" USED RHS\n"); for (i1; i<=m; i++){ Printf ("7.2d)", i); for (j1; j<=n; j++) if (pos[jJ == i) printf("1"); else printf("l."); printf("I %4d 7.3d\n", used~il, b~i]); } } /* end of printsol * 29

The function checkincumbent () first computes the objective value of the input position vector (feasible solution). If it is better than incumbent, updates incumbent and fathoms inferior dangling nodes. void checkincumbent(pos) int posE]; int j, vpos = 0; void removeinferiorso; for (j1; j<=n; j++) vpos += c~pos~ji] jl; if (vpos > incumbent) { incumbent = vpos; for (j1l; j<=n; j++) xinc[j] = pos[j]; removeinferiorso; } } /* end of checkincumbent */ void main(argc, argv) int argc; char *argv[j; char inpname[50]; void readino, BABO, setparamso, printsol(); bool fathomrootO; setparams(); if (argc == 2) strcpy(inpname, *++argv); else inpname[0] = readin(inpname); /* Solve the root problem; if duality gap, continue with branch-and bound */ if (!fathomrooto) BAB(); printf("\nOptimum = %d\nNo of branches = %d\n", incumbent, nbranch); if (PRINTOPTIM) printsol(xinc); } /* end of main */ 30

References [1] Bean, J. C., "A Lagrangian Algorithm for the Multiple Choice Integer Program," Operations Research, 32 (1984), 1185-1193. [2] Fisher M. L., "The Lagrangian Relaxation Method for Solving Integer Programming Problems," Management Science, 27 (1981), 1-18. [3] Fisher M. L., R. Jaikumar, and L. N. Van Wassenhove, "A Multiplier Adjustment Method for the Generalized Assignment Problem," Management Science, 32 (1986), 1095-1103. [4] Guignard, M. and M. B. Rosenwein, "An Improved Dual Based Algorithm for the Generalized Assignment Problem," Operations Research, 37 (1989), 658-663. [5] Karabakal, N., J. C. Bean, and J. R. Lohmann, "A Steepest Descent Multiplier Adjustment Method for the Generalized Assignment Problem," Technical Report 92-11, Department of Industrial and Operations Engineering, The University of Michigan, Ann Arbor, 1992. [6] Martello, S. and P. Toth, Knapsack Problems: Algorithms and Computer Implementations, John Wiley and Sons, 1990. [7] Toth, P., "Dynamic Programming Algorithms for the Zero-One Knapsack Problem," Computing, 25 (1980), 29-45. 31