THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING APPLICATION OF MATHEMATICAL PROGRAMMING TO THE SELECTION OF EXPERIMENTAL FACTOR ARRANGEMENTS WITH RESOURCE CONSTRAINTS John B. Neuhardt A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Industrial Engineering 1967 July, 1967 IP-785

PREFACE I wish to acknowledge the exceptional guidance provided by my chairman, Professor Hugh Bradley, and the continuing efforts of my committee members to provide help whenever it was asked. I further wish to thank Professor Paul Dwyer for initial reading of these ideas, and to acknowledge the helpful suggestions by Professor E. Peterson and Professor Bernard Galler. I would like to express gratitude to Professor Richard Wilson, Professor Dean Wilson, and Dr. Wilfred Bychinsky for direction and influence in my student graduate career. I also wish to thank those responsible for my having been chosen recipient of a Rackham Pre-doctoral Fellowship, and a General Motors Student Fellowship. Finally, I thank my family for confidence, understanding, and infinite patience. ii

TABLE OF CONTENTS Page PREFACE................................... ii LIST OF TABLES............................. v LIST OF ILLUSTRATIONS..................... vi CHAPTER I. INTRODUCTION........................ 1.1 The general problem.............. 1.2 A simple formulation............ 2 1.3 Several definitions of the information of a design............... 3 1.4 A typical problem................ II. A MINIMUM EIGENVALUE CRITERION........ 10 2.1 Detailed problem statement.......... 10 2.2 The constraint region and properties of the objective function.......... 12 2.3 The one-factor problem.......... 14 2.4 The gradient search procedure........ 16 2.5 Obtaining an initial feasible solution... 18 2.6 A gradient approximation............ 20 2.7 Three search illustrations......... 21 mI. A SUM OF SQUARES CRITERION.... 31 3.1 Sum of squares criterion defined; rationale............ 31 3.2 A property of the objective function... 33 3.3 The noninteger quadratic program...... 35 iii

CHAPTER Page 3.4 An algorithm of Kunzi and Oettli....... 38 3.5 Witzgall's parabolic constraint algorithm................. 46 3.6 A partial enumeration method of Lawler and Bell; an example....... 49 3.7 Use as a starting point for the minimum eigenvalue search problem..................... 53 WI. FURTHER EXAMPLES OF THE SUM OF SQUARES CRITERION......... 54 4.1 The equal-cost problem............ 54 4.2 Three preliminary observations........ 54 4.3 The one-factor problem............ 60 4.4 The latin square............... 62 4.5 The balanced incomplete block........ 63 4.6 Webb's problem................... 64 V. SUMMARY AND AREAS OF FUTURE INVESTIGATION...................... 67 5.1 Summary........................ 67 5.2 Areas of future investigation....... 69 APPENDIXES................................. 71 A. Eigenvalue Properties..................... 72 B. Obtaining a Feasible Solution.............. 75 C. Computational Considerations............... 94 REFERENCES............................... 102 iv

LIST OF TABLES TABLE Page 1. Integer Three-Factor, Two-Level Search Problem............................. 23 2. Noninteger Three-Factor, Two-Level Search Problem........................... 24 3. Integer Search, Minimizing Maximum Eigenvalue.......................... 27 4. Integer Search, Total Number of Observations Large........................ 29 5. Comparison of Eigenvalue Calculations for Various Matrix Dimensions............... 96 6. Experience in Linear Program Calculation Times.............................. 97 7. Optimistic Computation Time per One Step in Gradient Search..................... 97 V

LIST OF ILLUSTRATIONS FIGURE Page 1. Webb's solution for a 4x4x3 problem, with twelve observations.................... 26 2. The quadratic program surface............... 41 3. Illustration of the enumeration technique....... 52 4. Linear model........................... 77 5. Additive model.......................... 85 vi

CHAPTER I INTRODUCTION 1.1 The general problem Stemming from the early writings of R. A. Fisher, developments in the field of experimental design have advanced enormously. Because many authors have shown existence properties for various criteria of optimality, the scientist or engineer can now choose from an abundant storehouse of designs, many of which possess certain optimum properties. However, when the environmental constraints do not readily admit to the selection of any particular known "good" design, the next step is not always apparent. This may be particularly true if the experiment is constrained by cost considerations. Very often in industrial research, an experimental program must be planned for which... existing designs are inadequate. In such cases the consulting statistician may have a tendency to try and alter the thinking of the experimenter so that one of the standard published designs can be used. (51) Thus, as an expedient measure, the environment may be forced into a theoretical mold, a mold which may fit quite poorly. Even if the environmental constraints are consistent with some known "optimal" 1

