INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 165—188 (1997) OPTIMAL DESIGN WITH DISCRETE VARIABLES: SOME NUMERICAL EXPERIMENTS MIN-WEI HUANG* AND JASBIR S. ARORAs Optimal Design ¸aboratory, College of Engineering, ¹he ºniversity of Iowa, Iowa City, IA 52242, º.S.A. SUMMARY Continuous—discrete variable non-linear optimization problems are defined and categorized into six different types. These include a full range of problems from continuous to purely discrete and non-differentiable. Methods for solution of these problems are studied and their characteristics are catalogued. The branch and bound, simulated annealing and genetic algorithms are found to be the most general methods for solving discrete problems. After some enhancements, these and two other methods are implemented into a program for certain applications. Several example problems are solved to study performance of the methods. It is concluded that solution of the mixed variable non-linear optimization problems usually requires considerable more computational effort compared to the continuous variable optimization problems. In addition, there is no guarantee that the best solution has been obtained; however, good practical solutions are usually obtained. KEY WORDS: engineering design; optimization; discrete variables; numerical methods; test problems; evaluation of methods 1. INTRODUCTION In many engineering applications, the variables of the optimization problem cannot have arbitrary values. Instead, some or all of the variables must be selected from a list of integer or discrete values for practical reasons. For example, structural members may have to be selected from sections available in standard sizes, member thicknesses may have to be selected from the commercially available ones, the number of bolts for a connection must be an integer, the number of reinforcing bars in a concrete member must be an integer, etc. Such problems appear to be easier to solve due to a smaller solution set. In reality, they require more computations because many designs must be evaluated. The purpose of this paper is to define a general non-linear optimization problem with mixed discrete and continuous design variables, and to discuss various methods for its solution. Depending on the type of design variables and problem functions, the problem is categorized into six types. These categories are explained which cover a full range of problems; i.e., continuous, discrete, non-differentiable, and combinatorial. Based on a recent review on the subject,1,2 a few methods for solution of some of the problem types are selected for computer implementation. Certain enhancements of the methods are proposed and suitability of each method to solve different problem types is identified and cataloged. The implemented methods are evaluated using * Graduate Research Assistant s Professor, Civil and Environmental Engineering, and Mechanical Engineering CCC 0029—5981/97/010165—24 ( 1997 by John Wiley & Sons, Ltd. Received 5 February 1996 Revised 17 May 1996 166 M.-W. HUANG AND J. S. ARORA a set of structural design and other test problems. Most of the test problems have differentiable functions and the discrete variables can have non-discrete values during the solution process (i.e., Type 2 problems defined later). Conclusions drawn from the study are presented and discussed. It is important to note that many papers exist on the topic of optimal design with discrete variables. Only a few of these articles are cited in the present paper; for a more detailed review and citation of many papers, Arora et al.1 should be consulted. Section 2 defines the discrete-continuous variable non-linear optimization problem. Section 3 contains discrete variable optimization methods that have been selected for implementation into a computer program. Section 4 contains a description of the numerical examples and results for the examples. Finally, Section 5 contains some concluding remarks. 2. PROBLEM DEFINITION An optimization problem can be expressed as a general Non-linear programming (N¸P) problem of the following form: minimize f (x) subject to g (x)"0; j"1, . . . , p j g (x))0; j"p#1, . . . , m (1) j x )x )x ; i"1, . . . , n iL i iU where f and g are objective and constraint functions respectively, x and x are lower and upper j iL iU bounds for the variable x , and p, m and n are the numbers of equality constraints, total i constraints and design variables, respectively. The functions f and g are usually assumed to be j twice continuously differentiable. When some of the variables are discrete and others are continuous, we get the so-called Mixed Discrete-Continuous Non-linear Programming problem (MDNLP). The following formulation presents the mixed case of the nonlinear programming problem (1): minimize f (x) subject to g (x)"0; j"1, . . . , p j g (x))0; j"p#1, . . . , m (2) j x 3D , D "(d , d , . . . , diq ); i"1, . . . , n i i i i i1 i2 $ x )x )x ; i"n #1, . . . , n iL i iU $ where n is the number of discrete variables, D is the set of discrete values for the ith variable, q is $ i i the number of discrete values for the ith variable, d is the kth discrete value for the ith variable. ik In general, the number of available discrete values for each variable may be different. Since the variables having integer values can be considered as discrete, a distinction between these variables is not needed. Also 0—1 variable non-linear optimization (where a variable may or may not exist) is considered a special case of integer variables with only two allowable values, 0 and 1. It is important to understand some of the features of the MDNLP model (2): 1. The inclusion of equality constraints in the problem may result in an infeasible problem. This is due to the reason that for the allowable discrete values, it may not be possible to satisfy the equality constraints. In such cases, a solution closest to the feasible domain may OPTIMAL DESIGN WITH DISCRETE VARIABLES 167 have to be accepted, or the problem and the equality constraints have to be re-examined to provide a relaxed feasible domain. 2. No inequality constraint may be active at the optimum point because the constraint surface may not pass through any of the discrete points. So, in numerical calculations, only a point closest to the constraint boundary may be found. 3. There are no simple criteria to terminate the search process for optimal design. Thus, local optimality of the solution point cannot be assured unless an exhaustive search is performed. 4. The size of discreteness (i.e., how widely the allowable discrete values are separated) may govern the behaviour of some of the algorithms as well as the final solution to the problem. In many applications, optimization must be performed over a range of alternative standard designs available off-the-shelf, say, several types of engines, connections, crank shafts, electric motor types, etc. Each alternative has its own design variables and constraints which may have same physical interpretations but they act as different designs when optimization is performed. This is sometimes referred to as topological decision-making problem.3 Often the zero—one variables are introduced to represent the choice among alternative designs (zero for not used; one for used). For non-linear optimization problems, this approach amounts to essentially enumeration of all the possibilities. Zero—one variable approaches can be used along with sequential linear programming to solve MDNLP problems. These are discussed later in the paper. Depending on the type of variables and functions, the mixed continuous—discrete variable problems are classified into the following six types. Note that the convexity and linearity of the problem functions are not used in this classification since most engineering problems have nonlinear functions and are usually nonconvex (therefore only a local minimum can be expected). Problem ¹ype 1: Continuous design variables, functions twice continuously differentiable (standard NLP problem). Problem ¹ype 2: Mixed design variables, problem functions twice continuously differentiable; discrete variables can have non-discrete values during the solution process (i.e., functions can be evaluated at non-discrete points). Problem ¹ype 3: Mixed design variables, functions non-differentiable; discrete variables can have non-discrete values during the solution process. Problem ¹ype 4: Mixed design variables, problem functions may or may not be differentiable; some of the discrete variables must have only discrete value during the solution process. Problem ¹ype 5: Mixed design variables, some of the discrete variables are linked to others; assignment of a value to one variable specifies values for others; problem functions may or may not be differentiable. Problem ¹ype 6: Combinatorial problems. These are non-differentiable discrete problems. An example of such problems is the travelling salesman problem where the total distance travelled to visit a number of cities needs to be minimized. A set of integers (cities) can be arranged in different orders to specify a travel schedule (a design). A particular integer can appear only once in a sequence. A variation of this would be to select a unique set of n numbers from a much larger $ pool of contenders. Another example of this type of problem can be found in Reference 4 where a number of manufactured members need to be arranged to form a truss structure having minimum distortion. In these six type of problems, Type 1 can be formulated as problem (1), and Types 2—4 can be formulated as problem (2). Types 5 and 6 are two special classes of practical problems. Type 2 is the most commonly discussed problem in the literature and many methods have been presented to solve these types of problems. For Type 4, some of the problem functions can be evaluated only 168 M.-W. HUANG AND J. S. ARORA at the discrete values during the solution process. Examples for such variables are: number of strands in a prestressed beam or column, number of reinforcing rods, number of teeth in a gear, and the number of bolts for a joint. On the other hand, a problem is not classified as Type 4 if the effects of the non-discrete design points can be some how ‘simulated’. For instance, a coil spring must have integer number of coils. However, during the solution process, having non-integer number of coils is acceptable as long as function evaluations are possible (they may or may not have physical meaning). 3. MIXED VARIABLE OPTIMIZATION METHODS In recent years, there has been considerable interest in the area of discrete variable optimization due to the obvious reasons that most practical applications have such variables. Several review articles on the methods for optimization of nonlinear problems having mixed variables have been recently published.1,3,5~9 Among the methods for discrete variable non-linear optimization problems, the following techniques have been most commonly discussed: branch and bound methods (BBM), zero—one variable techniques (integer programming), sequential linear programming, rounding-off techniques, cutting plane techniques, heuristic techniques, penalty function approach, Lagrangian relaxation (duality) approaches, and sequential linear programming (SLP). More recently, simulated annealing and genetic algorithms have been discussed. Each of the methods has some positive and some negative features. Therefore, there are very few methods that can be applied effectively to all types of problems. In Sections 3.1—3.5, some of the more effective methods, along with their enhancements, are described. 3.1. Branch-and-Bound Method (BBM) Branch and bound is perhaps the most widely known and used method for discrete optimization problems. The method was originally developed for linear discrete problems; however, it has been used for mixed variable nonlinear problems as well. When applied to linear problems, the method can be implemented in a way to yield a global minimum point; however, for non-linear problems there is no such guarantee, unless the problem is convex. Here two forms of the method are explained in detail, one for linear discrete problems and the other for non-linear mixed variable problems where discrete variables can be treated as non-discrete during the solution process. More discussion about the method can be found in many other references, such as References 1, 10 and many references cited in there. The branch and bound is basically an enumeration method in which one systematically tries to reduce the number of trials to reach the minimum point (that is why it is sometimes called implicit enumeration method). The concepts of branching, bounding and fathoming are used to achieve these objectives. To explain these concepts, consider a linear problem of Type 4 as follows: f "!20x !10x 1 2 subject to g "!20x !10x #75)0 1 1 2 g "12x #7x !55)0 2 1 2 g "25x #10x !90)0 3 1 2 x 3M0, 1, 2, 3N, x 3M0, 1, 2, 3, 4, 5, 6N 1 2 x and x cannot have non-integer values during the solution process. 1 2 minimize OPTIMAL DESIGN WITH DISCRETE VARIABLES 169 Figure 1. A branch and bound method without solving continuous subproblems Since x and x can have 4 and 7 possible values respectively, the full enumeration will require 1 2 evaluation of 28 combinations; however, BBM can find the final solution in much fewer evaluations. For the problem, the derivatives of f with respect to x and x are always negative. 1 2 Thus, one can enumerate the discrete points in the descending order of x and x to ensure that 1 2 the cost function is always increased when one of the variables is perturbed to the next lower discrete value. The BBM is illustrated in Figure 1. For each point (called node), the cost and constraint function values are shown. From each node, two more nodes are branched out where functions are evaluated. If there is any constraint violation at a node, further branching is necessary. Once a feasible solution is obtained, the node need not be branched further since no point with lower cost is possible from there. Such nodes are said to have fathomed; i.e., reached their lowest point on the branch and no further branching is necessary from them. Nodes 6 and 7 are fathomed this way where the cost function has a value of !80. This value becomes an upper bound for the cost function, and any node having cost function value higher than the current bound is fathomed. Nodes 9, 10 and 11 are fathomed because the design is infeasible with the cost function value larger than or equal to !80. Thus the global solution for the problem is found in 11 evaluations. There are many variations of the BBM for linear discrete optimization problems. For example, the problem can be transformed to a linear integer programming problem by introducing 0—1 variables. It can be then solved by the BBM described in the foregoing paragraph. Many other integer programming methods are available to solve the problem.1 Several implementations of BBM are described in a book by Moré and Wright.11 For general MDNLP problems, the BBM has also been used in several different ways. For example, one can use the branching, bounding and fathoming procedures as described previously 170 M.-W. HUANG AND J. S. ARORA for the linear problems. Another way would be to linearize the problem at each iteration, if the functions are differentiable. The linearized problem can then be solved using either the BBM as described previously, or the 0—1 variables approaches. For non-linear problems where the discrete variables can have non-discrete values during the optimization process and the problem functions are differentiable (Problem Type 2), a variation of the BBM that is implemented in the present work is described below. In this procedure, one takes advantage of a local minimization procedure to reduce the number of nodes. Initially, an optimum point is obtained by treating all the variables to be continuous. If this solution is discrete, then the process is terminated. If one of the desired variables is not discrete, then its value lies between two discrete values; e.g., d (x (d . Now two subproblems are defined, one ij i ij`1 with the constraint x )d and the other with x *d . This process is called branching. The i ij i ij`1 two subproblems are solved again using the local minimization procedure, and the optimum solutions are stored as nodes of the tree containing optimum values of the variables, the cost function and the appropriate bounds on the variables. This process of branching and solving continuous problems is continued until a feasible solution to the MDNLP is obtained. The cost function corresponding to this solution becomes an upper bound on the optimum solution for the original MDNLP. From this point onwards, all the nodes of the tree that have cost function value (obtained by solving a continuous variable optimization subproblem) higher than the established upper bound are eliminated from further consideration; i.e., they are fathomed. The process of branching and fathoming is repeated from each unfathomed node. A better upper bound for the optimum cost function is established when a feasible discrete solution is obtained with the cost function smaller than the current upper bound. The nodes may be fathomed in any of the following three ways: (i) feasible discrete solution to the continuous solution with the cost function value higher than the current upper bound, (ii) infeasible continuous problem, and (iii) the optimal value of the cost function for the continuous problem higher than the upper bound. The search for the optimum terminates when all the nodes have fathomed. It is important to note that the BBM is guaranteed to find the global optimum only if the problem is linear or convex (assuming that fathoming is done in a valid way). In the case of general non-linear nonconvex problems, it is possible that a node is fathomed too early as one of its branches may actually contains the true global solution. The BBM has been used successfully to deal with problems with discrete design variables and has been proved to be quite robust. However, for problems with a large number of discrete design variables, the number of subproblems (nodes) becomes large, making the method inefficient. This drawback is more severe for non-linear optimization problems. Thus, improving performance of the method has been a major interest of many researchers.12,13 The Algorithm BBM (given later) is implemented, after several enhancements, for Type 2 problems only. These enhancements enable the method to create fewer subproblems (nodes) without neglecting those with possible lower cost function value. Others improve the efficiency and reliability when solving the subproblems. These enhancements are described as follows: 1. »ariables to be branched. After a subproblem is solved, two new subproblems need to be created if the solution is not acceptable (acceptable means that each of the discrete variables has a value from its allowable discrete set). The first step of creating new subproblems is to determine the variable for branching. For this, the following criteria are used: (i) If a variable has been previously branched, do not select it again (refer to ‘Re-branching the branched variables’ described later). (ii) If a variable is at its lower or upper bound of the current subproblem, do not branch it. A variable is considered ‘on the bound’ if the distance between the bound and the value of this variable is less than a small pre-set value (e.g., 10~6). The OPTIMAL DESIGN WITH DISCRETE VARIABLES 171 small value is used to prevent numerical difficulties. For example, the subproblem solver may find design variables with values ‘close to’ but never ‘exactly on’ the bounds depending on the convergence criteria used. If such close-to-bound variables are considered for further branching, an excessive number of subproblems will be created, and the same subproblems will be created repeatedly. (iii) Otherwise, choose a variable to branch according to the next enhancement. 2. Order of branching of the variables. An upper bound for the cost function is established once a feasible discrete solution is obtained. This upper bound is then used to eliminate other branches which have a larger cost function value or which have the same value but the solutions are not discrete. One can eliminate more branches if a smaller upper bound is generated as early as possible. This can be accomplished by choosing the right variable for branching at each step. Many algorithms with different branching orders for variables have been proposed by researchers. Hager and Balling12 and Tseng et al.13 have proposed a max-cost difference method for choosing the variable for branching. This method evaluates the cost function at new nodes directly by assigning each variable the lower or upper bound value. To explain the criterion, calculate the following quantities: Cost l"f (x , . . . , d , . . . , x ); Cost "f (x , . . . , d , . . . , x ); *Cost "DCost l!Cost D i 1 ij n i6 1 ij`1 n i i i6 where Cost l and Cost are, respectively, the estimated cost function values of the first and i i6 second nodes for the ith design variable, and *Cost is the estimated cost difference for the i ith design variable when it is assigned to its lower and upper discrete values. The variable chosen to be branched corresponds to max (*Cost , *Cost , . . . , *Cost ). 1 2 n 3. Re-branching the branched variables. After a variable is branched and two subproblems are created and solved, the two solutions will have discrete values for the variable just branched. However, a previously branched variable may not remain discrete. Although one can again put such variables among the choices for further branching, this approach has two drawbacks: (i) a branched non-discrete variable has a tendency to become non-discrete again even it has been branched several times (this results in creation of excessive subproblems), and (ii) a branched non-discrete variable is likely to be chosen again, so the other neverbranched variables stand less chance of being branched. This may delay the search for a feasible discrete solution which is needed as an upper bound to eliminate many branches. To overcome this difficulty: (i) we always first select the never-branched variables for branching (first level branching), and (ii) if all discrete variables have been branched, we then branch those variables that were branched before but do not have discrete values currently (second level branching). This procedure is allowed to be repeated for a pre-set number of times (see ILEVEL in Algorithm BBM) if some variables keep regaining the non-discrete values. Note that some researchers fix the values of the discrete variables once they have been branched. This may reduce the amount of subproblems created. However, it may also eliminate some promising branches prematurely since subproblems may become infeasible if some variables are fixed. 4. Order of branching of the subproblems. Two basic strategies, ‘depth first search’ and ‘breadth first search’ were discussed by Ringertz,14 each having its advantages and drawbacks. Our experience shows that the ‘best first search’ suggested by Tseng et al.13 often reaches the solution by creating and solving fewer subproblems. In this approach, the subproblem selected to be solved next was created from the solution of the subproblem that had the smallest cost function value. 172 M.-W. HUANG AND J. S. ARORA 5. Synchronizing the design bounds. For continuous subproblems, it is important that the imposed lower and upper bounds on each variable correspond to the smallest and largest values in the discrete set. Otherwise, unnecessary subproblems may be created. This is explained by the following two cases. (i) The lower bound is smaller than the smallest discrete value and the continuous solution of the parent problem falls in between: Only one subproblem can be created; however, there would be no need to create this subproblem if the lower bound had been set to the smallest discrete value. Furthermore, if another discrete variable was selected for branching before that one (which is in an infeasible area), the subsequently created subproblems would be affected (many of them would be infeasible or redundant). (ii) The lower bound is larger than the smallest discrete value: For many design problems, several variables reach their lower bounds at the continuous solution. If the lower bound had been set to the corresponding lowest discrete value, many variables would already have discrete values at continuous solution, and so fewer subproblems would be created. 6. Starting points for solving the subproblems. It is recommended that the solution of a solved parent problem (with the branched variable perturbed to the new bound) be used as the starting points of the two new created subproblems. This usually reduces the number of iterations needed for solving a subproblem. 7. Stopping criteria for solving the subproblems. When solving the subproblems at the initial stages, very precise solutions are usually not needed, therefore slight constraint violation is acceptable. This can reduce the number of iterations when solving many subproblems. Also many subproblems may be infeasible due to the newly imposed design bounds and so the subproblem solver may keep trying to find a solution until the allowable number of iterations is reached. Thus, a moderate value for the allowable number of iterations is recommended. Algorithm BBM: Branch-and-Bound Method Step 1: Assuming all design variables to be continuous, solve problem (1). Set IROUND"1. Step 2: If the solution is not feasible, stop. Otherwise, continue. Step 3: If the solution is discrete, save it and stop. Otherwise, continue. Step 4: Branch from the above solution and create two subproblem (nodes). The variable to be branched is chosen according to the following criteria: (i) If a variable has been branched, do not branch it. (ii) If a variable is on the bound, do not branch it. (iii) Choose a variable for branching using the max-cost difference method (i.e., Enhancement d2). If no variable can be chosen, then set IROUND"IROUND#1. If IROUND'ILEVEL, stop. Otherwise, treat all the non-discrete variables as never-branched ones and repeat this step. ILEVEL is the allowable number of branching levels (Enhancement d3). Step 5: Save the new nodes into the node-list which contains all of the unsolved subproblems. Step 6: If the node-list is empty, stop. Otherwise, select a node from the node-list according to the best first search criterion (i.e., Enhancement d4). Step 7: Solve the node. Step 8: If the solution is not feasible, delete the solved node from the node-list and go to Step 6. Otherwise, continue. Step 9: If the solution is discrete, then check the following two criteria: (i) If upper bound f * does not exist, then set f *"f and save this solution as the best discrete solution. (ii) If f * exists and f*f *, then delete the solved node from the node-list and go to Step 6. (iii) If f * exists and f(f *, then set f *"f and replace the best discrete solution by this one. Delete the solved node from the node-list and go to Step 6. OPTIMAL DESIGN WITH DISCRETE VARIABLES 173 If the solution is not discrete, then check the following two criteria: (i) If f * exists and f*f *, then delete the solved node from the node-list and go to Step 6. (ii) If f * does not exists or f(f *, then delete the solved node from the node-list and go to Step 4. The foregoing branch-and-bound algorithm solves each node as a continuous problem, and then creates two nodes based on the continuous solution. In our implementation, since an SQP method is used as the continuous problem solver, differentiable functions are required. 3.2. Simulated Annealing (SA) Simulated annealing is a stochastic technique to find a global minimizer for continuous—discrete—integer non-linear programming problems.15,16 The basic idea of the method is to generate a random point and evaluate the problem functions. If the trial point is infeasible, it is rejected and a new trial point is generated. If the trial point is feasible and the cost function value is smaller than the current best record, then the point is accepted, and record for the best value is updated. If the point is feasible but the cost function is higher than the best value, then the point is sometimes accepted and sometimes rejected. The acceptance is based on value of the probability density function of Bolzman—Gibbs distribution. If this function has value greater than a random number, then the trial point is accepted as the best solution even if its cost function value is higher than the recorded best value. In computing the probability a parameter called the temperature is used. For the optimization problem, this temperature can be a target value (estimated) for the cost function corresponding to a global minimizer. Initially, a larger target value is selected. As the trials progress, the target value is reduced (this is called the cooling schedule), and the process is terminated after a fairly large number of trials. The acceptance probability steadily decreases to zero as the temperature is reduced. Thus in the initial stages, the method is likely to accept worse designs while in the final stages, the worse designs are almost always rejected. This strategy avoids getting trapped at a local minimizer. The main deficiencies of the method are the unknown rate at which the target level is to be reduced, and uncertainty in the total number of trials and in the number of trials after which the target level needs to be reduced. Balling17 and May and Balling18 have proposed a filter for the method to improve its performance for structural design problems. As defined in Section 2, there are five types of mixed variable optimization problems (excluding Type 1), each having different types of design variables and/or functions. The implemented BBM can only handle Type 2 problems because its subproblem solver requires differentiable functions. On the other hand, the performance of the SA, or the genetic algorithm (presented in the next section), is not greatly affected by the variables or function type. Therefore, the method is applicable to all types of problems. Some implementation details and improvement of the method are described as follows: 1. Step size. The SA method generates a series of design points during its optimization process. Only one point is generated at a time, and the next point is generated within a certain neighbourhood of the current point. Thus, although SA randomly generates design points without the need for function or gradient information, it is not a pure random search within the entire design space. At the early stage, a new point can be located far away from the current point to speed up the search process and to avoid getting trapped at a local minimum point. After the temperature gets low, the new point is usually created nearby in order to focus on the local area. Thus, the step size is reduced as the temperature becomes smaller. In the implemented method, the next design point is calculated as follows: Randomly, select a variable, e.g., the ith variable x to be perturbed. The next design point will have a new value for the ith variable calculated as i follows (other variables will remain unchanged): 174 M.-W. HUANG AND J. S. ARORA Let z be a random number uniformly distributed in 0 and 1. (i) Continuous variables: x(k`1)"min (x , x(k)#a (x !x )), if z)0·5 iU i iU iL i x(k`1)"max (x , x(k)#a (x !x )), if z'0·5 iL i iU iL i where x(k`1) and x(k) denote, respectively, the new and the current values for the ith i i variable, k is the trial number within the same iteration (same temperature level), x and iL x are the lower and upper bound for the ith variable respectively, and a is the step size iU calculated as follows: a"max (0·01, 0·2(0·9)K~1), where K is the iteration number. (ii) Discrete variables: If the current design point has the mth discrete value for the ith variable (i.e., x(k)"d ), im i then the new value for the ith variable becomes x(k`1)"d , where j"min (q , m#J), if ij i i z)0.5 or j"max (1, m!J), if z'0·5. The integer J is calculated as J"max (1, INT(0·2(0·9)K~1q )), where INT(x) denotes the integer part of x, and q is the number of i i discrete values for the ith variable. 2. Constraint treatment. In the implemented SA method, the newly generated point has to be feasible. If it is not, another point is generated until feasibility is attained. Another method for treating constraints is to use the penalty function approach; i.e., the constrained problem is converted to an unconstrained one. However, the choice of a good penalty parameter varies from one problem to another, making this approach difficult to implement in a general purpose optimization program. 3. Reducing the temperature. For general optimization problems, temperature is an arbitrary parameter with the same units as the cost function. It depends upon the mean values of the cost function. This parameter is sometimes called the effective temperature. The temperature should be reduced gradually at each iteration. The following temperature reducing scheme is used in the implemented method: ¹ "max (10 000, TCOST), if K"0; and ¹ "0·9 ¹ , 0 K K~1 if K*1, where ¹ and ¹ are, respectively, the initial temperature and the temperature at 0 K the Kth iteration (an iteration here implies a design search process in which the temperature remains unchanged), K is the iteration number, and TCOST is the approximated value of the cost function which is estimated by the minimum cost of 10 randomly generated feasible points. 4. Stopping criteria. The following stopping criteria are used in the implemented SA method: (i) The program stops if the change in cost function value is less than 0·0001 for the last 4 consecutive iterations. (ii) The program stops if IS/ILIM (0·05, where ILIM is the number of trials (or number of feasible points generated) within one iteration, and IS is the number of trials which satisfy * f(0 (see Step 3 of Algorithm SA). ILIM " 100 is used in the current implementation. (iii) The program stops if K reaches a pre-set value. The following algorithm is based on the above improvements to the implementation by Kincaid and Padula.4 Algorithm SA: Simulated Annealing Step 1: Choose initial temperature T (expected global minimum for the cost function) and 0 a feasible trial point x(0). Compute f (x(0)). Select ILIM (e.g., 100). Initialize the iteration counters as K"0 and k"1. OPTIMAL DESIGN WITH DISCRETE VARIABLES 175 Step 2: Generate a new point x(k) randomly, using the step size criterion described earlier. If the point is infeasible, generate another random point until feasibility is satisfied. Generate a random number z uniformly distributed in [0, 1]. Compute f (x(k)) and *f"f (x(k))!f (x(0)). Step 3: If *f(0, then take x(k) as the new best point x(0), set f (x(0))"f (x(k)) and go to Step 4. Otherwise, calculate the probability density function: p (*f )"exp (!*f/¹ ). If z(p (* f ), then K take x(k) as the new best point x(0) and go to Step 4. Otherwise go to Step 2. Step 4: If k(ILIM, then set k"k#1 and go to Step 2. If k'ILIM and any of the stopping criteria is satisfied, then stop. Otherwise, go to Step 5. Step 5: Set K"K#1, k"1; set ¹ "0·9¹ ; go to Step 2. K K~1 3.3. Genetic Algorithm (GA) Genetic algorithms belong to the category of stochastic search methods. The basis of most of the contemporary developments of the GAs is found in Holland.19 In a GA, several design alternatives, called a population in a generation, are allowed to reproduce and cross among themselves, with bias allocated to the most fit members of the population.20 Combination of the most desirable characteristics of mating members of the population results in new designs that are more fit than their parents. Thus, if a measure which indicates the fitness of a generation is also the desired goal of a design process, successive generations produce better values of the cost function. An advantage of the approach is that gradients are not needed. In a GA, each design must be represented by a finite length string. Usually, binary strings are used for this purpose. Three operators are needed to implement the algorithm: (i) reproduction; (ii) crossover; and (iii) mutation. Reproduction is an operator where an old string is copied into the new population according to the string’s fitness which is defined according to the cost function value. More fit strings (those with smaller cost function values) receive higher numbers of offspring. There are many strategies to implement this reproduction operator. The next operator—crossover—corresponds to allowing selected members of the population to exchange characteristics of the design among themselves. Crossover entails selecting a starting and ending position on a pair of mating strings at random, and simply exchanging the string of 0’s and 1’s (for a binary string) between these positions on one string with that from the mating string. Mutation is the third step that safeguards the process from a complete premature loss of valuable genetic material during reproduction and crossover. In terms of a binary string, this step corresponds to selecting a few members of the population, determining at random a location on the strings, and switching the 0 to 1 or vice versa. The above three steps are repeated for successive generations of the population until certain stopping criteria are satisfied. The member in the final generation with the best fitness level is the optimum design. Some details and enhancements of the implemented GA, given as Algorithm GA, are described as follows: 1. Defining a schema. A schema is first defined to represent a design point. It is an ¸-digit genetic string with 0 or 1 for each digit, where ¸ is the total number of binary digits required to specify a design point. A genetic string is constructed by linking n binary strings together, where n is the number of design variables. A binary string of 0’s and 1’s represents a design variable. An m-digit binary string has 2m possible 0—1 combinations implying that 2m discrete values can be represented. The following method can be used to convert a 0—1 combination to its corresponding discrete value of a variable with 176 M.-W. HUANG AND J. S. ARORA n discrete values: Calculate m which is the smallest integer satisfying 2m'n . Calculate $ $ the integer j: A B n m $ j" + ICHAR(i) 2(i~1)#1; If j'n , j"INT ( j!n ) $ $ 2m!n $ i/1 (3) where ICHAR(i) is the value of the ith digit (either 0 or 1), and INT (x) is the integer part of x. Thus the jth discrete value corresponds to this 0—1 combination. Note that when j'n , $ procedures different from the one given in equation (3) can be used to adjust j such that j)n . $ 2. Random design generation. After the schema is defined, the first population with N members (genetic strings) has to be created. These members are generated randomly via 1 the use of a random number generator. Note that there is no need to use one random number to make decision (between 0 and 1) for each of the ¸ digits in a genetic string. Instead, the following procedure shows a way to produce a 30 digit genetic string: (i) Generate two double precision random numbers. (ii) Assuming these two numbers are 0·6832212546532534 and 0·1257213658623255, a string ‘68322125465325341257213658623255’ is then created. (iii) The first 30 digits of the above string are converted to ‘110000010110010000110001111000’ in which ‘0’ is converted from any value in 0 and 4 and ‘1’ from any value in 5 and 9. 3. Reproduction procedure. Reproduction is to select the more fit members of the population into the mating pool from which members for crossover and mutation transformations are selected. The ith member in the original population has a probability of selection P "F / + F ( j"1, . . . , N ), where F is the fitness of the ith design. Figure 2 shows how i i j 1 i a member is selected from the original population according to the probability of selection P . After the wheel in the figure is spun according to the value of a random number between i 0 and 1, the arrow at the starting position points to the member that is selected for reproduction. This member is copied to the mating pool (a new population of same size, N ), 1 but it still remains in the original population for further selection. Thus, the new population may contain identical members and may not contain some members that were in the Figure 2. Reproduction process in a genetic algorithm—selection of the mating pool OPTIMAL DESIGN WITH DISCRETE VARIABLES 177 previous population. The average fitness will be increased, however. The fitness function is defined as follows: F "(1#e) f !f (4) i .!9 i where f is the cost function for the ith design, f is the largest recorded cost function value, i .!9 and e is a small value (e.g., 2]10~7) to prevent numerical difficulties due to F "0. i 4. Amount of crossover and mutation. For each generation (iteration), three operators, reproduction, crossover and mutation, are performed. While the number of the reproduction operations is always equal to the size of the population, the amount of crossover and mutation can be adjusted to fine-tune the performance of the algorithm. In the implemented GA, the following procedure is used at each generation for crossover and mutation. First determine I and I which control the amount of crossover and mutation using the following scheme: .!9 . (i) The program starts with normal rates of crossover and mutation: I "1, and I "INT (P N ) (5) .!9 . . 1 where P , a fraction of the population selected for mutation, is set to 0·3 and N is the . 1 size of the population. It is concluded based on our experience with the test problems presented later that too much crossover results in poorer designs. Apparently, it produces new designs that are far away from the mating designs. Therefore, crossover is carried out only once between two members (I "1). The mutation, however, changes .!9 designs in a neighbourhood of the current design; therefore, a larger amount of mutation is allowed. Note also that the population size for the present applications is selected as twice the number of design variables with minimum and maximum limits as 100 and 300, respectively. (ii) Let f * denote the best cost function value among the entire population at the Kth K iteration. If the improvement on f * is less than e@ (e@"0·0001 is used) for the last two K consecutive iterations, then I is doubled temporarily. This ‘doubling’ strategy con.!9 tinues at the subsequent iterations and returns to 1 as soon as f * is reduced. The concept K behind this is that we do not want too much crossover or mutation to ruin the good genes (in genetic strings) if they keep producing better offsprings. On the other hand, we need more crossover and mutation to trigger changes when progress stops. (iii) If the improvement on f * is less than e@ for the last I (I "30 is used) consecutive K ' ' iterations, then P "0·6 is used. The crossover and mutation are performed as follows: . FOR i"1, I .!9 Generate a random number z uniformly distributed in [0, 1] If z'0·5, perform crossover as follows: Randomly select two mating parents from the mating pool, and randomly choose two sites on the genetic strings and swap strings of 0’s and 1’s between the two chosen sites; e.g., ‘010110’ and ‘110001’ become ‘010000’ and ‘110111’. If z)0·5, skip crossover. FOR j"1, I . Generate a random number z uniformly distributed in [0, 1] If z'0·5, perform mutation as follows: Randomly, choose a member from the mating pool and switch a 0 to 1 or vice versa at a randomly selected site; e.g., ‘110111’ becomes ‘110011’. If z)0·5, skip to next j. ENDFOR ENDFOR 178 M.-W. HUANG AND J. S. ARORA 5. ¸eader of the population. At each generation, the member having the lowest cost function value among the entire population is defined as the ‘leader’ of the population. If more than one member have the same lowest cost, only one of them is chosen as the leader. The leader is replaced if another member with lower cost appears. It is safe-guarded from extinction (due to reproduction, crossover or mutation). In addition, the leader is guaranteed at least a 10 per cent probability of selection for reproduction. One benefit of using a leader is that the best cost function value of the population can never increase from one iteration to another, and some of the best ‘genes’ can always survive. 6. Reducing the number of analyses. Among the three operators, the crossover and mutation change the designs (genetic strings), including their cost function values and feasibility. Reproduction, however, does not. It merely copies some of the designs into the new population. Thus, it is recommended that the cost function value and the maximum constraint value be saved for each member of the population. If, after performing all three operators, a member remains intact because it was not chosen for crossover and mutation, then the saved information can be used to avoid re-calculation of the constraints and the cost function when preparing the fitness values at the next iteration. Also when checking feasibility, a trial design is rejected as soon as any violation is encountered. 7. ºse of penalty function approach. As in the SA method, the penalty function approach for treating constraints is not recommended due to the reasons stated earlier. 8. Initial strings. In many cases, the designer already has some base designs or previous designs in hand. These can be used as part of the first population. The GA is then used to generate better designs using them. 9. Stopping criteria. If the improvement on f * is less than e@ for the last I consecutive iterations K ' after I "0·6 is used, then the program stops. e@"0·0001 and I "30 are used. . ' Algorithm GA: Genetic Algorithm Step 1: Define a schema, as previously described, to represent a design point. Step 2: Randomly, generate N genetic strings (members of the population) according to the 1 schema, where N is the population size. For constrained problems, only the feasible strings are 1 accepted. Set K"0. Step 3: Define a fitness function for the problem, as in equation (5). Step 4: Calculate the fitness values for all the genetic strings. Set K"K#1, and I "1. # Step 5: Reproduction: Select the more fit members of the population into the mating pool, from which members for crossover and mutation transformations are selected. Step 6: Crossover: Select two mating parents from the mating pool. Randomly, choose two sites on the genetic strings and swap strings of 0’s and 1’s between the two chosen sites. Set I "I #1. # # Step 7: Mutation: Choose a fraction (P ) of the members from the mating pool and switch . a 0 to 1 or vice versa at a randomly selected site on each chosen strings. If, for the past I consecutive generations, the member with the lowest cost remains the same, the mutation ' fraction P is doubled. I is an integer (e.g., 30) defined by the user. . ' Step 8: If the member with the lowest cost remains the same for the past 2 consecutive generations, then I "2. Otherwise, I "1. If I (I , go to Step 6. Otherwise, continue. .!9 .!9 # .!9 Step 9: If the following stopping criterion is satisfied, stop. Otherwise, go to Step 4. Stopping criterion: After the mutation fraction P was doubled, the best value of the fitness is not . updated for the past I consecutive generations. ' OPTIMAL DESIGN WITH DISCRETE VARIABLES 179 3.4. Dynamic rounding-up A simple approach for discrete variable optimization of Type 2 problems is to first obtain an optimum solution using a continuous approach. Then using heuristics, the variables are rounded-up to the nearest available discrete values to obtain a discrete solution. This is a simple idea, but it often results in an infeasible design for problems having a large number of variables. However, it is not necessary to round-up all variables to their nearest discrete neighbours. Some of them could be decreased while others could be increased. The main concern of rounding-off approach is the selection of variables to be increased and the variables to be decreased. The strategy may not converge, especially in case of high non-linearity and widely separated allowable discrete values. In that case, the discrete minimizer need not be in a neighbourhood of the continuous solution. The following dynamic rounding-up algorithm is based on the one described in Arora.21 It is implemented and compared with the other methods. Algorithm DRU: Dynamic Rounding-Up Step 1: Set n"number of design variable. Step 2: Assume all the design variables to be continuous and solve the NLP problem. An SQP algorithm is used to obtain the solution which gives values for design variables and Lagrange multipliers. Step 3: If the solution is discrete, stop. Otherwise, continue. Step 4: FOR k"1, n Calculate the Lagrangian function value for each k with the kth variable perturbed to its upper discrete neighbor. END FOR Step 5: Choose a design variable that minimizes the Lagrangian in Step 4 and remove that variable from the design variable set. This variable is assigned the selected discrete value. Set n"n!1 and if n"1, stop, otherwise, go to Step 2 The number of additional continuous problems that need to be solved by the above method is (n!1). However, the number of design variables is reduced by one for each subsequent continuous problem. A similar dynamic rounding-off method has been recently demonstrated by Chan et al.22 for discrete optimization of tall building steel frameworks. 3.5. Neighbourhood search method When the number of discrete variables is small and each discrete variable has only a few choices, the simplest way to find the solution of a mixed variable problem may be just to explicitly enumerate all the possibilities. With all the discrete variables fixed at their chosen values, the problem is then optimized for the continuous variables. This approach has some advantages over BBM: it can be implemented easily in the existing NLP solver, the problem to be solved is smaller and the gradient information with respect to the discrete variables is not needed. However, the approach is far less efficient than an implicit enumeration method such as BBM as the number of discrete variables and size of the discrete set of values become large. To reduce the number of enumerated cases, a neighbourhood search method is proposed which first obtains a continuous solution with all the discrete variables considered as continuous. Then a certain range of discrete values near the continuous solution is selected for enumeration. 180 M.-W. HUANG AND J. S. ARORA Table I. Methods appropriate for various problems types Branch and bound Simulated annealing Genetic algorithm Dynamic rounding-up Neighbourhood search Problem type solved Can find feasible discrete solution? Can find global minimum for convex problems*? Need gradient information? 2, (3, 4, 5, 6)s 2, 3, 4, 5, 6 2, 3, 4, 5, 6 2 2 Yes Yes Yes No guarantee Yes Yes Yest Yest No guarantee Yes° Yes/Nos No No Yes/Nos Yes/Nos * Continuous problem by treating all variables as continuous s Dependent on the type of continuous variable optimization solver used t If enough execution time and proper stopping criteria are provided ° If the number of examined neighbouring points (NP) is large enough 3.6. Dimension reduction strategy Some mixed variable optimization methods discussed in the previous sections use a local optimization method such as the SQP algorithm to solve the subproblems they create. This section presents a strategy which can further improve performance of the mixed variable optimization algorithms. In this strategy, the dimension of the QP subproblem at each iteration is reduced by eliminating the variables that are either on their bounds or are expected to be on the bounds. After solving the reduced QP subproblem, the Lagrange multiplier for each bound constraint is recovered using the optimality conditions with respect to the fixed variables. If the Lagrange multipliers are non-negative, then the solution of the reduced QP subproblem is accepted; otherwise the QP subproblem is re-solved after releasing the fixed variables. This strategy is more beneficial for the branch and bound and dynamic rounding-up methods since many design variables stay on their bounds due to the newly imposed lower and upper limits on the design variables. 3.7. Discussion of methods Features of various methods and their suitability for various types of non-linear problems are summarized in Table I. Note that Problem Type 1 is not included in Table I since it is the standard NLP problem. Several methods are available for this class of problems. The characteristics documented in the table are: type of problem solved, ability to find feasible discrete solutions (assuming that there are feasible solutions), ability to find a global minimum, and the differentiability requirement for the problem functions. Note that the BBM is capable of handling problems other than Type 2 (as explained previously). The implemented BBM, however, uses the SQP algorithm as the subproblem solver which requires differentiable functions. Therefore, the current implementation can solve Type 2 problems only. Note also that some of the methods for discrete variable optimization (not described here) use the structure of the problem to speed up the search for the discrete solution. Such methods are not suitable for implementation into a general purpose software. The BBM, SA and GA are the most general methods; however, they are also the most time-consuming ones. Traditionally, SA and OPTIMAL DESIGN WITH DISCRETE VARIABLES 181 GA have been used for unconstrained problems; however, in the present work, they are applied to constrained problems. 4. NUMERICAL EXAMPLES AND RESULTS 4.1. Introduction All algorithms described in Sections 3.1—3.5 are implemented. They all handle input and output in a similar way, so the I/O does not affect the comparison of execution times. Several example problems are solved and the results are summarized to assess performance of the methods. The methods are evaluated according to the execution time and the cost function values. An 80 MHz HP-UX 715 workstation is used for all calculations. All algorithms and test problems are coded in FORTRAN77. For the branch and bound, dynamic rounding-up and neighbourhood search methods, the program IDESIGN2 is used to solve the subproblems. It is noted that most of the test problems used are of Type 2 to compare the performance of the implemented methods. However, one test problem of Type 3 and one of Type 5 are also included since they can be solved by two of the implemented methods. Discussion and implementations of methods for other problem types can be found in recent papers by the authors: Huang and Arora23 for Type 1; Huang and Arora24 and Arora and Huang10 for Type 5 (design of steel frameworks using available sections); and Huang et al.25 for Type 6. Examples and solutions for problems of Type 4 can be found in Kocer and Arora.26 4.2. Test problems Fifteen example problems are solved using the selected algorithms. Detailed descriptions of the test problems can be found in several references,27 and are not presented here. Among these, Problem 2 is of Type 5, Problem 3 is of Type 3, and others are of Type 2. Results for all problems are summarized in the next section. Problem 1: Spring design problem-case 121 The problem is to design a coil spring loaded in tension or compression. The objective is to design a minimum cost (mass) spring to carry the given load without material failure while satisfying some performance requirements. For this Type 2 problem, the design variables are: wire diameter d, mean coil diameter D and number of active coils N. Constraints are imposed on deflection, shear stress, surge wave frequency, outside diameter, and lower and upper bounds on the design variables. The design variables must also be selected from the following sets: d3M0·05, 0·06, . . . , 0·10N; D3M0·3, 0·4, . . . , 0·10N; N3M2, 3, . . . , 15N Problem 2: Spring design problem-case 2 This is the same problem as the first one. However, one of the three available materials must be used for the final design. Each material has its own shear modulus (G), mass density (o), allowable shear stress (q ) and unit price (º ) as given in Table II. This is a Type 5 problem. ! 1 Problem 3: Non-differentiable function This Type 3 problem has only one design variable. The cost function is non-differentiable with respect to x (as shown in Figure 3): minimize f"(INT(x)!4)2 subject to x2)25 where INT(x) is the integer part of x. The design variable must also belong to the set M0·1, 0·2, 0·3, . . . , 9·9, 10·0N. 182 M.-W. HUANG AND J. S. ARORA Table II. Material data for spring design problem Material type 1 Material type 2 Material type 3 G (lb/in2) o (lb sec2/in4) q (lb/in2) ! U 1 11500000 12600000 13700000 0·000738342 0·000851211 0·000971362 80000 86000 87000 1·0 1·1 1·5 Figure 3. A non-differentiable function Problem 4: Three-bar space truss28 The objective of the problem is to minimize the weight of the structure. The 3 design variables (cross-sectional areas) must have one of the following allowable discrete values: M1·62, 1·80, 1·99, 2·13, 2·38, 2·62, 2·63, 2·88, 2·93, 3·09, 3·13, 3·38, 3·47, 3·55, 3·63, 3·84, 3·87, 3·88, 4·18, 4·22, 4·49, 4·59, 4·80, 4·97, 5·12, 5·74, 7·22, 7·97, 11·50, 13·50, 13·90, 14·20, 15·50, 16·00N. Three loading conditions are imposed. Stress constraint is considered for each member. Also displacement constraint is considered for the free node in two directions. Problem 5: Four-bar space truss29 There are 4 design variables which must have the same discrete values as in Problem 4. One loading condition is imposed. Stress and displacement constraints are considered. Problem 6: Design of a plate girder30 The problem is to design a plate girder for highway bridge. The dead load for the girder consists of weight of the pavement and self-weight of the girder. The live load with an impact factor of 1·33 consists of equivalent uniform load and a concentrated load based on HS20(MS18). The design variables are the height h, flange width b, flange thickness t and the web & thickness t . The objective is to minimize the total mass, or equivalently, the cross-sectional area 8 of the girder subject to the constraints on the bending stress, flange buckling, web crippling, deflection, shear stress, and size of the girder. The design variables must be from the following sets: h, b3M0·30, 0·32, 0·34, 0·36, . . . , 2·48, 2·50N t , t 3M0·010, 0·012, 0·014, 0·016, . . . , 0·098, 0·100N 8 & OPTIMAL DESIGN WITH DISCRETE VARIABLES 183 Problem 7: Mathematical problem31 This problem has seven design variables and four constraints. The first three variables must have discrete values. It is formulated as follows: minimize f"(x !10)2#5(x !12)2#x4#3(x !11)2#10x6 1 2 3 4 5 #7x2#x4!4x x !10x !8x 7 6 7 6 7 6 subject to g "2x2#3x4#x #4x2#5x )127 4 5 2 3 1 1 g "7x #3x #10x2#x !x )282 2 1 2 3 4 5 g "23x #x2#6x2!8x )196 6 7 2 3 1 g "4x2#x2!3x x #2x2#5x !11x )0 3 6 7 2 1 2 1 4 The first three design variables must be selected from the following sets: x 3M1, 2, 3, 4, 5N and x , x 3M0, 1, 2, 3, 4, 5N 1 2 3 Problem 8: 25-bar space tower-case 128 All design variables must have discrete values as for Problem 4. Two loading conditions are imposed. Stress, displacement and frequency constraints are considered. Problem 9: 25-bar space tower-case 232,33 This problem is the same as problem 8 except that the 25 members are grouped differently. Also different allowable discrete values are used: M0·1, 0·4, 0·7, 1·1, 1·5, 1·8, 2·0, 2·5, 3·0, 3·5, 4·0, 4·5N. Frequency constraint is not imposed. Problem 10: Cantilever beam34 A rectangular cantilever beam is subject to a point load of 50 000 N at its free end. The beam is 500 cm. long and is divided into 5 segments with different cross-section. The width (B) and depth (H) for each segment are defined as the design variables, giving a total of 10 design variables. The objective of this problem is to find the minimum volume subject to (1) stress constraints: maximum stress in each section must not exceed 14000 N/cm2; (2) depth/width ratio in each section must not exceed 20; and (3) the deflection at the free end must not exceed 2·7 cm. Young’s Modulus E is 200 GPa. All design variables must have integer values. Problem 11: 25-bar plane truss27 All design variables must have the same discrete values as in Problem 4. Two loading conditions are imposed. Stress and displacement constraints are considered. Problem 12: 49-bar plane truss-case 127 All design variables must have the same discrete values as in Problem 4. Two loading conditions are imposed. Stress constraint is not imposed. Problem 13: 49-bar plane truss-case 227 This problem is the same as problem 12 with one additional loading condition. Problem 14: 47-bar plane truss35 All design variables must have the same discrete values as in Problem 4. One loading condition is imposed. Stress constraint is considered. Problem 15: 200-bar tower28 This problem has 200 members divided into the following 96 groups. All design variables must have the same discrete values as in Problem 4. Three loading conditions are imposed for this structure. Stress constraint is considered for each member. Spring-1 (2, 3, 4) Spring-2 (5, 4, 4) Non-differentiable (3, 1, 1) 3-bar space truss (2, 3, 15) 4-bar space truss (2, 4, 7) Plate girder (2, 4, 5) Mathematical (2, 7, 4) 25-bar tower-1 (2, 7, 87) 25-bar tower-2 (8, 86) Cantilever beam (2, 10, 11) 25-bar plane truss (2, 13, 92) 49-bar truss-1 (2, 25, 98) 49-bar truss-2 (2, 25, 147) 47-bar truss (2, 47, 47) 200-bar tower (2, 96, 600) Problem name (PT, NV*, NC) Cost CPU Cost CPU Cost CPU Cost CPU Cost CPU Cost CPU Cost CPU Cost CPU Cost CPU Cost CPU COST CPU Cost CPU Cost CPU Cost CPU Cost CPU Simulated annealing 0·018 0·018 1·20 4·76 0·0156 0·016848 7·72 7·93 CH 0·0 —— 3·31 23·2164 23·2664 5·30 7·99 950·750 968·566 0·37 25·78 0·04524 0·04536 1·84 5·85 686·090 687·237 1·23 5·24 720·430 723·265 6·20 114·18 543·998 544·831 3·33 66·35 67800 69600 3·73 8·04 2183·98 2196·07 10·22 63·71 0·187053 0·194265 1·66 82·58 0·187053 0·194265 1·94 96·51 2412·69 2491·19 19·50 186·81 1·80589 2·20500 194·90 36980·28 Branch and bound 0·018 21·23 0·016848 18·63 0·0 1·99 23·2164 13·66 950·75 31·27 0·04532 27·61 686·956 28·59 732·343 67·15 544·831 66·19 69600 94·16 2310·62 2689·42 0·212634 91·72 0·212634 97·72 2665·41 517·14 2·37167 874·07 Genetic algorithm 0·01728s 0·31 CH —— CH —— 24·1649 0·30 950·750 0·27 0·04720 0·30 698·803 1·38 742·022 1·01 544·831 0·83 77000 0·68 2236·46 1·37 0·187053 1·95 0·187053 2·28 2412·69 5·91 NS —— Dynamic rounding-up NS 0·33 CH —— CH —— 24·0314 0·20 950·750 0·26 0·04680 0·22 686·090 0·61 NS 0·82 543·999 0·70 76200 1·03 2241·72 2·71 0·190912 7·61 0·187106 6·73 BIG —— BIG —— Neighbourhood search (NP"2) 0·0198 0·30 CH —— CH —— 23·9749 0·39 950·750 0·29 0·04604 0·34 686·090 8·68 NS 11·70 543·999 20·72 68700 498·37 BIG —— BIG —— BIG —— BIG —— BIG —— Neighbouhood search (NP"4) 0·012695 0·24 0·012695 0·35 CH —— 20·5435 0·19 676·460 0·20 0·0443517 0·12 683·981 0·33 590·723 0·73 534·807 0·18 63108·7 0·33 2167·90 0·69 0·15657 0·5 0·15451 0·6 1858·94 0·58 0·746626 55·17 Continuous solution PT: problem type, NV: number of design variables, NC: number of constraints, COST: cost function value, CPU: execution time (s), CH: cannot handle, NS: no solution found, and BIG: too many points * All design variables are discrete except for problem 7 in which only the first three design variables are discrete, others are continuous s With maximum constraint violation of 0·0313 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 No. Table III. Comparisons of the mixed variable optimization methods 184 M.-W. HUANG AND J. S. ARORA OPTIMAL DESIGN WITH DISCRETE VARIABLES 185 4.3. Numerical results Table III contains results for each of the 15 test problems with various algorithms. The presented data include cost function value (COST), and execution time (CPU). During the discrete optimization process, the BBM creates continuous subproblems (nodes) that are solved using the SQP method. The allowable constraint violation for continuous optimum solutions is set to 0·01 per cent for most problems and 1·0 or 0·1 per cent for some problems. The implemented BBM solves all Type 2 problems. As presented in Table III, this method always finds the solution with lowest cost among all test methods. Thus its results are also used to evaluate the performance of other methods. The BBM works best with small scale problems. For problems having many design variables, the number of subproblems created can be very large if a feasible discrete design cannot be found quickly (see Problems 11 and 15). A feasible discrete design serves as an upper bound for cost function value which eliminates many branches of the tree. For Problems 12 and 13, CPU is relatively small even with 25 design variables. This is due to the fact that many of the design variables reach their lower bounds at the continuous variable optimization stage, and need not be branched further. Note that the subproblem solver (the SQP method) of the implemented BBM requires differentiable functions. The functions of the second problem (Spring-2) are not differentiable with respect to one of the design variables: the material type. This problem is solved by the BBM incorporated with an explicit enumeration procedure. The third problem is also not solved since the cost function is not differentiable. For SA, the major stopping criterion is to quit the program if change in the cost function value is less than 0·0001 for 4 consecutive iterations. For the GA, the program stops if the change in cost function is less than 0·001 for 30 consecutive generations. These two methods do not require gradient information, and so the methods can handle problems other than Types 1 and 2. Table III contains results with the dynamic rounding-up method and two neighbourhood search methods that examine 2 and 4 discrete points near the continuous solution, respectively (NP"2 and NP"4). The dynamic rounding-up method may fail to find a feasible discrete solution. However, it is relatively fast when compared to other methods. The neighbourhood search method can be very effective for small problems if the best discrete solution is close to the continuous solution. However, for larger problems, when NP is set to 2 (or a larger number), the number of points to be examined become huge and large CPU time is needed. Both dynamic rounding-up and neighbourhood search method cannot handle problems other than Type 2, so the second and the third problems are not solved with these methods. The continuous solutions for all the test problems are also given in Table III for reference. The cost function for these solutions provides a lower bound for the cost function for the discrete solutions. It is also seen that the mixed variable optimization generally requires much more CPU time than the continuous variable optimization. Some algorithms, such as the BBM, obtain a continuous solution before going to the discrete optimization phase. Therefore, the CPU time required for such algorithms is expected to be much larger than those for the continuous problems. In summary, the implemented BBM performs better than other methods in the quality of the solution. The dynamic rounding-up method is the fastest method but the solution may not be as good or even feasible. The neighbourhood search method does not guarantee finding a feasible discrete solution if the number of examined neighboring points (NP) is small. However, there may be too many points to examine in a tolerable execution time if NP is set to a larger number. The SA method and GA usually need large execution times to find a global minimum. Although it is possible to find the best solution if temperature is reduced slowly and enough execution time is . 186 M.-W. HUANG AND J. S. ARORA allowed, the tested SA algorithm was terminated due to excessive amount of execution time for some of the test problems. This also points out a drawback of SA, that is, lack of an effective stopping criterion. It is difficult to tell whether a global or fairly good solution has been reached. It is also important to note that the CPU times for SA and GA can vary from one run to the next for the same problem. This is due to the random processes used in the algorithms. In addition, other parameters of the algorithm (rate of reduction of the temperature in SA; population size, and amounts of crossover and mutation) can also affect the total CPU time. 5. DISCUSSION AND CONCLUSIONS A general mixed variable non-linear optimization problem is defined. Six types of such problems are identified, each having its own characteristics and requirements. This classification covers the full range of non-linear problems, from continuous to purely discrete. The classification is important because it can guide in the selection and development of methods for each problem type. Several methods for solving various types of mixed variable problems are reviewed. The methods are cataloged depending on their features and the type of problems they can solve. Methods for solution of the six types of problems are identified. Five of the methods are implemented into a general purpose software. They are: the branchand-bound method, simulated annealing, a genetic algorithm, dynamic rounding-up, and neighbourhood search method. Many enhancements of the methods are proposed and implemented. Fifteen test problems are solved to assess performance of the methods—13 Type 2, 1 Type 3 and 1 Type 5. Based on the test results, the following conclusions are made for each implemented method. For convex problems, the BBM has a solid theoretical basis and a fathoming rule which guarantees a global minimum. For non-convex problems, it may end up fathoming nodes that should not have been fathomed. However, a useful solution can usually be obtained for such cases. The most obvious disadvantage of this method is the computational effort when the number of variables is large; however, with the proposed enhancements, it appears to be quite good and versatile. For problems having a reasonable number of design variables, the method is recommended. It also has an advantage that the number of possible discrete values for a variable does not significantly affect its performance. The SA and the GAs require many more function evaluations than other methods. However, they have some advantages: gradients of functions are not needed, convexity of the problem is not required for determining a global solution, and there is no subproblem to be solved. For large scale problems, if the analysis part can be speeded up (i.e., with the use of parallel computers and approximate methods), and the search space can be reduced, they can be faster than BBM. In comparing these two methods, the GA appears to use larger CPU time than the SA for most of the test problems. This is mostly because the implemented GA uses random access files to store various generations requiring more I/O which is inefficient. The algorithm actually uses fewer function evaluations than SA for most problems.27 This could be due to the reproduction process used in GA. Both methods generate random points around the existing point(s), but the reproduction process helps in speeding-up the search for better points without the need for function evaluations. Another advantage of GAs is that, since many trial points are generated simultaneously, their analyses can be divided into several parallel computing jobs. Note that, to avoid excessive CPU time, the stopping criteria used in the two methods are relaxed. As can be seen in Table III, these two methods never gave results that were better than the BBM. A much larger CPU time would be needed even if slightly better results were desired. Thus, the CPU times OPTIMAL DESIGN WITH DISCRETE VARIABLES 187 given in Table III may not be considered as performance benchmarks. They, however, do show that a larger CPU time is required by these two methods compared to the others. The main advantage of the dynamic rounding-up method is that the number of subproblems to be solved only grows linearly with the number of design variables. However, with a large number of design variable, selection of the variable for rounding-up becomes more difficult as some design constraints always tend to be violated. Thus, this method does not always find a feasible solution. Fortunately, for structural design problems such as truss design problems, rounding-up to a larger size for a member usually produces a feasible design. Thus, this method can be a good choice for some structural design problems; however, it needs to be fully developed and evaluated. The neighbourhood search method can be very effective for small problems. However, it is difficult to determine the range for the neighbouring points to be examined. The best solution may not be in the small neighbourhood, and, if more than two neighbouring points are examined for each variable of a large problem, the execution time can be very large. An advantage of this method and the dynamic rounding-up method, is that the number of subproblems to be solved can be estimated. Therefore, the required CPU time can also be estimated and can be used to judge the effectiveness of these methods. There are other techniques that were not selected for implementation. Zero—one variable technique is suitable for problems with alternative design selections; e.g., choosing from several types of materials. However, the problem will become large when the number of selections is large since each selection (alternative design) introduces its own set of design variables. For MDNLP, sequential linearization and integer programming techniques are not implemented for Type 2 problems. Based on the study, it is concluded that mixed variable non-linear optimization problems require substantial more computational effort compared to the continuous variable optimization problems. In addition, optimality of the final solution is not ensured because no optimality conditions can be verified. It is also difficult to specify a precise stopping criterion for such problems. This study does not intend to suggest a particular method over others for solving a wide variety of problems. It, however, suggests that the problem size and problem type be first identified and then a suitable method be chosen for its solution. The study also suggests that more comparative evaluation of the numerical methods for discrete variable problems needs to be carried out to better understand the methods so that they can be developed further. It is important to note that the use of parallel processing capabilities can reduce the CPU times dramatically for most discrete optimization algorithms, such as branch and bound, simulated annealing and genetic algorithms. This is particularly true for larger scale problems. ACKNOWLEDGMENTS This paper is based on a presentation made by the authors at the 36th AIAA/ASME/ASCE/AHS Structures, Structural Dynamics and Materials Conference, April 10—12, 1995, New Orleans, LA. The authors acknowledge partial support for this research under the project, ‘Design Sensitivity Analysis and Optimization of Non-linear Structures’, from the NSF of U.S.A., Grant No. CMS-93-01580. We are also thankful to the reviewers for providing us with useful comments and suggestions. REFERENCES 1. J. S. Arora, M. W. Huang and C. C. Hsieh, ‘Methods for optimization of nonlinear problems with discrete variables: a review’, Struct. Optim., 8, 69—85 (1994). 2. J. S. Arora, T. C. Lin, O. Elwakeil and M. W. Huang, IDESIGN ºser’s Manual »ersion 4.2, Optimal Design Laboratory, College of Engineering, The University of Iowa, Iowa City, Iowa 52242 U.S.A., 1994. 188 M.-W. HUANG AND J. S. ARORA 3. E. Sandgren, ‘Nonlinear integer and discrete programming for topological decision making in engineering design’, J. Mech. Des. ASME, 112, 118—122 (1990b). 4. R. K. Kincaid and S. L. Padula, ‘Minimizing distortion and internal forces in truss structures by simulated annealing’, Proc. 31st AIAA SDM Conf., Long Beach, CA, 1990, pp. 327—333. 5. P. B. Thanedar and G. N. Vanderplaats, ‘A survey of discrete variable optimization for structural design’, J. Struct. Eng. ASCE, 121, 301—306 (1994). 6. H. T. Loh and P. Y. Papalambros, ‘A sequential linearization approach for solving mixed—discrete nonlinear design optimization problems’, J. Mech. Des. ASME, 113, 325—334 (1991a). 7. H. T. Loh and P. Y. Papalambros, ‘Computational implementation and tests of a sequential linearization approach for solving mixed—discrete nonlinear design optimization’, J. Mech. Des. ASME, 113, 335—345 (1991b). 8. M. Bremicker, P. Y. Papalambros and H. T. Loh, ‘Solution of mixed-discrete structural optimization problems with a new sequential linearization algorithm’, Comput. Struct., 37, 451—461 (1990). 9. E. Sandgren, ‘Nonlinear integer and discrete programming in mechanical design optimization’, J. Mech. Des. ASME, 112, 223—229 (1990a). 10. J. S. Arora and M. W. Huang, ‘Discrete structural optimization with commercially available sections: a review’, J. Struct. Earthquake Eng., Japan Society of Civil Engineers (1996) to appear. 11. J. J. Moré and S. J. Wright, Optimization Software Guide, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1993. 12. K. Hager and R. J. Balling, ‘New approach for discrete structural optimization’, J. Struct. Eng. ASCE, 114, 1120—1134 (1988). 13. C. H. Tseng, L. W. Wang and S. F. Ling, ‘Enhancing branch-and-bound method for structural optimization’, J. Struct. Eng. ASCE, 121, 831—837 (1995). 14. U. T. Ringertz, ‘On methods for discrete structural optimization’, Eng. Optim., 13, 47—64 (1988). 15. G. L. Nemhauser and S. J. Wolsey, Integer and Combinatorial Optimization, Wiley, New York, N.Y., 1988. 16. E. Aarts and J. Korst, Simulated Annealing and Boltzmann Machines: A Stochastic Approach to Combinatorial Optimization and Neural Computing, Wiley-Interscience Series in Discrete Mathematics and Optimization, Wiley, New York, 1989. 17. R. J. Balling, ‘Optimal steel frame design by simulated annealing’, J. Struct. Eng. ASCE, 117, 1780—1795 (1991). 18. S. A. May and R. J. Balling, ‘A filtered simulated annealing strategy for discrete optimization of 3D steel frameworks’, Struct. Optim., 4, 142—148 (1992). 19. J. H. Holland, Adaptation in Natural and Artificial System, The University of Michigan Press, Ann Arbor, MI, 1975. 20. D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine ¸earning, Addison-Wesley, Reading, MA, U.S.A., 1989. 21. J. S. Arora, Introduction to Optimal Design, McGraw-Hill, New York, 1989. 22. C.-M. Chan, D. E. Grierson and A. N. Sherbourne, ‘Automatic optimal design of tall steel building frameworks’, J. Struct. Eng. ASCE, 121, 838—847 (1995). 23. M. W. Huang and J. S. Arora, ‘A self-scaling implicit SQP method for large scale structural optimization’, Int. j. numer. methods eng., 39, 1933—1953 (1996a). 24. M. W. Huang and J. S. Arora, ‘Design of steel structures using standard sections’, ¹echnical Report No. OD¸-96.04, Department of Civil and Environmental Engineering, The University of Iowa, Iowa City, IA 52242, 1996b. 25. M. W. Huang, C. C. Hsieh and J. S. Arora, ‘A genetic algorithm for sequencing type problems in engineering design’, ¹echnical Report No. OD¸-96.05, Department of Civil and Environmental Engineering, The University of Iowa, Iowa City, IA 52242, 1996. 26. F. Y. Kocer and J. S. Arora, ‘Design of prestressed concrete poles: an optimization approach’, J. Struct. Eng. ASCE, (1996), to appear. 27. M. W. Huang, ‘Algorithms for mixed continuous-discrete variable problems in structural optimization’, Ph.D dissertation, Civil and Environmental Engineering, The University of Iowa, Iowa City, IA, U.S.A., 1995. 28. E. J. Haug and J. S. Arora, Applied Optimal Design, Wiley-Interscience, New York, 1979. 29. V. B. Venkayya, ‘Design of optimum structure’, Comput. Struct., 1, 265—309 (1971). 30. J. S. Arora, ‘What is structural optimization’, Proc. Structures Congr. ’93, ASCE, Irvine, CA, 1993. 31. W. Hock and K. Schittkowski, ‘Test examples for nonlinear programming codes’, Lecture Notes in Economics and Mathematical Systems 187, Springer, Berlin Heidelberg New York, 1980. 32. G. R. Olsen and G. N. Vanderplaats, ‘Method for nonlinear optimization with discrete design variables’, AIAA J., 27, 1584—1589 (1989). 33. D. K. Shin, Z. Gürdal and O. H. Griffin Jr., ‘A penalty approach for nonlinear optimization with discrete design variables’, Eng. Optim., 16, 29—42 (1990). 34. G. N. Vanderplaats, Numerical Optimization ¹echniques for Engineering Design: ¼ith Applications, McGraw-Hill New York, 1984. 35. E. Salajegheh and G. N. Vanderplaats, ‘Optimum design of trusses with discrete sizing and shape variables’, Struct. Optim., 6, 79—85 (1993). .