2 design, the experimenter may not have access to a description of this design that meets his requirements. One alternative to this problem is to enumerate all possible designs feasible within the constraints. However, such a list could become extremely large, even in the smaller experiments. Another approach would be to develop an iterative procedure that ended in a reasonably good design within the resource constraints. Other than an integer programming formulation for the problem of finding orthogonal latin squares (16, 3), the only work known to this author involving the mathematical programming approach is that of Webb (52), who describes a linear programming formulation for the construction of orthogonal designs. Here, the objective was to minimize the total number of experimental observations. 1.2 A simple formulation Suppose a factorial design problem has potentially p unique factor combinations that could be observed. Denoting Euclidean P-space as E(P) all vectors X6E(P) which have integer-valued nonnegative components (each component corresponding to a particular A simple 23 factorial design with eight factor combinations that has a constraint of a maximum of five observations per combination could conceivably have 58 different possibilities.

3 factor combination) are possible solutions to the design problem. Defining I(X) to be the information associated with the vector X, and U the subspace of E( representing the resource constraints, the problem is then stated: Maximize I(X) Subject to X U XC > ~0 integer (;,,..., p) Here, Xr is the it component of X. In all past work in the literature,I(X) was a nondecreasing function of X. The objective functions considered in Chapters II and III also satisfy that condition. 1.3 Several definitions of the information of a design The information of a design has been defined in many ways. To discuss some of these briefly, we consider the factorial model for the expected value of the observation vector Y e E( E(Y) = ZQc with oc the vector of unknown constants and Z a matrix of known constants related to X over which the experimenter has control. The exact relation of Z and X is shown in Chapter II.

4 ~tf\ t Least squares estimates,LoC of the linear combinations L;oc are usually of interest, and most optimality criteria are directed toward the variances of these estimates. It is also assumed that the covariance matrix of Y is I, where I is te identity matrix of suitable dimension. For the moment it will be assumed that S=Z Z is nonsingular. In general, however, the matrix S will be singular in this dissertation. When S is nonsingular, there exists a set of LoC< which have variances V; =- S;^/det(S)) L- Is..., p, s;:is the cofactor of S; in S-(sij), and "'et" denotes the determinant operator. Under certain restrictions on Z, Plackett and Burnam (42) develop designs which simultaneously minimize the VY. (Of course, simultaneous minimization of the V; is generally not possible, but the authors do establish this property for a class of additive designs for which the number of factor-levels are all equal.) If the generalized variance is defined as the product of the V;, this quantity (except for the multiplicative constant T ) is equal to the product of the eigenvalues of the inverse Later, the one-factor case will be an exception, in that the nonsingular matrix will greatly facilitate discussion. Although the oC and Z may be transformed in any model to achieve a nonsingular S, interpretations become more difficult and hence this transformation will be avoided.

5 of S, or S, which in turn equals the determinant of S. This latter quantity divided by the maximum obtainable was defined by Wald as the design efficiency; latin square arrangements are shown to have efficiency equal to one (50). Defining average variance to be the sum of the vr, Kempthorne (35) shows the relationship of this quantity to the harmonic mean of the eigenvalues of S. Kshirsagar (36) extends these results, showing the relationship of the generalized and average variances, and discusses optimal properties of the balanced incomplete block and Youden Square designs. Turning to the case of a possibly singular S, Ehrenfeld (19) defines the efficiency of a design to be the ratio of the minimum nonzero eigenvalue of S to the maximum such value obtainable over all feasible arrangements. It is this criterion that will be used in Chapter II. Webb (51) defines an "estimation index" proportional to the generalized variance and a "fitting index" inversely proportional to the average variance. He calculates these indexes for various designs.

6 1.4 A typical problem For illustrative purposes, a typical problem will be posed and its characteristics determined throughout the general problem development. Suppose a manufacturer of ceramic parts is interested in maintaining the average level of an observable characteristic y, a quality measure of the strength of the manufactured material. Suppose that three production variables (controllable) are known to affect the strength, and these are (a) the pulverizing process of the raw materials (process 1 or 2), (b) the pressure (PSI) under which the parts are formed, and (c) the oven temperature at which the parts are "baked." Assuming that only two values (or levels) of each of the three independent variables are of interest, an additive linear model is assumed: E(yj) M u + vj wK (ti,,kv 1IZ) where E is the expectation operator, i is a constant strength level uT is the differential effect on strength of process one or two ( = 1, 2) V, is the differential effect on strength of pressure level low or high (j = 1, 2) Ow is the differential effect on strength of oven temperature low or high (k= 1, 2)

7 All the independent differential effects are assumed unknown parameters. It is then desired to estimate linear relations in the parameters of the type (U1-, V2-v) wofw), here the differential effects of the independent variables on strength. If several observations of yi;, are to be taken, we label such a experiment in the following tabular manner: VI V2 vI V2. W wq WI W2 In this table, X.(i = 1,2,...8) represents the number of observations that are to be taken with the corresponding settings of the independent variables A, v, w, as noted. For instance, x5 represents the number of observations with the (u,)V, ) 2) factor combination, or variables u, v,w, at levels 1, 1, 2, respectively. Defining the vector X: ( xX - )..x) E(, for any information function 1(X), nondecreasing in X, we would like to increase X without limit to accomplish a maximum in information, I(X). However, restrictions of some form always exist on X. Classical techniques treat this problem (of the general class called 2 factorials) with the restriction: ~x; 8S for =' 2', (here, rnO, I). For

8 rn I, these are fractional factorials in common use (34, 15, 1, 14, 11). In the example, the classical design with known optimal properties is r 1/2, an orthogonal one-half replicate. There are two such arrangements, XO= (1,0,O,1,0,I,1,o) and E - 0, where E5 is a vector with all elements unity. Suppose, however, that different types of constraints were present in the form (x + x2 x; + _ 3), and (t'2 x5 ) a < _2 These types of restrictions will be called resource constraints, the two listed representing an upper limit on the total observations available for variable w at level 1, and variable A at level 1, respectively. Also, if we attach a positive number c; to the ith factor combination, representing the cost or time to take an observation (increment X; by 1), the cumulative experimental cost or time may be constrained by the value Co. It is these types of problems that will be of interest during the next chapters. It is noted that the constants (uA; vj, w) imply a "fixedeffect" model, assumed throughout. Unless specified, it will generally not be assumed that the model is additive (separable in the i, j, ). Also, the central problem is that of selecting a layout, combination, or arrangement of factor levels which will ultimately be observed in a completely randomized design, hence the errors will always be assumed uncorrelated. For example, we consider

9 properties of latin square "arrangements" rather than latin square "designs."

CHAPTER I A MINIMUM EIGENVALUE CRITERION 2.1 Detailed problem statement The linear estimation problem is examined in more detail in this section. The linear model E(Y):Zo is postulated with Y C(. a random observation vector with covariance matrix T C I, oc (mx a vector of unknown constants, and Z(Nxm) the design matrix over which the experimenter has control. Given Z, an (rmxm) information matrix S:Z Z is defined, and relative to the set of linearly estimable functions V=lLi; L; \\}, Ehrenfeld (19) shows that the t A variance of least squares estimators of elements of V, L.oC, attains a maximum value T/A, with A the minimum nonzero eigenvalue of S (provided the rank of S is a maximum r <_ r, taken over X;, for a given row space of 2). The eigenvalues of S are taken to be A\ >^2>-.., _>A. When r=v, S is nonsingular, again not generally assumed here. i t: t L Oc is linearly estimable if there exists an estimate L;c, linear in the observations, such that E~L.T^~o- Lifo 10

11 The nature of the factor arrangement problem is now considered. The linear model includes a definition of several factors, each of which can be controlled at one of several "levels." A particular factor-level combination can be thought of as corresponding to an element of the index set = 1) 2,...,p. A set of column vectors {i; IEE(m) iE. is defined whose elements are 0 or 1. Each; becomes a row of Z, specifying which elements of o are included in an observation of Y. Further, the selection by the experimenter of the integers x;i E I} represents the number of observations that will be taken at each factor combination. Thus, the factor arrangement to be observed is the specification of the vector {x; = X 6 (P) The matrix Z may be structured as follows: z-= (1, 1, ),..S 0,2, *'",,-"'"'~ fp.,p)t with each;i repeated x; times so that It is seen that the minimum nonzero eigenvalue, of is then a It is seen that the minimum nonzero eigenvalue, Air of Sis then a function of X and the information of such a design is then defined as A (X). The mathematical problem to be solved is stated:

12 Maximize /r(X) Subject to XEU Xi > O0 integer (i=l,...,p) 2.2 The constraint region and properties of the objective function Several pertinent properties of the function /2(X) are now listed. The eigenvalues of S are defined to be Al > A?> -. r, >..?Ap_ with A,> if S is of maximal rank r, and A=O(=r+l,..,]r in any case. (Note that S is a positive, real symmetric matrix). As is well known, Ar(X)is a continuous function of X (54). Other properties include: (a) Ar(X) is nondecreasing in elements of X; i.e., if Xi is element by element greater than X2, then Ar ( XI). (X,) (b) If sj is the (i j) element of S, then <(X) < Minimum {si,; (c) (X)=< X, for scalar k< (d) A; = Trace(S) v x;, where V is the (constant) number of l's in any (e) A = Trace( Se ) 2 sIj (f) If the linear additive model has g factors, each factor having nr' levels ( i X,,...,g) and if h is the maximum

13 of the r;, then /r(X)_ N//, where N = Z % (g) /Ar(X _E([ x;Y)...2A( 0] \ 2 (Reference 12) (Properties b, d, and f are discussed in Appendix A.) The constraint region U will be assumed to be defined by constraints linear in X. Two types of constraints are considered: resource limitations and an upper bound on total experimental cost or time. The resource limitations are of the form AX b, with elements of A zero or one and b integer-valued. To each factor combination is attached a cost or time per observation, C;, and total experimental cost or time expenditure is assumed to be t7xc;, which has an upper bound constraint Co. Consider the example of section 1.4. Suppose that the total observations of factors represented by L, V, and W, at their "low" and "high" levels is bounded by (b, b),(b, b.), and (b b6, respectively. Given the {c;( and the constant Co, the constraint region U would be defined by the inequalities: Xi8+ X27 + 8 ba (3 + 5x + x7+8 Xt + 3 tx5 +X7 3 ( X4 ++x6 N +8 X <w, +'<2 +3 6+ 7 bW

14 Each resource constraint represents the total number of observations allowable for each factor level. 2.3 The one-factor problem This case may be solved by inspection. If the original model E-(y)-q has a singular S-Z:Z, this corresponds to a parameter vector oC= (/fi>,B,) ^ 0 ~**, **p ). Defining: A'+ Ja 2j;= 2 ) *3 ~.'- ar... * ) p and Z to be the matrix Z with the first column removed, we have: (y)Z Zoc =Z O Further, S'=(Z )t(Z, ) will be a diagonal, nonsingular matrix since now each row of Z' has one unit entry with the remaining components zero. The element s$; corresponds to the number of observations of the ith level (or the number of measurements involving +-A^? As mentioned in Chapter I, only in the one-factor problem will the concept of this transformation be utilized. Although textbooks treat such a case for an arbitrary number of factors, the general case of a singular S will be treated for two or more factors.

15 and is exactly the eigenvalue Ai of S (ir,,.. p). As mentioned above, the resource constraints are assumed to be limitations on the total observations for each parameter or factor level, so that this problem may be stated: Maximize W Subject to W < X; (- )2).., p) x >O0nteer(i,2,..,p) x<62; (Al,2,..' adp; Integer) With d the minimum of the J,}, k the largest integer for which ZC {< Co/K,2 the minimum of (dk) and Xothe vector with all components unity, a solution is n ~. It is noted that this solution is not necessarily unique since, for a solution {X?3 XO that has LC( X~ strictly less than C,, some of the {X} could possibly be increased within the constraints even though the objective function, W, remained unchanged. Essentially, all the resources might not be "used up." This shortcoming could be remedied somewhat by considering the problem in stages, the nr stage being the above problem with the additional constraint, Z Xi- n. By solving the stage with resultant A\r, n could be increased until the feasible set was empty, selecting the solution of the nm stage as optimal, where

16 m is the largest integer such that At > A4, for all stages q in the sequence. 2.4 The gradient search procedure The problem of section 2.1 is not a typical mathematical program. The objective function is not known to be convex, and the feasibility region consists only of vectors with integer coordinates. Standard mathematical programming techniques would appear to fall short in solving this type of problem, with the possible exception of a search technique or enumeration methods of integer programming (38, 7). The latter will be explored in a later chapter. For the present, the gradient search is discussed. Since the feasibility region is integer-valued, standard approaches are not available. Only in the one-dimensional case are such integer or "lattice" search problems considered (53). However, the gradient technique might be applicable if the objective function were relatively "linear" for integer steps. Presuming that this might be the case, the search procedure is reviewed. The well-known ingredients of a gradient search procedure are the following (see, for instance, ref. 23): (a) An initial feasible solution (b) For every X, an approximation of the gradient O7 t(X)

17 (c) Calculation of a feasible direction vector dX, where V7r' cX>O and X^)SdXA U, for some scalar S. (d) Stopping rule(s). The general premise in the gradient approach is that the step size S, when confined to values for which cdX is integer, is small enough so that ^(X)+^V,54Xis approximately equal to )r(X+S4), As will be seen, this is generally not true and the gradient search procedure confined to integer X has limited value. Relative to the stopping rules, since the objective function is not known to be convex, an Xo which is a good local optimum need not necessarily be a solution to the global problem. However, if A,(X) attains a known upper bound, obviously the corresponding X is optimal. For instance, the upper bound of property (b) in section 2.2 might be attained. Also, presuming that we are at a feasible point X, that is where AX b, CX co, and given the gradient u=VAX)r(X)=D\r (X),..,^D((r()Hadley (32) suggests solving the following nonlinear program for a "step" y given some k: Maximize t y Subject to A(X+Y) b ctx*+ ) ( co

18 If the solution A 0o is zero, search could be halted with this local maximum. 2.5 Obtaining an initial feasible solution For convenience in presentation in the remainder of this chapter, it will be assumed that the additional constraint, /r(X) >0, is placed on the solution X. Because it is desired to keep the time-consuming eigenvalue calculations (Appendix C) to a minimum, methods for obtaining a "good" initial solution are now discussed. The problem of obtaining a feasible solution is to find an X so that the four constraints are satisfied: (a) AX0Ocb oj=oJI (; i-);jl,..,p), X (b) CtX~< C CJ > O (j-., 1P), X~ O (c) X e r (i-,, ) (d) t (X ) > 0 Should either constraint (a) or (d) be eliminated, advantage can be taken of the resulting problem structure. For instance, consider This constraint is superfluous in the maximization problem, but is necessary in the gradient search method since the gradient approximation that will be used is valid only when Ar >0

19 the problem of obtaining an X satisfying constraints b, c, and d. Note that requirement d is satisfied only if the space generated by the set {; X, ois of dimension r, i.e., there must be r linearly independent vectors corresponding to the nonzero Xi. Therefore, if we can choose a minimum cost collection,', of r vectors from the; and set X = if;67, with Xj'O otherwise, the resultant X> (X~.. X) is a feasible solution (if C XO Co) or no feasible solution exists (C X > o) since X is the minimal cost vector subject to constraints b, c, and d. A method for finding the set r may be found in Appendix B. Suppose now that constraint (d) is ignored and again we seek a minimum cost solution and formulate the problem (inserting equalities in the resource constraint for simplicity): (2.5.1) Minimize CtX Subject to AX =b X;_O,integer (i-,)2,..)p) If, for the solution X to this integer program, CX > Co, obviously no feasible solution exists. On the other hand, if C X~Co and )((X~)>Oone has obtained a feasible solution. One final observation is made relative to this last approach. If the model E(&O):Zo is assumed additive, each vector of the set

20 {;} (potentially a row of Z), corresponds precisely to a column of the matrix A above. This is due to the way the resource constraints were defined. As a result, the constraint Ar(X)>0can be adjoined in (2.5.1) by requiring instead that the solution be nondegenerate. It is possible that an algorithm could be devised to enforce such nondegeneracy, although no such treatment is known to this author. Algorithms for the other extreme, enforcing a high degree of degeneracy, have been devised for the fixed-charge linear programming transportation problem. 2.6 A gradient approximation For a gradient calculation, Goertzel and Tralli (26) consider the matrix S with an eigenvalue )>0 and corresponding eigenvector P of unit length, perturbed by bS, resulting in a corresponding eigenvalue A= A A and eigenvector P +hP, so that: (2.6.1) S P A P (S + S) (P+ aP)= (A +A)0 (P+P) Expanding the latter expression, simplifying due to equation (2.6.1), and multiplying by P results in: Pt SaP)PtchS)+ p t( ^SXP)=pth t( )* pt(ft. 4 If

21 The leading terms on each side cancel (since, from 2.6.1, P lSP\); utilizing the fact that P P- I results in: ^d pt(S)P+p (AS)^P + P (AP) Since AP involves information we do not have until ^ s A+MA is calculated, we approximate the above by:.t _p t(^5) P Defining = A X ( g ) then: A^ pt( pt)P v - - t with;=5 P. The change inA is approximated by the square of t a scalar, the inner product of g and the eigenvector associated with \. This inner product is simply the sum of those elements of P which correspond to the nonzero elements of i (zero or one). 2.7 Three search illustrations A computer program was written to generate these gradients for a given point X, and the "linearity" approximation investigated for various design configurations. Only for designs which allowed P the total xX to be greater than p did the optimization seem to merit the gradient approach over, say, pure random search. This

22 is obviously a serious deficiency, if true in general, since much of the interest in experimental design lies in problems where the total number of observations must be far less than the total factor-level combinations, p. After generating several series of these computations and investigating several heuristic rules unsuccessfully, work was halted along this line of investigation. Three search histories are next illustrated, one of which resulted in a local optimum solution. In each of the following, the problem to be solved is: Maximize,(X) Subject to ZX; \ X i_ Ointeger (i=-I, ) The resource constraint region was simplified to more easily investigate the properties of )r(X). Since ArX) is nondecreasing in X, the inequality X; N may be replaced by the equality Z X= N. Illustration one. —A simple three-factor, two-level experiment is considered, the type described in (1.4), with the maximum number of observations N, equal to four. The solution to this problem, noted in (1.4), is either Xo=-(\,,0 Io0 \,O)or E-Xo, where E1 is a vector with all elements unity. Here S has maximal rank four,

23 so r = 4. The following rules are defined to increase (and decrease if necessary to insure feasibility) one component at a time, with the minimum step, 1 or -1, for a maximum change in A+: (a) Given a feasible solution, X2 x(, calculate the values of,) such that dj =^ Sjt is maximized. In this example, Sj= if 7x;- +, otherwise 8$ 0. Let the new X have elements i X +O, where rK = t- k _t3wsj svS ZX 4- L+ 0, oihtkeyse (b) Stopping Rule: If )h= N/ A4/., stop. Table 1 illustrates the history for one initial feasible solution. D represents the maximum dlj or predicted change in Am TABLE 1 INTEGER THREE-FACTOR, TWO-LEVEL SEARCH PROBLEM K_________________^ D0 ).... (1,1,0,0,1,0,1,0)................ 0.59 1.2 (6,1) (0,1,0,0,1,1,1,0)........... 0.30 1.2 (4,5) (0,1,0,1,0,1,1,0)............ 0.30 1.8 (1,2) (1,0,0,1,0,1,1,0)..P............. 200a 0.0 Computation stopped with optimal A).

24 Although the optimum was found in the last row, the gradient technique was not responsible since the increase expected in As at each step (D) is seen to be far in excess of what was obtained. Relaxing the integer requirement, and introducing steps of size 0.5 resulted in the history shown in Table 2. TABLE 2 NONINTEGER THREE-FACTOR, TWO-LEVEL SEARCH PROBLEM X As_ D(,j) (j) (1,1,0,0,1,0,1,0).................... 0.59 0.6 (6,1) (0.5,1,0,0,1,0.5,1,0)................. 0.82 0.98 (4,2) (0.5,0.5,0,0.5,1,0.5,1,0)............... 1.0 0.5 (3,5) (0.5,0.5,0.5,0.500.5.55,1,0)............ 1.2 0.8 (8,7) (0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5).......... 2.0 0.0 Again, the optimal solution /^- is obtained, but with noninteger components in X. The estimate D is still much higher than the actual change in As. It was found that "steps" in CXi) had to be reduced to 0.1 or 0.2 in order that Cilf be approximately linear in aX and hence a proposed step direction to maximize a4d be successful. However, as seen in the second case of this example, with no provision to "force" an integer solution, the search converges toward the point (X{ = constant for all i), in this case.

25 Illustration two. —Two fundamental problems associated with the search technique are illustrated by the previous example. These problems also exist in the example in this section, so a different approach is shown. This problem has three factors at 4, 4, and 3 levels. We suppose the number of observations is limited to twelve. A design for this problem was proposed by Webb (51), although no comments as to its properties were given. The proposed solution was Xt 1, for = 1, 9, 11, 16, 20, 24, 27, 29, 31, 38, 42, 46 and X = 0, otherwise (Figure 1). In this problem r= 9, and the value for A? for the above values of the (Xi( is Aq= 2.0. Since 9 ~ ~/ 3 = 3, it might be expected that a better design exists for Z X = 12. It is the opinion of this author after many computer runs that no better design does exist and 2.0 is the maximum value attainable for Xq. Integer search similar to that in the previous example resulted in little success. A different approach was tried as follows: Rather than begin with a point X at the boundary of the feasible region x; = 12, begin with an interior point. Since A9= 0 for many of these points it was decided to take a unit step (when in the interior) which minimized A, the maximum eigenvalue. It The gradient approximation for Ar(X) as developed earlier is not valid for Aq = 0.

(paoZPTOua uolnlos) SUOltWAJ9sqo 9AT9a& qlT&'mtaqoJd gxpx, B Joj (1g) uoFlnTos sqq^ —q t'lTJ 0' 9~x, ZIx (hm Cox I zx 5 seX 8e X 9I f (b) u x e 6~ X ('"x) (X h^ ( <x %X JJ ^x I Ix 1X eh 6 si ja 92x zx Ax'AX 8, si 9~, 2,9^

27 was hoped that by suppressing A, the increase in X; would be more apportioned to the lower-valued eigenvalues. Table 3 is a history of this approach. The initial value of X was (35 = 1, and all other XI = 0. The column labeled 4X indicates at each step which X; was increased from zero to one. TABLE 3 INTEGER SEARCH, MINIMIZING MAXIMUM EIGENVALUE (initial Xs- 1) 4.0 0 X = 1 5.0 0 X = 1 6.0 0 X4^ =1 7.6 0 X2= 1 9.8 0 xz= 1 11.4 0 X2 = 1 13.2 0 X o= 1 14.8 0 X3 = 1 16.7 0 X16= 1 18.7 0 =3.= 1 20.4 0.14 Xo= l1a 22.5 1.32 Computation halted at boundary.

28 Another initial value, XI = 1, was tried, and the value of \q = 1.32 was (with the same approach) also obtained. Since 1.32 was quite a bit lower than 2.0, this approach was for the most part abandoned for others mentioned in section 2.5, and Chapter III. It is noted that gradient calculations for this problem were found valid for a range in step size of 0.2 to 0.25. Illustration three. —The search approach has one area of application, that in which the total number of observations allowable is greater than the number of factor combinations. This is illustrated again merely with several gradient calculations and resultant values of Ar for the 4x4x3 factor-level problem of the last illustration. However, it is assumed that the constraint region is given by 2. 5f(I. Searching on the boundary gave results shown in Table 4. In this example, the initial solution was as follows: X = 5, t = 1, 9, 11, 16, 20, 24, 27, 29, 31, 38, 42, 46, and X = 0 for all other Because of the severe nonlinear character of the minimum eigenvalue function in areas of many experimental arrangement problems, no obvious extension of the gradient search method is evident. Several extensions were attempted, however.

29 TABLE 4 INTEGER SEARCH, TOTAL NUMBER OF OBSERVATIONS LARGE Predicted Step Sugest Change in "9 Change 1 4.4 ), 0 -1 1.8 2 6.2 X, 0 - 1 1.4 X2, 5 — 4 3 7.5 X38, 5- 4 1.0 47, 0 -- 1 4 8.5 ~, 5- 4 0.6 X45 9.1 5 9.0

30 A heuristic rule which allowed noninteger steps (steps which preserved approximate linearity of the objective function), but attempted to weight heavily directions toward integer points, met with little success. In illustration two, section 2.7, one result of the method of beginning the search with an interior point of the resource constraint region was given. This led to a specification of an alternate objective function (the maximum eigenvalue, which was minimized) during the initial search. Although this approach did not yield direct results of any merit, the idea of an alternate but related objective function pointed to a host of possibilities. The result of following one such direction is next discussed in Chapter III.

CHAPTER m A SUM OF SQUARES CRITERION 3.1 Sum of squares criterion defined; rationale The minimum eigenvalue criterion gave rise to solution methods that were far from adequate. An alternate criterion is now proposed which appears similar to the minimum eigenvalue formulation and is much more attractive relative to methods of solution. Let x?} = X be a solution to the following minimum eigenvalue problem: Maximize A r( X) Problem (I) Subject to X E U X ~ 0, integer-valued It is assumed that U is bounded. Let -{1X;3 zxim Xo_ i-,.,p} Consider the following problem: Minimize Z ^(X) Problem P(m) Subject to X E U X 2 0) integer-valued X~? Y,

32 From property (g), section 2.2, we have: (3.1.1) (X) - (X), Xe i where k,, ka are constant provided ZX; is fixed (in this case, equal to rn). Now define, from the optimal solution XK Xf} of Problem (I), the quantity Mo=0X,. Provided I\0 were known, Problem P(Mo) would then seek to maximize an upper bound on A (X) through the inequality (3.1.1) by stipulating that XEY. We are attempting to replace Problem (1) entirely, however, so that Mo will generally not be known (unless of course it is specified in the constraint region). This obstacle is circumvented somewhat by considering a (finite, if U bounded) sequence of problems {P(m);ml.., l. Let the solution to the problem P(M) be X, so that a finite sequence {Xmtis generated. The solution to the problem is then defined to be X2 for which it is true that A.(X )_A,(XJj for all X in the sequence {XJ}. The resultant solution will be called the sum of squares solution, satisfying the sum of squares criterion. That is, if X solves problem P(M) above, then the sum of squares solution, X)< satisfies: nAtX t)= tax n[rnf[A ( XI)

33 It has been impossible to establish that the minimum eigenvalue and sum of squares criterion yield the same solution for an arbitrary feasibility region U, defined by linear inequalities. Such is the case for specific instances, however, and these are discussed in Chapter IV. In addition to inequality (3.1.1), it is pointed out that for Xx/=m7, the minimization of 2At' attempts to "equalize" all the fA,] and hence in some sense maximizes At. The most appealing characteristic of the sum of squares criterion is that algorithms are available for its solution. These will be the topic of the remainder of this chapter. 3.2 A property of the objective function From section 2.2, property (e), we have:? 2 Sly where {s;j3==x(- SZ X<(f). The quantity V (property d of that section) was defined to be the constant number of ones in any of the vectors 1;. Since the elements of each'; are zero or one, elements of'S will also be zero or one. We then have: SS=(2 X<<tt)i: )2, >^ t +t 2 t INK -XQX ~~x'Q x~

34 where q[( ] The basic fact used in this development is trace AB = trace BA K - and is applied to: trace ( K )= trace (; K ) = i 1^ = A balanced model is defined to be one in which all row (or column) sums of q are equal. The additive model is always balanced. In case h-factor interaction terms (nonadditive model, parameters involving h factors) are included, the model is balanced if all h factor interaction terms possible (all factor-level effects taken h at a time) are included in the vector oC. For instance, in a two-factor, two-level factorial with all two-factor interaction terms present, (OC= X ) (/, k2 OD(; ^2, / t92/ ~ /222) so that the four vectors t;3 are: a = (1,1,0,1,0,1,0,0,0).= (1,0,1,1,0,0,1,0,0) = (1,1,0,0,1,0,0,1,0) 4= (1,0,1,0,1,0,0,0,1) and the model is balanced. In this case,

35 16 4 4 1 4 16 1 4 Q - i 4 1 16 4 1 4 4 16 For an example of an unbalanced model, suppose c(=(],(,2 Then, = (1,1,0,1,0,1) 2= (1,0,1,1,0,0) 4= (1,1,0,0,1,0) o (1,0,1,0,1,0) and in this instance, F16 4 4 1 4 9 1 4 Q= 4 1 9 4 1 4 4 9 3.3 The noninteger quadratic program As mentioned previously, the advantage of the quadratic objective function is that problem P(m) of section 3.1 is a quadratic program for which solution methods are available. One approach would be to ignore the integer requirement, with P( )then being

36 stated as (for simplicity, we define the vector X so that the constraint region is specified by linear equalities): Minimize X Q X Subject to A X b x,>0O, /=/-.,P Here, Q is positive semidefinite or positive definite so that X QX is convex and sufficient conditions for an optimal solution may be derived through the Kuhn-Tucker Conditions (32): X is optimal if A X - b, x, >O, /= /,. p, and there exists a row vector v2 with V'>O and 4 of appropriate dimension so that: ZQX*- AtaL- V =O V XJ = = - 2., p Algorithms for the definite and semidefinite cases are known when the integer requirement on the [X/} is ignored (32, 17, 8, among others). An optimal solution is now illustrated for the three-factor, two-level problem of Chapter I (page 6), in which the constraints are defined to be:

37' X +X2 +X5+X6 -)2 X3+XU f X7 f*X8 2 AX-b /y x, tX3* * -2 )r5 + ~$ 4,'7 ),x8 =2. x/ X3 (5 7 W 2 i /2 + Xs v v6 v Xa =2 )(, #K3 SX~' XL2> X,...>o. i =/.,.., e Assuming the additive model as stated in section 1.4, the matrix Q in the quadratic objective function is: 16 9 9 4 9 4 4 1 9 16 4 9 4 9 1 4 9 4 16 9 4 1 9 4 4 9 9 16 1 4 4 9 9 4 4 1 16 9 9 4 4 9 1 4 9 16 4 9 4 1 9 4 9 4 16 9 1 4 4 9 4 9 9 16 The initial Tableau of Dantzig's Method (17) would appear: I tZ X3 X4 X'6 7\ V y l A2 V3 a4 I5 Ab v, VV' v. s vGS V8 o | Q -|At zs b |4 O O Q _______^_____________o ^

38 It may be seen that the following solves the problem: X = (1,0,0,1,0,1,1,0) A - (20, 20, 20, 20, 16, 16) V = (0,0,0,0,0,0,0,0) Since the vectors v and X are orthogonal, and satisfy the conditions in the tableau, they are optimal. Fortunately in this case, the solution was integer-valued, which generally may not be true. This solution is also not unique. Since the above quadratic algorithm does not guarantee the integer property of the solution, we will investigate several methods which do, two algorithms which take advantage of the quadratic objective function and an enumeration method which takes advantage of its nondecreasing characteristic. 3.4 An algorithm of Kunzi and Oettli In a recent paper (37), Kunzi and Oettli developed an algorithm for the integer quadratic program for the case when Q is positive definite: Minimize X QX (I) Subject to AX b ^ ^0^ integer ( ^^-^)Je

39 It is assumed that the feasibility region is bounded away from zero, so that the trivial solution X, iO ( 1/ " P) is not feasible. The form of the proof of this algorithm is changed slightly for present purposes, with the addition of the case when Q may be positive semidefinite. First consider a solution to problem (I) without the integer constraint, not necessarily unique if Q is positive semidefinite (uniqueness is not needed as the authors point out). Let the objective function of the solution have the value f()= tQ y d. An equivalent problem (to be shown) to that of (I) is the following: Minimize 1 (II) Subject to: AX ( xi - O) integer ( - I..) p) rr 2X^K/xtQR where 5 is any point for which X-Q= d, i.e., any point on the curve X )QKd Later, a constructive procedure which generates solutions to (II) is defined, but first it is shown that a solution to (II) also solves (I). This is done in four steps following the next definition. Define, for any point, the point X)= 4,^ where rA.E| d/XtQXa Note that,Q <d, so that Xqlies on the

40 intersection of X QX"X and a line through the origin and *A (Figure 2). Lemma 1: Given Xa, the point Xa solves: Maximize 2X'QX (Qpositive definite or positive semidefinite) Subject to X Q X = X >0 Proof: (i) Let L(X )/ - XQX-(;(QX-) be defined over the convex set (X 0 _O ^ 0). (ii) L( X, 4) is the sum of three concave functions, 2 X)cX XVt(-q) X,, and 4d, hence L (XK ) is concave. Therefore to illustrate a global maximum, it is sufficient to show that V ~ L(X*/)=G (32) ) (for feasible X). (iii) Let X = — /a,W/. Then, 7X L(,X /ro) = 20X, -(2/~y(Q(-0 The Lemma is true for positive definite or positive semidefinite Q. Theorem 1: If X,,solves (I), then (XaTOT) solves (II) where rr= o Q X / Proof:

41 x3 / f~~~~~~~~~~~~~,J~~~~~~~~ J~~~~~~~~~~~~~~~~~~~~ // xtJx: /~~~~~~~~~~~~~~~~~~~~~~ /~~~~~~~~~~00 Fig. 2. —The quadlratic program surface. X~~~~~~~^ I Fig.~^ 2lTh udai pormsrae

42 (i) Let YCY {X X b,^ integer; i1, P (ii) By definition, TrT = 2 Xr Xa Q IX Q X (iii) From Lemma 1, 1T 2a 2 QX/d, for any X satisfying t;X (iv) Suppose (Tr,^ X) do not solve (II). Then there exists t am some X for which: TTT2 b (v) By definition of TT, we have: t t eta= - 2X > 2,^ t, X Since this implies that X kQ)oX^QX^ a contradiction of the minimum property of X in the hypothesis (X (fQ < X QX,) the proof is complete. Lemma 2: Let (Xb77) solve (11). Then T1b=2/rb Proof: (i) Since (Xb, r) solve (II), dJ 7b> X X for any X satisfying ) Q X =. (ii) Since Tb is minimum, and2 t"' is maximum at,;)Q (lemma 1) it must be true that:d-bT=2XQ X2/' so that Tl _ 2/ rb

43 Theorem 2: If (Xb)Tb) solves (H), then Xb solves (I). Proof: (1) Let a solution to (I) be Xa. Then 1TF lTb since (Xa, T,) solves (. ~b(ii) ( ) ==(Xo4 so that Xb solves ([). A constructive procedure for a solution to (D), and hence (I), is then described as proceeding in stages, the k stage represented by the mixed-integer program: Minimize 1T Problem S() Subject to A X b X O) integer-valued -r 2 Xi QX/ i- K. Let the solution to this stage be: ( XK)IT). Suppose, for the first' t time, XK also satisfies: T't d _ 2 Xw(QKX.We then have, by the following lemma, that (Xtrw<)solves (II) and hence X is a solution to (I). Lemma 3: Let the sequence 1S(K) be defined as above in the constructive procedure, with IC the smallest integer such that the solution (XK^)TI) to problem S(K) also satisfies T> 2 tQ X,.

44 Then ( X<KT) solves problem (II). Proof: (i) Suppose that (XK)1TT) as defined in the hypothesis does not solve (I), but that (X9)1g)does. (ii) Then fT< f7K, and further lT is not a member of the solutions P'r 1 J.T2.,1TTto the problems S(i)j S (2z)). -, 5(), since d Tr 2 X QXand this property was true only for T, in the sequence P. (iii) There exists a "v" such that 1TV,, T +,, 6 P and rv w < 7V+, (7j''TV J4t (iv) By lemma 1, we have ciT =ZX QX>ZXX for every X that satisfies X so that ( X, ) solves 5(V&)0 rather than ( Xv+1 7j,/ ). This contradiction completes the proof. An example of the use of this method is shown by a simple two-factor, two-level additive model represented in tabular form: FACTo0R 2. LEV EL LEVEL2. FACToR I 3 LEVEL2-. k. ~ ^ ^H

45 The constraints are assumed to be: XI+X 3 -2 AX=- X3 * -4-2 c < c, (,+ ZXf+3)X3+4)44I The matrix ZQ is: 18 8 8 2 8 18 2 8 ZQ 8 2 18 8 2 8 8 18 A solution to the noninteger quadratic program: Minimize X X Subject to A X = C X Co X > 0 is ('/2,3/2,'3/2)which can be illustrated using methods of section 3.3. The first step in the Kunzi-Oettli algorithm, given this noninteger solution, is then to solve:

46 Minimize 17 Subject to A X b CtX < Co 80 r~ 2 28 x, + 44x28 x~.4 X Z O) integer-valued. A solution to this problem is 1Tr2, X (,O,) I), with Xl-X- T. Since 1T,(d/2)=8o is less than tQX = 82 g — o 8A9, we proceed to problem S(2): Minimize Tf Subject to AX= C^<<Co 8on- 28x, + 44 x2.28 x3 444 8o0 X/ (26x,* XZ 7'Dx3+4ZX4) ~ 0 integer-valued. A solution to problem 5() is XZ (0 2, / 2),42 02S. Further, 7T'Z(d/z)-"8O9'Z CtZ-O. satisfying the requirement that 2 is optimal. 3.5 Witzgall's parabolic constraint algorithm An algorithm which also handles the ase whe c hn g is positive definite or positive semidefinite is reported by Witzgall (55). Here,

47 problem (I) in section 5.2 is transformed to: Minimize Z Subject to Z XtQ X AX~ b X>O), integer (1=\,2,. ) with the objective function linear and the constraint region nonlinear. The technique is shown to resemble the integer linear programming algorithm of Gomory, except that the nonlinear constraints are handled slightly differently. The method requires that the function X X be transformed into: o- o(X-l()-...- (), where each L (X)is linear in X and all are linearly independent of Lo(X). This transformation could be unattractive as evidenced by the following two-factor, twolevel problem. Here, -= (1,1,0,1,0), = (1,0,1,1,0), o (1,1,0,0,1), (1,0,1,01), so that:

48 K X2B S tg ^X, +X.'3 +X4 X.X3 X2+-(4 X,Xz ^x,+z 0 x X ~ 5= 1 23 X4 0 X3+4 X3 do I X,+X3 Xl X3 +X3 0 X2+ X4 )y2'W 0 X4b For this problem, the Q matrix is then: 9 4 4 1 4 9 1 4 4 1 9 4 1 4 4 9 and XtQX is equal to q(x, +x+x2+x 2 I)+8(x +yw243)+2(tV S) Finally, "completing the square" as Witzgall suggests, results in: XQX = 9 L(X) + (65/9) (X) + (Q/b) LO(X) + (c/d)L4(() where: a 4-144 b 585 C _405,076 d= \ 515 L(X) = xl + (4/)(X2 4 x X4/4) L2(X) = X2+(9/65) X3 + ( 36/65) X4 L3(X)= X3 + (585/1036) X L,(X)= X4

49 This transformation to the L; (not necessarily unique) could also be accomplished by computing the eigenvectors of Q. The inner product of these vectors with X would provide the set{L;. 3.6 A partial enumeration method of Lawler and Bell; an example The two previous algorithms utilized the quadratic character of the objective function in their solutions. An enumeration technique of Lawler and Bell (38) is now used which does not require such a form, but takes advantage of the monotone property of X(QX,XO, This algorithm could also be used for the eigenvalue formulation since that objective was also monotone in X. Because of the calculations involved, the quadratic form will be used to illustrate this method. The problem stated in the previous section is stated with a feasibility region: Minimize 2 5Z = o(X) Subject to ( X)= X 2 4= b, 92(X)= x3 + X42 = b 33Ox) x, X3 <3 - b3 g(X)= K. Z + 3 * 4-= m= bM 9s(0 = X, + zxz+ 5X3 3Y43 \43(X) = azl +< —3 =b6 x >o), In,,ter(= 1,2z,3,4)

50 Since 6 is the maximum value of rn consistent with the constraints, this is the value that will be assumed (not necessarily optimal over all m). Further, with this value of hi, the inequalities become equalities (with the exception of the cost constraint). The objective function may be reduced so that the problem is stated: 2.. 2. 2. Minimize Xt+ X X3+ x + (constant, ignored here) Subject to ( 4 9(X) = 3 S9CX) 6 95(x <_ 4 Xi C0 integer (f = 2, 3,4) Since X1_4) X24, X3 2 ) X C2., following (38) we define: 1 \ Zo+ aZZ + 2Z1L ) X24 z + Za +22 z X3 Z30+ 2Z31 X4- Z40 2 wert Z 0o1 The algorithm is then as follows (Let Z- Zi.4)Z,Z3,Z3oZ2)22,Z7,ZOZ/Z/LZ/o):

51 The vector Z is considered in increasing numerical order. Define Z to be the first vector following Z in that order which has the property that Z, element by element, is not less than or equal to Z Let gZ)be the minimum value of o(Z) to date. The following rules are defined indicating which vectors may be "skipped" over (note the j(X) are nondecreasing in X orZ): Rule 1: If So(Z) co(ZO skip to Z since no vector between Z and Z can have a lesser value than 90(r Rule 2: If Z is feasible, skip to Z since no vector between Z and Z can have a lesser value of the function o than o (Z) Rule 3: If g;{ZY l)<t;,{5;, skip to Z since if 9;(Z- I)4b; no vector less than Z 1 but greater than Z can satisfy (Z - b 5, since g is nondecreasing. Rule 4: If 9(Z)>b, skip to Z; since q is nondecreasing (t0). Figure 3 is a partial history of the solution to this problem. The optimal solution is X = (2,2,1,1), with 9= 10. At this point, the numerical "value" of Z is obtained by identifying Z with the integer h (Z)= z'2 a. *8. z,0o2?o Z.2z,, 2z ~. In the rules 1-4,Z refers to this binary number.

52 X4 X3 X2 XI Z41 Z40 Z31 Z30 Z22 Z2, Z20 Z,2 Z, Zo RULE INFEASIBLE? TO Z? O O O 0 0 0 0 0 0 1 3 X O O O O O O O O 1 O 3 X O O O O O O 0 1 0 0 3 X O O O O O 0 1 O O 0 3 X O O O O O 1 0 O O O 3 X O O O O 1 0 O O O O 3 X 0. O 1 O O O O 0 3 X O 0 1 0 O O 0 1 0 O 4 X O 0 1 0 0 0 1 O O O X O O 1 o 0 0 1 0 0 1 X O 0 1 0 0 0 1 0 1 0 4X O 0 1 0 0 0 1 1 0 0 4 0 0 1 0 1 0 0 1 0 2 130\ X 0 1 0 1 0 1 0 1 0 0 4 X 0 1 0 1 1 0 0 0 0 0 4 0 1 1 0 0 0 0 0 0 0 4 1 0 0 0 0 0 0 0 0 0 X 1 0 O O O O O O 1 X 1 0 0 0 1 O 0 1 X 1 0 O O O 0 0O 0 O 1 0 0 O O O O 1 1 g311 X 1 l 0 lo0 l0 1 1 0 0 4 - I 40 0 0 a ~ I I, Fig. 3. —Illustration of the enumeration technique.

53 3.7 Use as a starting point for the minimum eigenvalue search problem In section 3.1 it was pointed out that a relation does exist between the minimum eigenvalue and the sum of squares criterion. It has not been proved that both criteria yield identical results in all cases, however. The solution, X, to the minimum sum of squares problem could provide an initial starting point for the minimum eigenvalue gradient search procedure, probably a very good one. If it is true for X that A(X) is very close to an upper bound (say, the one provided in section 2.2), then X might be considered adequate for practical purposes. In this chapter, an alternate formulation to the problem of Chapter I was accomplished. This alternate approach was shown to be (a) computationally attractive, (b) heuristically related to the Chapter II formulation, and (c) related to improving an upper bound for the Chapter nI objective function, with some constraint on the feasibility region. The two criteria are next shown to be more closely related when, in the next chapter, the constraint region is further restricted to one for which constant cost of experimentation (all c; equal) exists, and the total number of observations is specified.

CHAPTER IV FURTHER EXAMPLES OF THE SUM OF SQUARES CRITERION 4.1 The equal-cost problem The sum of squares is further substantiated as a reasonable objective in this chapter for the case when all C; (cost or time per unit increase in Xi) are equal so that the constraint CX<co may be replaced by 2Xi (. It will be shown that in three important cases this approach yields precisely the solution that would have been derived by the minimum eigenvalue criterion. In a fourth example, the arrangement proposed by Webb (51) is also derived. To facilitate these results, three observations are now made, one which describes the optimal solution of a simple noninteger quadratic program, and the other two, sufficient conditions for optimality of the integer quadratic program for the additive linear model. 4.2 Three preliminary observations Observation one. — )X t-(lil..,' I) for scalar k)O not necessarily integer, solves: 54

55 Minimize (X) Subject to X; - N X > where the model from which Q is derived is balanced (section 3.2). (i) Define E to be a (pX ) vector with all components unity, and let F-f(X) +T(ZX;- Since Q is balanced, let Q;~E for some scalari. (ii) Since XQX is convex, a sufficient condition for optimality of nonnegative X is the existence of T* (and let:s /p to insure feasibility) so that: 7VX (X",*')=oor 2QX + T*EO (iii) Now X(= (/p), and t= -2 p satisfies the above condition, since c;i' } for ( i, p, true for a balanced model. Observation two. —Assume the three-factor model, E(Y)ZZO, or relative to an element ijK E b ) E(jlK)-/A+T;+ Sj4~ fAJ:-)..nf The extension of the following to an arbitrary number of factors is obvious. In this example, the total number of factor-level combinations

56 is p-n 2n. It is desired to establish sufficient conditions for X to solve: Minimize (4.1.1) j Subject to j XK' K iiK X 2 kOj integer for all i ^ j where it is also assumed that N p and XiLj denotes the number of observations to be taken of the random variable ij.' The (e) subscript denotes summation over that index, so that j. = X, X.oK= X K',etc. For a given X=(f}, the matrix S is represented: SC W it I D D 3 t3 Y43 DIS Di3 D3 where: \'J( xn)=-(xt.. X2I,,... )r,..) W2(I X=(X.. )X^.,. " ~ ) W3 (IX t).x2.J), 6..) -n.) DI (nxn,)=dg(x1. xz.). *..' xN,..) Dt (n~,x n3)'~ (X *"- ) x *X..*-,-.

57 { 3(n n3)= {X; j-.~ ^ D23(nrxn {X-j S From observation one, it is known that the optimal solution without the integer constraint is CXxI=N/p(I, 1I..,1) It is desired to find sufficient conditions in order that an integer-valued X= xl.] also be a solution to this problem. These conditions are obviously: (4.1.2) X.k r.i-n N/P \t i X.;n 0 n%)/e-) n3J X;jK 20, n er ) Observation three. —Finally, we suppose that such an optimal integer-valued X does not exist which satisfies conditions (4.1.2). Let: Tixjt=enh 3 opn$j, th er n 1on. I2 5 {( t)j); ( ij)Ion iSj- )z,...I+ni+ntrnoll The matcix S then has components S.j, with(i'j) either in I1 or I2.

58 With (X).=.2 and h(X)=-S2., the problem (4.1.1) now becomes: Minimize 9(X)+h(X) Subject to ZX JN XiK >Ointeger We seek sufficient conditions for a vector )X to simultaneously minimize 9(X) and h(X), and hence minimize the sum. The function ^(X) may be written: h\ (X) = h2 3 { X0, + Xj.i + 1 xe A sufficient condition for X =X to minimize h(X) is that the components for each of the three sums be equal, i.e.: X.o= N/ (4.1.3) N/ X..- N n3 Since g(:-2 X{.+.~+ ZXKZ, a sufficient condition that 9(X) be minimized subject to the same conditions is that:

59 axj - WVnn -f3N/ p (4.1.4) nO/) VOJ % msete,r Suppose 4A/nInt)N/nAI or NItnl^ are noninteger, less than one. Then no integer-valued XO can satisfy conditions (4.1.4). By taking any one of the three elements of 9((X)-zX-iZ.+ZXI )-X(, the problem becomes, say: Minimize! Xj. I.j Subject to X. =, to J integer Since we are requiring the sum of squares of nfl, nonnegative integer elements to be minimized while constraining their sum to a number less than nfli, an obvious sufficient condition for optimality is: (4.1.6) I ILK l K O, i Here, the last two equations are obtained in a manner exactly as the first.

60 These three observations are now used in further illustrations of the sum of squares criterion. For brevity, a solution using the sum of squares criterion will be called L-optimal if it is a member of the set of solutions to the minimum eigenvalue criterion with the same constraint region. 4.3 The one-factor problem The problem to be solved according to the second criterion is (we may assume each X/ is equal to an eigenvalue of S in this case): Minimize A-= x i Subject to X< I ( ).., p) i 0, integer ( i, -,p) It is assumed that qll, 4 Cot... 4 qp, and that a feasible solution exists, i.e., l,; >h. It is desired to show that for a solution to this problem X = iX}, it is true that X' = Minimum U ['m/pl}, )b] denoting the largest integer less than or equal to b. In this case the Minimum XX has attained its largest possible value, and hence X is L-optimal.

61 Let X Ix-xij\;X ) Xj integerJ i,...p}. If X (rE, then Xi is L-optimal since the value of minimum {Xt } cannot be increased subject to L x! = m. Suppose now that Xi?7, and that X,< Q. Two cases are considered (let = maximum {X. ). Case 1: X< f. Then there exists a "i" such that Xj > X+2 (since X',). Let X=; be defined as follows: Xl o XA I. w x<K) I < ) <3.(Kfeasible) By optimality of X, (xi()2,x 2(,1)2 I-Z+ (XuI (x in or kj i(K+ 2 <X j^^ -XJ) 2I so that X -X <, a contradiction. Hence X 6 F or X i r but X 4i1. In either case, X is L-optimal. Case 2: X - f. Here there must exist XJ such that )X Z X. +2. Interchanging the role of X and Xj in case 1, the argument is identical, noting

62 that the X in this case is feasible since a % as, Kz, p- 4.4 The latin square The notation of section 4.2 will be used with the model E(Y) 7being written: E (;) =/L t=;+ s'+ T(ijK i-,d'). The number of observations of the random variable YI, is again given by X;j, The problem is stated: Minimize eS ij J Subject to X cj K =ad )(jK>Ointeger ({j', l.., a) In section 4.2, we have the sufficient conditions (4.1.2) for XIK - X to solve this problem (here p: 3 n,,n:-n.id, v' 0 X.' X?. \l \ ij:,i (the three implying) XjK I orO X

63 These conditions precisely define a latin square arrangement, and since Ehrenfeld (19) has shown that /A(X)is maximized for this factor arrangement, X, this solution is L-optimal. 4.5 The balanced incomplete block The model for this problem is: (Yj) ^ k+; + ( i-,,.. t; j).6 ) The element ij is the number of observations to be taken with factor one (usually called treatments) at level i, and factor two (blocks) at level j. The problem is written: Minimize S.j Subject to X.. < bt Xij O,integer ( = W~,t j -,-,b) It is supposed that a balanced incomplete block exists so that N b and N/t are integers. It is now desired to examine the sufficient conditions for XO to be a solution. We utilize the results of section 4.2, observation three, specialized to this two-factor case. Since Nlbt< I (here nI=b,n2-t in notation of section 4.2), sufficient conditions for optimality of VXX^ =X are obtained from (4.1.3) and (4.1.5):

64 X,. =N/b (assumed integer), i = lI, t X\ N /t (assumed integer), j I.i These conditions define a balanced incomplete block and hence X0 is L-optimal since Ehrenfeld (19) has shown that A/(X)is maximized with such an arrangement. 4.6 Webb's problem Webb (51) reported a factor arrangement for a three-factor problem with 4, 4, and 3 levels. This problem was referred to in section 2.7 (Figure 1), and it was noted there that the gradient search technique failed to obtain the design proposed by Webb. The number of observations was limited to twelve (here p-n, nzn3=4) or one-fourth the number of factor-level combinations. We shall see that the arrangement proposed by Webb (the encircled elements of Figure 1) satisfies the sufficient conditions for optimality relative to the sum of squares criterion. The model is a three-factor additive one so that the notation of observation three, section 4.2, may be used where now n,=4-, nr^4, N=33, and NZ —. The minimum sum of squares problem is:

65 Minimize 2 S Subject to X.. =12 Xi'j'kOinteger i,:,..4; j., 2, 3 Let X be an optimal solution to this problem. Since N/1,3= I and NIJ^^rtl, (4.1.5) yields the conditions: I,, \ hJ 60 X' J 1K^1,2, 3 Since N/nl, t.<, 4.1.6 yields: x~j.=o, \ Property (4.1.3) finally yields: X-. =3, -,, 4 X~,- j-. t. 2,3 We now observe that Webb's arrangement clearly satisfies all these criteria as shown by the matrix 3:

66 12 3 3 3 13 3 3 3 4 4 4 3 3 0 0 0 1 0 1 1 1 1 1 3 0 3 0 0 0 1 1 1 3 0 3 0 1 1 O 1 3 0 00 3 1 1 0 1 1 1 3 1 0 1 1 3 0 0 0 1 3 1 1 1 3 0 0 1 3 0 1 1 0 00 3 1 3 1 1 0 1 0 0 0 3 1 1 311 4 3 1 1 1 1 1 4 4 1 I 1 11 1 1 110 4 0 411 1 I 1 1 1 1 0 4

CHAPTER V SUMMARY AND AREAS OF FUTURE INVESTIGATION 5.1 Summary In structuring the problem of experimental factor arrangement as a formal mathematical optimization problem, two objective functions were used. The first, a minimum eigenvalue criterion, was a conservative approach described possibly as a MAX-MIN method in that the minimum informational content (over normalized linear combinations of parameters to be estimated) was maximized. This criterion led to a multidimensional lattice search problem, useful only when the method of "rounding off" noninteger solutions could be justified allowing conventional continuous parameter search methods to be used. It was empirically determined that such approaches might be valid for problems which allowed the total number of observations to be at least as great as the number of factor combinations. An alternate criterion was introduced which was related to the minimum eigenvalue criterion. The second objective function, 67

68 the sum of squares of the eigenvalues of the information matrix, was seen to yield identical solutions to the minimum eigenvalue problem in three important cases: (a) the one-factor problem, (b) the latin square, and (c) the incomplete block arrangement. Furthermore, this second criterion yielded a quadratic form in the decision variable so that existing algorithms for the integer quadratic program could be utilized. Computational feasibility of both criteria is marginal (Appendix C), with problems on the order of five factors and four to eight levels per factor probably requiring minutes to hours of computation time on the 7090. Investigative effort in programming efficiency could possibly cut these times by a factor of ten. Originally, work on this dissertation began with a costminimization problem. That is, with the environmental constraints and a lower bound for information, the objective was to seek a minimum cost factor arrangement. With information defined as the minimum positive eigenvalue, t, the problem became an exceedingly difficult one since the constraint region was no longer defined by linear relations. It should be noted that with the sum of squares definition of information, this leads to a quadratic form in the constraints, the problem that Witzgall's algorithm (55) does handle. It was eventually decided, however, that probably most sizable experimental factor arrangement problems have a fixed upper bound on

69 total cost, and it is desired to find the arrangement of maximum information content; hence the formulation of Chapter I. 5.2 Areas of future investigation The areas most suggestive for future effort are most easily categorized as (a) theoretical and (b) operational. Both are now discussed separately. Theoretical areas of investigation. —One extension of this dissertation would be to include random effects models and mixed models. Another would be to alter the conservative objective function by considering the minimization of a particular variance with other variances being constrained. Effects of "decomposing" the larger problems into smaller ones and heuristically "tying" them together might prove fruitful. There remains the question as to what sufficient conditions on the constraint region must be in order that the quadratic formulation yield the same solutions as the minimum eigenvalue problem. Finally, the question of improving existing partial enumeration techniques for this specific problem has possibility. Operational questions. —The first operational question to be answered would seem to be whether "typical" problems which involve

70 constraints of the type mentioned can really be solved. A particularly promising area could be that of marketing. "The high cost of retail experimentation discourages the consideration of more than three or four alternatives" (25). Because the determination of a feasible solution involves existing computer programs, the question of how "close" such methods do get to optimal solutions would be very useful information. Finally, due consideration should be given to the vast existing knowledge in experimental layout. If the constraints readily admit a factor arrangement of known optimal properties (presuming such information could be stored in a computer, such as characteristics for an incomplete block), then a computer program should test for such a condition. Such investigations could obviously become quite exhaustive.

APPENDIXES 71

72 APPENDIX A EIGENVALUE PROPERTIES Below we prove properties b, d, and f of section 2.2, Chapter H. Property b: Ar(X) Minimum { 5;;! (i) Let AI ~. — " Ay' be the nonzero characteristic roots of the real, symmetric (rxmv) matrix S of rank r. Let R(mwn) be such that RtS R: pd0 g A;} 2, RtR - I (ii) Then'IAr- Ma ie SX:O XX^ X ioe, since, X K 5X for X' Y, Y { Y;Y Xt =I Ye Y C (iii) Therefore, letting X; be a vector of all zeros except unity in the jti position, it must be true that Ar xsX i = So -or all ior I minlOwts

73 Property d: 2 A; = V Z),, where V is the (constant) number of nonzero (one) entries in any'. Recall, S=. X g 3. where elements of T/are 0 or 1 (and hence so are./.t) and each first component of l is 1. We have 2A; =trace S = Ix; race,t^7 = Zx; [ race S = i X; Property f: For the additive model with g factors, each factor having n; levels, Ar(X) /ii, where ~X;I, and t is the maximum of the V;. The additive model introduces a constraint on the S; to the effect that the components of 5{ relating to factor K ( nK in all, one component for each of the nK levels) may have only one nonzero entry. The singular S model assumes all the S; have "1" in their first position. Then s5; represents the total number of times parameter ""' (pertaining to a level of one of the factors) is observed. With Zx{=-J, we have: n,I =; 2. I::n~+2 t+c.

74 Therefore, with the definitions Ei = maximum {nrj, S = minimum s5 1, j I the largest value that S can assume is N/n. From property (b) above, Ar(X)< S, so that the quantity N/A represents the maximum value attainable for ^A(X)with fixed N.

75 APPENDIX B OBTAINING A FEASIBLE SOLUTION B.1 General approach If the gradient procedure of Chapter I is used, an initial feasible solution is needed. Further, every attempt to minimize costly computation of eigenvalues of that method should be investigated by looking for less costly methods of obtaining "good" initial feasible solutions. It is also possible that such a solution would provide a better algorithm for the integer quadratic programming problem of Chapter I. Additional reasons for such investigations are: (a) This approach could yield an optimal solution in itself. (b) Smaller problems might be solved heuristically. (c) Greater insight into the entire problem might be gained. The general approach will be to ignore certain constraints so as to take advantage of the resulting problem structure, and if the solution satisfies the previously ignored constraints, a feasible solution has been found. By judiciously choosing the constraints to ignore, this strategy could prove successful. The problem of obtaining a feasible solution is to find an X, so that the following four constraints are satisfied:

76 (a) A X0 b,a j o, I (i= s,,); j) X~ > (b) CX <Co) C O ('j-'> l )0 ~ > (c) X4 integer, ( — ) P) (d) Ar((X) >0 As in Chapter II, requirement (d) is added. In the general problem of the maximation of /r(X), the requirement that this function be positive is superfluous and this restriction further introduces the undesirable property that the feasibility region is an open, rather than closed, set. For purposes of this section, however, the requirement that Ar be positive has the advantage that should an initial solution actually be used, A>0o insures the maximum number of estimable linear combinations of the unknown parameters. Recall that the matrix A has elements zero or one and each row of A X0 corresponds to the total number of observations available for a particular factor level. Requirement (d) is satisfied only if the space generated by the vectors {g X? >0} is of dimension r, i.e., there must be r linearly independent vectors corresponding to the nonzero Xi. The problem of obtaining a feasible solution is approached by solving Problem I (Figure 4), or, if the solution to that problem

77 Feasible Solution Problem: Find X so that (1) AX b, X > 0 (2)CtX < Co, X (3) X integer valued (4) r (X) > 0 Let = set of feasible solutions Problem I Min CtX Subject Xi = 0,1 Xr(X) >0 Solution = XA ctA>C? Yes Empy Axs xA AXA b? Yes XA ^ ^p No Problem IT Maximize V (Ai denotes subject v Ai X bi, for all i ith row of A) CtX < Co X>O, integer valued Solurion = X, Vb Ar(XB)>O? >Yes{1 XE / No INDETERMINATE Fig. 4. —Linear model.

78 produces indefinite results, Problem II is then solved (Figure 4). Problem I essentially seeks a minimum cost factor arrangement subject to r>)O. If, for the solution, the (minimum) cost does not exceed the allowable upper bound, co, but the resource constraints AX<bare violated, Problem II is solved. This latter problem essentially ignores the constraint /)rO, and seeks to maintain the cost and resource constraints. The function Or is not completely ignored in Problem II, however, since there is a set ~K~ such that S'i = AX ~j I K for every diagonal element Ai;. Here, AL denotes the LAtrow of A. Since V is maximized in Problem II, and Ar(X) 4 minimum {S;' (property b, section 2.2), Problem II maximizes this upper bound of /.(X). This procedure does not guarantee that ^ (X0 )> 0 for a solution, X to Problem I, and the result of solving the two problems may still be an indeterminate situation. The problems are now examined in more detail. B.2 Solution of Problem I, linear model Ignoring constraint (a), the following procedure is defined to obtain a "least cost" set of r linearly independent 5i which in turn define a set of least cost observations satisfying (d): Proceed by stages, considering vectors A, in order of nondecreasing cost,

79 Ct. At each stage choose the vector if it is linearly independent of the previously chosen ones; otherwise discard it from future consideration. This process is ended when r vectors are chosen. Corresponding to a chosen S, set X; equal to one. Otherwise A set X equal to zero. Let the resultant vector be denoted by X As a result of this process, three cases are possible. If C X> Co, the problem has no solution since C tX is the minimum value of C X subject to /\r()>. The proof, due to E. Lawler, is found in section B.7. If C:X co, and for X,, X )it is true that Z aCj Xg bi t(,..) p, then X is an initial feasible solution. Otherwise one or more resource constraints are violated. This last condition is now considered. B.3 Problem II, linear model The problem is altered slightly to accommodate the last case above. Find X O,y j Oj 0 so that: jx.+y-i (+ -y I X') CX -+VI " -Co A,(x)> o The third case in the previous example corresponds to a set ^ 20 2 Nw 7 0,but for which one or more Y are negative.

80 Ignoring the constraint Ar(X)>O, one could seek to "drive" those negative elements of ys toward positive values while maintaining the already feasible elements by solving the problem (equivalent to Problem II): Maximize V Subject to V\ < ( il.-e\)'a; jxi y = b; (i-l,.,) Xjy;,w 7~ (ill)-^^, -j v;p) If, for a solution X, it is true that /r(X )>Q, X is a feasible solution. Otherwise, we have an indeterminate situation. We now consider the special case of an additive model, with the general approach depicted by Figure 5 and discussed in the next section. B.4 The additive linear model; analogy to the transportation problem At this point, for reasons presently obvious, we briefly consider the characteristics of the solution to the generalized transportation problem of linear programming. The primal and dual problems (28) consist in finding a solution (X (~) so that the following conditions hold:

81 Primal Dual CtX is a minimum bt U is a maximum AX~= b AtUO Co X~>o The matrix A has a special structure, and the equations AX=t) are illustrated for the case when elements of X have three subscripts: Z, - b(~) IjK z X b(J) (j t n ) L_ < = i (i-,,~, kc) Z,:., nK) ijK " (I-$' In the transportation analogy, one would have three shipping "stages" labeled,JJ K, each with n Z., Vn3 substages respectively. Then XiKi refers to the amount of goods shipped from the (t- substage of I, through the jet substage of J, to the K% substage of K. An optimal solution to the noninteger problem (the solution is automatically integer-valued for two shipping stages, but for no greater such number) which minimizes CtX consists of a "basic" set {Xi} of nonzero elements (the remaining being zero) and further, the columns of A corresponding to the basic elements are linearly independent. The vectors X) and U are orthogonal, i.e., U XO,

82 and for XiJ r0, the corresponding value of UjK inserted in Cj<Ci i-Ui gives the marginal cost CiI in introducing a unit increase in the variable X ib (zero, if XK is in the basic set). The additive experimental design model E(Y)=Zo( is considered, and assume for simplicity that a three-factor model exists so that X j must be determined for i.:,n,3j,,,^^),\& =L)) * ), 3 3 to) (0 The vectors ~ (or rows of Z) are of the form e'+ e) e se,, where the vector e is defined with n, +n2r + 3 components all zero except component o+b. It may be seen that every t corresponds to exactly one column of A in the transportation problem structure so that specification of a basic solution Xo, and a linearly independent set of columns of A then determines a linearly independent set of r (here rr:n,+nI.* n- I ) vectors, t }. The above analogy for the additive model enables a different, possibly faster and better approach to the problem of finding a "good" initial feasible solution. B.5 An approach to the three-factor additive model The problem now considered is again for three factors, the generalization to nY factors being obvious. With the three factors at nl, nr, and n levels (or with this number of parameters associated with each, plus one constant parameter) let X;j again

83 represent the number of observations taken with factor one at level i, factor two at level j, and factor three at level K. Now define: tP. {X; (X) I> 2 e { X; AX< b, XO Bnd^ ) = {EY AXr A - dfr 0d <a b< X2 0 tL4.-{XOClt< NOTE: C3(C2 The inequality AX ~ b is defined in further detail:. c)b(,) { —<blb ) (i -...n,) x. <b) (j - I, n) The subscript ) again indicates a summation over that index, i.e., Xa = X F XKj = X, KS etc. Two approaches are now considered, one preceding the other. They attempt to find a "good" initial feasible solution, taking advantage of existing linear programming techniques for their solution. Recall that i)(X) < Minimum {5i i-. Here, V = minimum i1..) X.j, x... An attempt to obtain an initially good solution is to make V as large as possible by solving the following integer

84 linear program (Problem II, Figure 5): Find X, so as to Maximize V Subject to V <X ( I..In,) v Xj. (j -,., 3) X integer valued If, for a solution (v, X)) it is true that Vo 0, the original problem has no solution since it follows that r() = 0. If V0 O, and Ar(Xo) >, then Xo is the initial solution sought. Otherwise, the next problem is suggested. Let the elements of the previous X be denoted by xy^, and define a vector d with elements: (0) o d(') =x?.. (- =.. _^j. j =.,) )-.K -,.K.. Now consider the problem: Minimize CtX Subject to X 3(d) With X integer-vaiuea

85 Maximize V Problem IL Subject VAi X bi, for all i Ct X < Co X >0, integer Solution= X8,VB V=O 07 Yes -- s Empty No Problem In Minimize CtX Problem ~ | Subject AX = AX X > O, integer Solutions = XC(note:XCt Co) (Any) Solutions XC non-degenerate Yes XC No >No > INDETERMINATE Cor all i? Yes INDETERMINATE Fig. 5. —A dditive model. enforce non-degeneracy so s3 xDe INDETERMINATE Fig. 5.-Additive model.

86 It is observed that since )XoE w and Vo >, then O<d<, a requirement that (d) be nonempty. This problem is precisely the three-dimensional transportation problem for which a procedure for approximating an integer solution exists (Dwyer, 18). If, for the solution (with the integer approximating technique) Ct O>C, this approach has been fruitless. If CXo <, and Ar(Xo) >, an initial solution has been found. Finally, if C Xo C, but /rXo)-0, some of the elements of the basic solution in Xo are zero, or one has a degenerate solution (otherwise, the columns of A corresponding to nonzero X; would span a space of dimension r). Let these zero elements of the basic solution be the set J. One possible next step would be to increase each element of J to one, and if the resultant ^ y 4, it would be the feasible solution (clearly R't ). Of course, it may be possible that even though a degenerate solution is obtained, one finds multiple solutions one of which is nondegenerate. Necessary conditions for such a situation could be easily developed. B.6 The two-factor additive model: illustration It is now assumed a two-factor problem ended in the first phase with \/o, but with the corresponding Xo, Ar(Xo)( 0. The second problem above would then be stated:

87 Minimize CtX Subject to X.. d( x. =d ) (i= -...r,) - da^ (J0 I- ^,n2) X'j- J It is known that solutions to this problem, the two-dimensional transportation problem, are integer-valued. The solution, Xy, is assumed to satisfy CX = Co- o, for /o >, but that r ( Xo) 0. The latter means that X) is a degenerate solution, but the cost constraint to the original problem is inactive. We seek a solution Xi, such that C0 C X, C for which the "row" and "column" sums, Xl and <, remain unchanged. The matrix C'{ } is defined as the dual solution of )X (see for instance, 49, 31). It represents the marginal cost of introducing an element into the solution which is not in the present optimal solution. The vector X= fXtZ } has tn,+ ol members of the basic solution, the indices of which may be identified by elements of the "tree": K= [(,jl)) (i2,j )) (:,1j~)'( [3IJ-).".)(ir) ij)] where one or more fXij; (ij)E K are zero. Notice that consecutive elements of K cannot correspond to zero t. for otherwise an X; or Xy would be zero, contrary to the assumption that vo was at least two. The elements (ilj,),

88 (iOJ).)..., etc., are sometimes called "odd" with the remaining elements of K called "even" (28). Next consider all elements ( j) not in the basic solution and from these set G X ij < C Nj I. Since C represents the cost of adding Xij, we could not introduce this element if the resultant total cost were greater than Co. For every XAdG, it is known (49) there exists a "loop" or sequence of positions L= (,j )(i( ji ) (kj ))...) with initial and final elements (I j) and intermediate elements those of the basic solution where each such entry (c b) is preceded by (OK) or ( Ob) for some k. For any Xi 6, a sufficient condition for a solution to our original problem (that of feasibility) is that for some L; j, all degenerate basic elements are in the "odd" group. The resulting Xl is to increase all odd elements of L ij by one, with the even elements decreased by one. Then it follows that the columns Al (and hence the f; ) corresponding to the nonzero X/ have t linearly independent vectors among them. If C X, Co with X,-d 6, we then have an initial feasible solution. An extremely simple example of this two-factor problem is now given. Let the costs, Cy, be represented in the following array with the original resource constraints indicated at the side and bottom of this array:

89 Problem (1) Maximize V Subject to < Xi (V=il..,n,) V < x.j (j -,.), nD AX< b ctx< 21 X>o, Integer I 11 6 <6 3 2 14 (6 l 1 4 1 5 6 where bt= (6,6,6,2,2,2) Ct= (1,11,6,3,2,14,12,4,5). A solution to this problem (V\lo,o)is indicated: 2 0 0 0 2 0 e O _ 0 z i.e., VO =- 2, C Xo \6 xI I = 3, h -- 0 i j33 Since the value of,.(Xo) is zero in this case, the problem is:

90 Problem (2) Minimize CtX Subject to AX(d where d= (6,2,2,2,2,2,2). The transportation algorithm then yields the following solution, the exponents corresponding to the basic solution, and the main elements representing the dual variables or marginal costs. o II 5 00 0 II Cost = 16. 7 0~ 0Z Since there are two elements which are zero in this basic solution, (2, ) and (3)2), we seek an entry which has C~ <2 /-/6 =5= /o and for which a loop can be constructed as described above. It is seen that (/,3) conforms to these requirements with L/3 = ( I3)) (3) 3)) (3)2)) (2 2), (2) I)) () /)) (/)3)2. Making the adjustment indicated above, the initial feasible solution now becomes: 0 \ \ O O o Cost = 16 + 5 = 21. 0 \.

91 It is interesting to note that, with the proper interpretation of the factors, this is a simple case of an incomplete block arrangement, one of optimal Ar(X), provided the constraint 2,j 6 is added (15) (also a solution to the problem of Chapter III). Of course, full advantage may be taken of this two-factor, transportation problem analogy by utilizing existing algorithms for the solution to the latter. The additional constraint that exists is that the solution be nondegenerate. Such a constraint might be ignored initially to obtain what might be an excellent first step toward a solution to the factor arrangement problem. This would allow more complex cost functions and resource constraints. For instance, if each cell, or factor combination has an upper bound and if a convex piecewise linear cost function, rather than strictly linear, exist, the capacitated transportation algorithm of Graves and Thrall (28) could be used. It might be quite realistic that increased marginal costs occur for certain factor combinations, and such an effect should be included. B.7 Proof of the minimum cost algorithm of section B.2 The proof of the algorithm attaining a minimum cost solution in section B.2 is embraced by the following Lemma, due to E. Lawler. Private communication.

92 Lemma: Given a set of vectors {oC E()' i=l,' with associated costs C;)O, the minimum cost algorithm of section B.2 results in a minimum Z= /C,, the sum taken over all vectors chosen for which the chosen set is a basis for E (the original set spans E(N). Proof: (i) Let U.{o (,,oda.o^N be a solution using the algorithm, with associated cost CI. Let U2- {,) A2,, /} be a set of the original vectors which is also a basis for E, and has optimal cost C2 (any optimal set has precisely N vectors). Suppose C 2 C. (ii) Let &o be the minimum cost vector, not necessarily unique, in UI but not in UJ. One exists since UI~t. 2 may be represented as a unique linear combination of vectors in UJ. This nontrivial combination must include at least one vector not in U 1 since that set is a basis for E (. Let V/ be the set of vectors in this unique linear combination which are in UJI, and \Z be the set of vectors also in this linear combination in Up2, but not (Jl

93 (iii) Let K be the vector in V2 of maximal cost, C.,. If Co is the cost of o(, Co cannot be less than C,, for otherwise ( could replace s in )2 resulting in a basis at lower cost than optimal. But Cm cannot be A less than Co since by construction, X was the lowest cost vector, independent of \V. Therefore CMCo.. (iv) Replace g by o, denoting the set by WI. Repeating this process until no vector exists in U! which is not in W9, we have U1: Jk, each of which have costs CI, C2, unequal by assumption. By this contradiction it follows that all sets of optimal cost must have cost equal to C, and hence the algorithm results in one such solution.

94 APPENDIX C COMPUTATIONAL CONSIDERATIONS C.1 Problem sizes In order to assess present-day computational capability in solving the two problems of Chapters II and mI, and the set proposed in Appendix B, the problem size relations are first defined. Let a linear model have "-" factors at n, n,.., nr levels. Define pl+ En;i, and pi+r(~lZl. Then (PxP represents the dimension of the information matrix in the additive model, and (P2.XP) that size with the two-factor interaction terms included. No higher order interaction terms will be considered for purposes of computation calculation. Thus, in a 2 factorial arrangement, with K= 6 (six factors each at two levels), the information matrix (and also the Q matrix of the quadratic program) is (13x13), or written as (13,13), for the additive case and (28,28) if the two-factor interactions are included. It will be assumed that the number of inequality constraints is equal to, + I (one per factor level, one for the total observations, and one for the cost constraint). Furthermore, there are Pf ( Pa with the two-factor interactions) variables and all but the cost constraint coefficients will be (0, I). The percent of nonzeroes in

95 the coefficient constraint matrix is approximately P- n,, where N N is the total entries in the matrix, equal to,(ptI) or pl(pil). All the computation times reported refer to a 7090 processor unless otherwise stated. C.2 The eigenvalue search problem In the search technique of Chapter II, each step required (a) calculation of the minimum eigenvalue, and (b) calculation of a "step' direction. In the following, it will be assumed that the integer requirement may be ignored. Calculation of the eigenvalues of a real symmetric matrix reportedly depends very much on the technique used. Estrin and Viswanathan (22) predict several computation times for various matrix dimensions. Included below in Table 5 is the maximum time of a method by Corbato (faster than Jacobi's method) labeled column 2, and the minimum times, predicted by a special-purpose computer configuration and algorithm, column 3. Here, h refers to the number of rows in the matrix. The second phase in section 2.4 is the solution of a convex program. Actually, the objective function is linear, and all the constraints are linear with the exception of the simple form, yei < K. Being optimistic in that this program could be solved

96 TABLE 5 COMPARISON OF EIGENVALUE CALCULATIONS FOR VARIOUS MATRIX DIMENSIONS (1) (2) (3) (1)' C ~Corbato Method Special (sec.) (sec.) 10 1 0.3 25 18 2.8 50 145 12.9 100 1150 123.0 in approximately the same time as a quadratic program, Wolfe (56) reports that the latter involves about the same computing time as a linear program of equivalent dimension. On this basis, Smith (48) quotes such times for various dimensions and density of zero elements in the constraint matrix (for linear programs) as shown in Table 6. Using the most optimistic times above, Table 7 shows calculations for one step in the gradient search technique for typical problem configurations noted. Here, T1 is the eigenvalue calculation time estimate, and T2 is the step calculation time estimate,

97 TABLE 6 EXPERIENCE IN LINEAR PROGRAM CALCULATION TIMES Percent Zeroes T. Dimension in Constraint (e Matrix 5x10 1.8 27x20 50 24 118x225 95 258 TABLE 7 OPTIMISTIC COMPUTATION TIME PER ONE STEP IN GRADIENT SEARCH DimenFactor- Dimen- sion of Pct T1 T2 Factor- Pet. T1 12 Levels Description sion of Con-.) (sec.) LevelS straint 0's (sec.) (sec.) S straint Matrix 23 additive 7x7 8x7 79 1 2 23 2-factor 10x10 8x10 85 1 2 interaction 26 additive 13x13 14x13 87 1 2-24 26 2-factor 28x28 14x28 94 4 12-24 interaction 3323 additive 18x18 19x18 88 2 12-24 4x5x6 additive 16x16 17x16 72 2 12-24 4x5x6x7x8 2-factor 40x40 17x40 70 13 24-48 interaction

98 It is seen that with the best possible present-day specialized programming, computing, and assumptions about applicability of other reported results, the time per step is several seconds to one minute for the problems listed. If a 7090 were used and the Corbato method programmed, a factor of eight to ten increase could easily result. Actually, the latter effect could possibly be offset by requiring only the minimum eigenvalue and eigenvector to be calculated. The reduction in computing time for this case is not clear. A problem involving sixty steps could take an hour of computing time in the best of conditions. C.3 The quadratic integer program The three algorithms are now investigated as to their computational feasibility. Only in the case of the partial enumeration technique is there such experience, so that estimates will have to be made in the case of the algorithms specifically involving the quadratic objective function. The Kunzi-Oettli algorithm. —Problem II of section 3,2 is a mixed-integer program with possibly a noninteger valued objective function, which is pointed out by the authors (37) and Balinski (4), and convergence is not guaranteed with the mixed integer algorithm proposed by Gomory (27). The authors also point out that replacing

99 the objective function TT by N1 for large N might suffice, where now fT J would be required to be integer so that an all-integer algorithm could be used. Assuming that this could be accomplished, computational experience on integer programs varies widely in performance, depending upon the code and problem structure. Balinski (4), in attempting to consolidate the widely varying reports, mentions that Martin, using a special algorithm (41), solves scheduling problems of the order (70)160) < (rmny) ( (600,2200) regularly, approaching one hour each in computing time of the 7090 processor. In these problems, the constraint matrices had about 85 percent zeroes. The remaining nonzeroes had about 85 percent ones. In Table 7, section C.2, the values there are not inconsistent with Martin's. Witzgall's parabolic constraint algorithm. —Since the method of Witzgall (55) parallels that of Gomory (27) with the former method involving extra computation only in one row (corresponding to the single quadratic function) we might assume that this method would be approximately the same as Gomory's for similar size problems of a linear integer program character. Experience, however, is meager with Gomory's algorithm, although ten problems of the order (10,50) have been reported solved in presumably "reasonable" lengths of time (27).

100 Partial enumeration technique. —More computational data are available for this algorithm, used in section 3.6. The authors (38) report times of 56 to 716 seconds for a problem with twenty-seven (0 I ) variables and a constraint matrix involving integers having values up to 9. "Covering" problems were solved for twenty-five (0,j ) variables in a maximum of 2.5 minutes, with a problem of thirty-five variables unsolved after 10 minutes. One interesting fact reported is that computation time may vary by 10 or 20 to 1, depending on the initial ordering of the variables. Also, the extreme simplicity of the algorithm allows rapid modification of the rules to suit the problem (Rule 3, section 3.6 was a trivial addition due to the equality constraint). C.4 General remarks Probably the only remark of substance that could be made concerning computational feasibility of the methods in previous chapters is that for moderate-sized problems (4x5x6x7x8, say), considerable computational effort and skill might result in operational methods of solution. Smaller problems (25) could undoubtedly be solved with existing algorithms in reasonable time (minutes, not hours) provided the variables did not exceed a value of one or two.

101 As in most specific application areas of integer programming, it would seem that most significant improvement would come in investigating and taking advantage of, the specific problem structure. For instance, if the cost constraints were not present, the entire constraint matrix would consist of O's and l's, a fact leading possibly to advantageous algorithm modification. Toward this end, it would also appear tha a primary consideration would be simplicity of the initial algorithm which was being modified.

REFERENCES 1. Addleman, S. "Irregular Fractions of the 2n Factorial Experiments," Technometrics, Vol. 3, No. 4, November 1961. 2. Anderson, T., An Introduction to Multivariate Statistical Analysis, J. Wiley and Sons, 1958. 3. Balinski, M., "Notes on Operations Research," 1962 Engineering Summer Conference, University of Michigan. 4. Balinski, M., "Integer Programming: Methods, Uses, Computation," Management Science, November 1965. 5. Barankin, W., and R. Dorfman, "On Quadratic Programming," University of California Publication in Statistics, Vol. 2, 1958. 6. Bellman, R., Introduction to Matrix Analysis, McGraw-Hill, 1960. 7. Benders, J. F., A. R. Catchpole, and C. Kuiken, "Discrete Variables Optimization Problems," Rand Symposium on Mathematical Programming, March 16-20, 1959. 8. Boot, J. C., Quadratic Programming, North-Holland Publishing Company, 1964. 9. Bose, R. C., I, M. Chakravarti, and D. E. Knuth, "On Methods of Constructing Sets of Mutually Orthogonal Latin Squares Using a Computer," Technometrics, Vol. 2, No. 4, November 1960. 10. Box, G. E. P., and D. W. Behnken, "Some New Three Level Designs for the Study of Quantitative Variables," Technometrics, Vol. 2, No. 4, November 1960. 102

103 11. Box, G. E. P., and J. S. Hunter, "The 2k-P Fractional Factorial Designs," Technometrics, Vol. 3, No. 3, August 1961. 12. Chakrabarti, M. C., "On the C-Matrix in Design of Experiments," Journal of the Indian Statistical Association, 1959. 13. Charnes, A., and W. Cooper, Management Models and Industrial Applications of Linear Programming, Vol. 11, J. Wiley and Sons, 1961. 14. Connor, W., and S. Young, "Fractional Factorial Designs for Experiments with Factors and Two and Three Levels," National Bureau of Standards Applied Mathematics Series 58, U.S. Government Printing Office, Washington, D.C., September 1, 1961. 15. Cox, P., Planning of Experiments, J. Wiley and Sons, 1958. 16. Dantzig, G. B., "On the Significance of Solving Linear Programming Problems with Integer Variables," Econometrika (28), 1960. 17. Dantzig, G. B., "Quadratic Programming, a Variant of the WolfeMarkowitz Algorithm," Research Report 2 of the Operations Research Center of the University of California, 1961. 18. Dwyer, P., "The Direct Solution of the Transportation Problem with Reduced Matrices," Journal of Management Science, Vol. 13, No. 1, September 1966. 19. Ehrenfeld, S., "On the Efficiency of Experimental Designs," The Annals of Mathematical Statistics, Vol. 26, No. 2, June 1955. 20. Ehrenfeld, S., "Complete Class Theorems in Experimental Design," Third Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1, 1956. 21. Elston and Bush, "The Hypothesis That Can Be Tested When There Are Interactions in the Analysis of Variance Model," Biometrics, December 1964.

104 22. Estrin, G., and C. Viswanathan, "Organization of a'Fixed-Plus Variable' Structure for Computation of Eigenvalues and Eigenvectors of Real Symmetric Matrices," Journal of the Association for Computing Machinery, January 1962. 23. Falk, J., "A Constrained LaGrangian Approach to Nonlinear Programming," Technical Report 06122-8-T, University of Michigan Office of Research Administration, March 1965. 24. Fisher, F. A., Statistical Methods for Research Workers, Hafner Publishing Company (13th Edition, 1963). 25. Gatty, R., "Statistical Models for Experiments in Merchandising," 1965 Proceedings, American Statistical Association, Business and Economic Statistics Section. 26. Goertzel, G., and N. Tralli, Some Mathematical Methods of Physics, McGraw-Hill, 1960. 27. Gomory, R., "An Algorithm for the Mixed Integer Problem," RM-2597 Rand Corporation, July 7, 1960, 28. Graves, G. W., and R. M. Thrall, "Capacitated Transportation Problem with Convex Polygonal Costs," Rand Memorandum RM-4941-PR, April 1966. 29. Graybill, F., An Introduction to Linear Statistical Models, I, McGraw-Hill, 1961. 30. Graybill, F., and W. Pruitt, "The Staircase Design: Theory," Annals of Mathematical Statistics, (29) 1958. 31. Hadley, G., Linear Programming, 1962, Addison-Wesley. 32. Hadley, G., Nonlinear and Dynamic Programming, 1964, AddisonWesley. 33. Johnson, N., and F. Keone, Statistics and Experimental Design, Volume II, J. Wiley and Sons, 1964. 34. Kempthorne, O., Design and Analysis of Experiments, J. Wiley and Sons, 1952.

105 35. Kempthorne, O., "The Efficiency Factor of an Incomplete Block Design," The Annals of Mathematical Statistics, Vol. 27, 1956. 36. Kshirsagar, A., "A Note on Incomplete Block Designs," The Annals of Mathematical Statistics, Vol. 29, 1958. 37. Kunzi, H., and W. Oettli, "Integer Quadratic Programming," Recent Advances in Mathematical Programming, Edited by R. Graves and P. Wolfe, McGraw-Hill, 1963. 38. Lawler, E., and M. Bell, "A Method for Solving Discrete Optimization Problems," Operations Research, November-December 1966, Vol. 14, No. 6. 39. Lax, M., "Localized Perturbations," Physical Review, Vol. 94, 1954. 40. Mann, H., Analysis and Design of Experiments, Dover, 1949. 41. Martin, G., "An Accelerated Euclidean Algorithm for Integer Linear Programming" in (37). 42. Plackett, R., and J. Burman, "The Design of Optimum Multifactor Experiments," Biometrika (33), 1946. 43. Rao, C., "A Note on Balanced Designs," Annals of Mathematical Statistics, 1953. 44. Rao, C., Linear Statistical Inference and Its Applications, J. Wiley and Sons, 1965. 45. Rhode, C., and J. Harvey, "Unified Least Squares Analysis," Journal of the American Statistical Association, June 1965. 46. Roy, S., and Gnanadesikan, "Some Contributions to ANOVA in One or More Dimensions: I," Annals of the American Statistical Association, Vol. 30 (1959). 47. Shah, B., "On Balancing in Factorial Experiments," Annals of Mathematical Statistics, Vol. 29, No. 3, September 1958.

106 48. Smith, D., and W. Orchard-Hays, "Computational Efficiency in Product Form LP Codes," Recent Advances in Mathematical Programming, McGraw-Hill, 1963. 49. Thrall, R., "Lecture Notes on Operations Research," Engineering Summer Conferences 1966, University of Michigan. 50. Wald, A., "On the Efficient Design of Statistical Investigations," The Annals of Mathematical Statistics, Vol. 14, No. 2, June 1943. 51. Webb, S., "Construction and Comparison of Non-Orthogonal Incomplete Factorial Design," Proceedings of the Eleventh Conference on the Design of Experiments, Army Research Development and Testing ARD Report 66-2. 52. Webb, S., "Design, Testing and Estimation in Complex Experimentation," ARL 65-116, Part 1, June 1965, Aerospace Research Laboratories. 53. Wilde, D., Optimum Seeking Methods, Prentice-Hall, 1964. 54. Wilkinson, J., The Algebraic Eigenvalue Problem, Clarendon Press, 1965. 55. Witzgall, C., "An All-Integer Programming Algorithm with Parabolic Constraints," Journal Society Industrial and Applied Mathematics, Vol. 11, No. 4, December 1963. 56. Wolfe, P., "Methods of Nonlinear Programming," p. 67 in Recent Advances in Mathematical Programming, McGraw-Hill, 1963.