Parallel optimization algorithms for high performance computing : application to thermal systems

UPC CTTC Parallel optimization algorithms for High Performance Computing. Application to thermal systems. Centre Tecnològic de Transferència de Calor Departament de Màquines i Motors Tèrmics Universitat Politècnica de Catalunya Imanol Aizpurua Udabe Doctoral Thesis Parallel optimization algorithms for High Performance Computing. Application to thermal systems. Imanol Aizpurua Udabe TESI DOCTORAL presentada al Departament de Màquines i Motors Tèrmics E.S.E.I.A.A.T. Universitat Politècnica de Catalunya per a l’obtenció del grau de Doctor per la Universitat Politècnica de Catalunya Terrassa, February 28, 2017 Parallel optimization algorithms for High Performance Computing. Application to thermal systems. Imanol Aizpurua Udabe Director de la tesi Dr. Assensi Oliva Llena Tribunal Qualificador Dr. Antonio Pascau Benito Universidad de Zaragoza Dr. Cristóbal Cortés Gracia Universidad de Zaragoza Dr. Francesc Xavier Trias Miquel Universitat Politècnica de Catalunya To my family. Wayfarer, the only way is your footprints and no other. Wayfarer, there is no way. Make your way by going farther. By going farther, make your way till looking back at where you’ve wandered, you look back on that path you may not set foot on from now onward. Wayfarer, there is no way; only wake-trails on the waters. Antonio Machado (translated by A. Z. Foreman) i A mi familia. Caminante, son tus huellas el camino y nada más; caminante, no hay camino, se hace camino al andar. Al andar se hace camino y al volver la vista atrás se ve la senda que nunca se ha de volver a pisar. Caminante no hay camino sino estelas en la mar... Antonio Machado iii x Contents 2 Implementation of a new optimization library: Optimus 57 2.1 Introduction ..................................... 58 2.1.1 General design requirements of a library ............... 58 2.1.2 Specific design requirements of the new optimization library . . . 59 2.1.3 Computing facilities at CTTC ...................... 60 2.1.4 Concluding remarks ............................ 61 2.2 State of the art of optimization libraries .................... 62 2.2.1 Description of the project’s needs .................... 62 2.2.2 Free open-source software ........................ 64 2.2.3 Proprietary commercial software .................... 66 2.2.4 Concluding remarks ............................ 67 2.3 Main features of Optimus ............................. 69 2.3.1 Development strategy ........................... 69 2.3.2 Definition of the optimization problem ................. 69 2.3.3 Single-objective vs. Multi-objective optimization ........... 71 2.3.4 Genetic operators ............................. 73 2.3.5 Hybrid methods .............................. 78 2.3.6 Continuation criteria ........................... 79 2.3.7 Statistics .................................. 81 2.3.8 Parallelization ............................... 82 2.3.9 User interface ............................... 88 2.3.10 Other features ............................... 91 2.3.11 Optimus vs. Paradiseo .......................... 92 2.4 Validation tests ................................... 92 2.4.1 Benchmark mathematical functions .................. 93 2.4.2 CFD & HT tests .............................. 96 2.5 Conclusions ..................................... 110 References ...................................... 111 3 Load balancing methods for parallel optimization algorithms 115 3.1 Introduction ..................................... 116 3.2 Approach to the load balancing problem .................... 118 3.2.1 Definitions ................................. 118 3.2.2 Factors affecting parallel performance ................. 124 3.2.3 Applications using task schedulers ................... 127 3.2.4 Methodology for developing load balancing algorithms ....... 128 Contents xi 3.2.5 State of the art of algorithms for solving the combinatorial scheduling problem ................................. 132 3.2.6 State of the art of time estimation techniques ............ 138 3.3 Load balancing strategies ............................. 142 3.3.1 Overview .................................. 142 3.3.2 Task management algorithms ...................... 149 3.3.3 Task scheduling algorithm ........................ 166 3.3.4 Task assignment algorithm ....................... 172 3.4 Theoretical case study of load balancing strategies ............. 179 3.4.1 Design of the experiments ........................ 179 3.4.2 Short cases with linear scalability ................... 186 3.4.3 Long cases with linear scalability .................... 198 3.4.4 Long cases with non-linear scalability ................. 206 3.4.5 Concluding remarks ............................ 216 3.5 Implementation and testing of time estimation techniques ......... 218 3.6 Illustrative example ................................ 221 3.7 Conclusions ..................................... 232 References ...................................... 233 4 Conclusions and future work 237 4.1 Conclusions ..................................... 238 4.2 Future work ..................................... 241 xii Contents Abstract The need of optimization is present in every field of engineering. Moreover, applications requiring a multidisciplinary approach in order to make a step forward are increasing. This leads to the need of solving complex optimization problems that exceed the capacity of human brain or intuition. A standard way of proceeding is to use evolutionary algorithms, among which genetic algorithms hold a prominent place. These are characterized by their robustness and versatility, as well as their high computational cost and low convergence speed. Such drawbacks are usually tackled by hybridizing them with local search methods, e.g. gradient methods, in order to obtain significant speed up. The use of multiple levels of fidelity of the objective functions is also a common practice. Many optimization packages are available under free software licenses and are representative of the current state of the art in optimization technology: single-objective and multi-objective search techniques, global and local search methods, plenty of genetic operators, mixed integer optimization capacity, fitness landscape analysis techniques, parallelization strategies, etc. However, the ability of optimization algorithms to adapt to massively parallel computers reaching satisfactory efficiency levels is still an open issue. Even packages suited for multilevel parallelism encounter difficulties when dealing with objective functions involving long and variable simulation times. This variability is common in Computational Fluid Dynamics and Heat Transfer (CFD & HT), nonlinear mechanics, etc. and is nowadays a dominant concern for large scale applications. Current research in improving the performance of evolutionary algorithms is mainly focused on developing new search algorithms. Nevertheless, there is a vast knowledge of sequential well-performing algorithmic suitable for being implemented in parallel computers. The gap to be covered is efficient parallelization. Moreover, advances in the research of both new search algorithms and efficient parallelization are additive, so that the enhancement of current state of the art optimization software can be accelerated if xiii xiv Abstract both fronts are tackled simultaneously. The motivation of this Doctoral Thesis is to make a step forward towards the successful integration of Optimization and High Performance Computing capabilities, which has the potential to boost technological development by providing better designs, shortening product development times and minimizing the required resources. A generic mathematical optimization tool has been developed for this aim, applicable in any field of science and engineering. Nevertheless, being this research activity hosted by the Heat and Mass Transfer Technological Center (CTTC), a special focus has been put on the application of the library to the fields of expertise of the Center: Computational Fluid Dynamics and Heat Transfer (CFD & HT), multi-physics simulation, etc. This document is structured in four chapters. A thorough state of the art study is conducted in the first chapter with the aim of obtaining a global scope of the mathematical optimization techniques available to date, as well as their most remarkable virtues and shortcomings. After classifying optimization problems according to their principal characteristics, such as the number of objective functions (single-objective vs multi-objective) or the nature of the optimization variables (real vs discrete), an insight of constraint handling techniques, surrogate-based optimization and hybrid optimization methods is provided. The most widespread global search algorithm is then introduced, namely the genetic algorithm, followed by a description of the main concepts and shortcomings of the standard parallelization strategies available for such population-based optimization methods. Finally, an overview of the most common test suites for optimization algorithms is given together with some remarks related to random number generators. The second chapter explains how the implementation of the new optimization library Optimus has been carried out based on the research on optimization theory conducted in the first chapter. The first step has been the definition of the design requirements of the library, including both general code development patterns and specific requirements for the optimization tool. A state of the art study of currently available optimization libraries is then included, taking into consideration open-source and proprietary software. Once the development strategy of the new library is fixed, the main features of Optimus are introduced. Finally, several validation tests are performed in order to demonstrate the suitability of the new library for solving benchmark mathematical optimization tests and real-world CFD & HT optimization problems. The third chapter contains the main contribution of this Doctoral Thesis. It is started with an approach to the computational load balancing problem detected in the first chapter when the state of the art study on the parallelization of genetic Abstract xv and other population-based optimization algorithms was carried out. The core of the problem is that processors are often unable to finish the evaluation of their queue of individuals simultaneously and need to be synchronized before the next batch of individuals is created. Consequently, the computational load imbalance is translated into idle time in some processors. This fact was identified as the key point causing the degradation of the optimization algorithm’s scalability (i.e. parallel efficiency) in case the average makespan of the batch of individuals is greater than the average time required by the optimizer for performing inter-processor communications. According to the methodology defined for developing load balancing algorithms, the load balancing problem is split into two sub-problems: the estimation of the time required to process each individual, and the subsequent resolution of a combinatorial task scheduling problem in order to map tasks to processors with a certain precedence relation and in the most efficient manner with the aim of reducing the evaluation makespan of each batch of individuals. Several load balancing algorithms are proposed and exhaustively tested by means of 3 theoretical case studies using the genetic algorithm. Being the latter the most widespread optimization heuristic, the impact of the research is expected to be maximized. Note that the proposed algorithms and the reached conclusions are extendable to any other population-based optimization method that needs to synchronize all processors after the evaluation of each batch of individuals. Since time estimation techniques have not been studied in detail due to lack of time, the availability of perfect individuals’ evaluation time estimations has been assumed. However, a first implementation of time estimation techniques together with some preliminary tests is included towards the end of the chapter. Finally, a real-world engineering application that consists on optimizing the refrigeration system of a power electronic device is presented as an illustrative example in which the use of the proposed load balancing algorithms is able to reduce the simulation time required by the optimization tool. The fourth chapter gathers the main conclusions of the conducted research and outlines the next steps to be followed in this approach towards the integration of Optimization Techniques and High Performance Computing. xvi Abstract 1 State of the art of optimization algorithms Abstract. A thorough state of the art study is conducted in this first chapter with the aim of obtaining a global scope of the mathematical optimization techniques available to date, as well as of their most remarkable virtues and shortcomings. Being the application field of interest that of Computational Fluid Dynamics & Heat Transfer (CFD & HT), which is characterized by the simulation of computationally expensive non-linear equation systems, special emphasis is put on the parallelizability of the optimization algorithms. 1 2 §1.1 Introduction 1.1 Introduction Every human activity is characterized by the search of Best: sport, social relations, work. . . There are plenty of things in our daily life which we try to optimize, such as the time spent going from home to our work place, the time spent cooking, the budget for our holidays, our personal appearance, etc. However, perfection is a quite philosophical concept. Although it guides human actions, reality cannot be understood without constraints. In real life, we never have the option of satisfying all our wishes at the same time: getting closer to perfection in some aspect pushes us to imperfection in some other. Hence, finding a compromise solution among all real options becomes our only feasible goal. In this sense, optimization could be defined as the art of making the most of something rather than the search for perfection. Once the labyrinth formed by the objective (or objectives) and the constraints is defined, the next question arises: how to get to the goal? We know our starting point, we know our destination and we know the rules of what is feasible and infeasible. Assuming an explorer’s role, we lack the map and the compass that can guide us on the way. This orientation technique is what Optimization Theory takes care of. Optimization is the technique that brings explorers to their destination through the shortest way, avoiding that they ran out of resources before reaching their goal. Using a more rigorous language, it can be said that Optimization Theory encompasses the quantitative study of optima and methods for finding them [1]. Given life’s and world’s complexity, optimization techniques require a thorough study. The search for the right way to the goal is not easy, and the explorer will encounter plenty of “distractions” that may mislead him in his mission. He might think his search has finished, just because he has found some interesting place during the journey. Usually the goal is only one, and the explorer’s orientation technique should be sophisticated enough to let him know the improvability of his findings, until he reaches the final Goal. Leaving the explorer’s analogy aside, but following the same reasoning, let’s move to the field of engineering and applied mathematics. In any design process, we have a device, mathematical process or experiment which we want to optimize according to some criteria. The optimality measure that will allow us to say if a solution is better than another one is the objective function , also called cost function or fitness function . As previously said, there may be one or multiple objective functions, and the engineer will want to maximize or minimize them. The output of the objective function is dependent on its inputs, which correspond to the characteristics of the §1.1 Introduction 3 device, experiment, etc. These inputs will be called optimization variables from now on. Moreover, the optimization problem may be constrained, meaning this that each optimization variable cannot be modified independently from others. As stated in [2], several categories may be distinguished in optimization: • Single-variable / multi-variable: Optimization becomes increasingly difficult as the number of dimensions (number of optimization variables) increases. • Static / dynamic: Dynamic optimization means that the output is a function of time, while static means that the output is independent of time. • Constrained / unconstrained: Constrained optimization incorporates variable equalities and inequalities into the objective function, whereas unconstrained optimization allows the variables to take any value. A constrained variable may often be converted into an unconstrained variable through a mathematical transformation. Consider the simple constrained example of minimizing f ( x ) over the interval − 1 ≤x≤ 1. The variable x may be converted into an unconstrained variable uby letting x=sin(u) for any value of u. • Discrete / continuous: Discrete variables (combinatorial optimization) have only a finite number of possible values, whereas continuous variables have an infinite number of possible values. • Single-objective / multi-objective: Single-objective optimizations seek to improve a unique objective, whereas multi-objective optimizations search for a compromise optimal solution for a set of objectives. Since usually a fitness surface has many peaks, valleys and ridges (see Fig. 1.1), a difficulty with optimization is to determine if a given optimum (minimum or maximum) is the global optimum or a local optimum (suboptimal value). In this sense, some optimization methods have explorative nature, whereas other methods have exploitative nature. Exploration refers to the ability of keeping a global view of the solutions space, and allows distinguishing the global or local nature of an optimum. However, a drawback of explorative methods is their low performance when converging to the optimum in the neighborhood of a known solution. On the other hand, exploitation refers to the ability of converging fast and accurately to a local optimum starting from some solution in the neighborhood of that optimum. Nevertheless, the limitation of exploitative algorithms is their incapacity to know if the solution found is a local or global optimum. 10 §1.2 Single-objective vs. Multi-objective optimization maker”. Thus, the goal is to optimize (maximize or minimize) k objective functions simultaneously. A global MOP problem can be formally defined as in [31]: Definition 1.2 General Multi-Objective Optimization Problem (MOP): A general MOP is defined as minimizing (or maximizing) F ( x ) = ( f1 ( x ) ,..., fk ( x )) subject to gi ( x ) ≤ 0, i={ 1 ,...,m} , and hj ( x ) = 0, j={ 1 ,..., p}x∈Ω . An MOP solution minimizes (or maximizes) the components of a vector F ( x )where x is a n-dimensional decision variable vector x= ( x1,..., xn )from some universe Ω . It is noted that gi ( x ) ≤ 0and hj ( x ) = 0 represent constraints that must be fulfilled while minimizing (or maximizing) F ( x )and Ωcontains all possible xthat can be used to satisfy an evaluation of F(x). The main conceptual difference between the single-objective and multi-objective optimizations is the difficulty when comparing two possible solutions in the latter case. A single-objective optimization seeks to maximize or minimize a unique objective, so two objective values can be easily compared in order to decide which solution is the best. On the other hand, and following the same logic, when more than one objective is present a solution may beat another one according to some objectives, but not according to others. If such a situation arises, how should a decision be made? In order to answer this question, three important concepts of multi-objective optimization are introduced next: fitness assignment, diversity preservation and elitism. Fitness assignment Fitness assignment schemes may be classified into four different categories [32]: •Scalar approaches , where the MOP is reduced to a single-objective optimization problem, for instance by means of a weighted-sum aggregation. •Criterion-based approaches , where each objective function is treated separately. In VEGA (Vector Evaluated Genetic Algorithm) [33], for instance, a parallel selection is performed where solutions are discerned according to their values on a single objective function, independently from the others. •Dominance-based approaches , where a dominance relation is used to classify solutions. The main techniques are i) dominance-rank techniques, which compute the number of population items that dominate a given solution, ii) dominancecount techniques, where the fitness value of a solution corresponds to the number of individuals that are dominated by this solution, and iii) dominance-depth §1.2 Single-objective vs. Multi-objective optimization 11 strategies, which classify a set of solutions into different classes or fronts. Hence, a solution that belongs to a class does not dominate another one from the same class. •Indicator-based approaches , where the fitness values are computed by comparing individuals on the basis of a quality indicator I . The chosen indicator represents the overall goal of the search process. Examples of indicator-based EAs are IBEA (Indicator-Based EA) [34] or SMS-EMOA (S-Metric Selection Evolutionary Multi-objective Optimization Algorithm) [35]. The scalar and criterion-based approaches are the most simplistic strategies, being dominance-based and indicator-based approaches usually preferred. These two techniques provide a set of optimal solutions, not just a single solution. But no matter the selected strategy, some criterion is needed in order to select a set of solutions rather than another one. Dominance-based approaches commonly use the so called Paretodominance criterion (see [36] for more details on Pareto’s Optimality Theory), although some new techniques appeared recently based on other dominance operators such as ε-dominance [37] or g-dominance [38]. Some definitions related to Pareto’s Optimality Theory are introduced next, assuming a multi-objective minimization problem: Definition 1.3 Pareto Optimality [31]: A solution x∈Ω is said to be Pareto Optimal with respect to (w.r.t.) Ω if and only if (iff) there is no x’ ∈Ω for which v=F ( x’ ) = ( f1 ( x’ ) ,..., fk ( x’ )) dominates u=F ( x ) = ( f1 ( x ) ,..., fk ( x )). The phrase Pareto Optimal is taken to mean with respect to the entire decision variable space unless otherwise specified. In other words, x* is Pareto optimal if there exists no feasible vector x which would decrease some criterion without causing a simultaneous increase in at least one other criterion. Definition 1.4 Pareto Dominance [31]: A vector u= ( u1,...,uk )is said to dominate another vector v= ( v1,...,vk )(denoted by u¹v ) if and only if u is partially less than v, i.e., ∀i∈{1,...,k},ui≤vi∧∃i∈{1,...,k}:ui<vi. Definition 1.5 Pareto Optimal Set [31]: For a given MOP, F ( x ), the Pareto Optimal Set, P∗, is defined as: 12 §1.2 Single-objective vs. Multi-objective optimization P∗:={x∈Ω| ¬∃ x’ ∈ΩF(x’)¹F(x)} Summarizing, Pareto optimal solutions are those solutions within the decision space whose corresponding objective vector components cannot be all simultaneously improved. When plotted in the objective space, the non-dominated vectors are collectively known as the Pareto front. Although single-objective optimization problems may have a unique optimal solution, MOPs usually have a possibly uncountable set of solutions on a Pareto front. Each solution associated with a point on the Pareto front is a vector whose components represent trade-offs in the decision space. This is the reason why defining an MOP’s global optimum is not a trivial task as the “best” compromise solution depends on the preferences of the (human) decision maker. Thus, the Pareto front represents the “best” solutions available and allows the definition of an MOP’s global optimum. Diversity preservation Approximating the Pareto optimal set is not only a question of convergence. For a good characterization of the true optimal Pareto front, the final approximation must also be well spread over the objective space. Therefore, a diversity preservation mechanism is usually integrated into the algorithm to uniformly distribute the population over the trade-off surface. Popular examples of Evolutionary Multi-objective Optimization (EMO) diversity assignment techniques are sharing and crowding [32]. Sharing consists on estimating the distribution density of a solution using a so-called sharing function that is related to the sum of distances to its neighborhood solutions. Crowding allows maintaining diversity without specifying any parameter. It consists in estimating the density of solutions surrounding a particular point of the objective space. Elitism Elitism [32] consists on maintaining an external set (archive) that allows storing either all or a subset of non-dominated solutions found during the search process. This secondary population aims at preventing the loss of these solutions during the stochastic optimization process, and is continuously updated with new potential non-dominated solutions. Even if an archive is usually used as an external storage only, archive members can also be integrated during the selection phase of an EMO algorithm, leading to elitist EMOs. Keeping in mind the concepts introduced so far, and being a Pareto front the solution provided by the optimization algorithm to the engineer, the main goals of a multiobjective optimization algorithm may be established as follows [32]: §1.2 Single-objective vs. Multi-objective optimization 13 • Convergence, i.e. the distance of the resulting non-dominated set to the true Pareto front should be minimized. • Diversity, i.e. the Pareto front must be well characterized in its whole length, avoiding a high concentration of solutions in one area and lacking solutions in another area. • Elitism, i.e. non-dominated points in the objective space and associated solution points in the decision space must be preserved. Several ways have been proposed in order to achieve the goals mentioned above. The following classification of multi-objective optimization algorithms, very similar to the one proposed for single-objective methods, gathers the most widely spread algorithms to date: •Enumerative algorithms: This is the family of the inefficient brute-force algorithms, such as the exhaustive search. •Deterministic algorithms: The range of available algorithms for multi-objective optimization is more limited than for single-objective optimization, or the use of such methods is at least less spread [5]. Depth-first search (hill-climbing) [39] and gradient-based search methods could be mentioned as the most usual strategies. The latter are only applicable to MOPs consisting of continuous variables, being the derivatives obtained based upon a specific direction from selected non-dominated points in the known Pareto front with the aim of moving a point towards the true Pareto front. Note however that derivatives tend to be quite noisy in such situations. •Stochastic algorithms: Natural optimization methods (mainly evolutionary algorithms [5]) have remarkably succeeded in solving multi-objective optimization problems. A key advantage has been that the population-based nature of EAs allows the generation of several elements of the Pareto optimal set in a single run. This field is now called Evolutionary Multi-objective Optimization (EMO), which refers to the use of evolutionary algorithms of any sort. The most notorious algorithms currently available are the following: MOGA [40], NPGA [41], VEGA [33], NSGA [42], PAES [43], NSGA-II [44], SPEA [45], SPEA2 [46], PESA [47], ε-MOEA [37,48] and IBEA [34]. They use techniques going from a simple linear aggregating function to the most popular Multi-Objective Evolutionary 14 §1.2 Single-objective vs. Multi-objective optimization Algorithms (MOEAs) based on Pareto ranking. Other natural optimization algorithms that are used for multi-objective optimization are the multi-objective particle swarm optimization (MOPSO) , cultural algorithms , differential evolution , predator-prey algorithm and ant colony optimization . Some additional strategies that do not belong to natural optimization methods are also used, e.g. simulated annealing [6,20] and tabu search [19], but there seems to be a clear preference for the use of evolutionary algorithms. A brief overview of some stochastic algorithms mentioned above is provided hereafter (based on [49]): •Pareto ranking (MOGA): Fonseca and Fleming [40] proposed a variation of Goldberg’s fitness assignment where a solution’s rank corresponds to the number of solutions in the current population by which it is dominated. •Pareto sharing: A fitness assignment like the previous one tends to produce premature convergence, what does not guarantee a uniformly sampled final Pareto approximation set. To avoid that, Fonseca and Fleming [40] modified the strategy above by implementing fitness sharing in the objective space to distribute the population over the Pareto-optimal region. •Non-dominated Sorting Genetic Algorithm (NSGA): Srinivas and Deb [42] introduced this algorithm which classifies the solutions into several classes (or fronts). A solution that belongs to a class does not dominate another one from the same class. Logically, the best fitness value is assigned to solutions of the first class, because they are closest to the true Pareto-optimal front of the problem. Diversity is preserved by means of a fitness sharing procedure. •NSGA-II: This is a modified version of NSGA introduced by Deb et al. [44]. The algorithm is computationally more efficient, uses elitism and keeps diversity by means of a crowding technique. •Strength Pareto Evolutionary Algorithm (SPEA): This elitist algorithm was proposed by Zitzler and Thiele [45]. It maintains an external population (an archive) that stores a fixed number of non-dominated solutions found during the optimization process in order to define the fitness of a solution based on these archive members. §1.2 Single-objective vs. Multi-objective optimization 15 •SPEA2: It is an improved version of SPEA, introduced by Zitzler et al. [46]. In comparison to its predecessor, SPEA2 includes an improved fitness assignment technique, a density estimation technique and an archive truncation method. •Indicator-Based Evolutionary Algorithm (IBEA): Introduced by Zitzler and Künzli [34], it has the characteristic to compute fitness values by comparing individuals on the basis of an arbitrary binary quality indicator I (also called binary performance metric). Thereby, no particular diversification mechanism is necessary. The indicator, determined according to the decision maker preferences, denotes the overall goal of the optimization process. Two binary quality indicators commonly used are the additive ε -indicator [50] and the IHD -indicator [50] that is based on the hypervolume concept [45] (see [50] for an overview about quality indicators). IBEA is a good illustration of the new EMO trend dealing with indicator-based search that started to become popular in recent years. Once the goals for multi-objective optimization methods have been established, some performance assessment method is to be developed in order to know which algorithm works best. The existing performance metrics can be classified into three classes [51]: • Convergence metrics: They evaluate how far the known Pareto front is with respect to the true Pareto front. • Diversity metrics: They evaluate how scattered the final population of the Pareto front is. • Metrics for both convergence and diversity: They evaluate both the distance to the true Pareto front and the dispersion of final solutions. Various quality indicators have been proposed in the literature for evaluating the performance of multi-objective search methods [32,50,52 – 54]: entropy, contribution, generational distance, spacing, coverage of two sets, coverage difference, S-metric, Dmetric, R-metrics, hypervolume metric (both in unary and binary form), additive and multiplicative ε -indicators, etc. However, none of the metrics can be considered the best, being usually necessary to use more than one metric to evaluate the performance of the multi-objective evolutionary algorithms. The reader is referred to [50] for a general review. 16 §1.4 Surrogate-based optimization 1.3 Constraint handling Realistic engineering problems are always multidisciplinary. Consequently, constraints of both equality and inequality type are very likely to appear. A set of design variables that does not violate any constraints is said to be feasible, while design variables that violate one or more constraints are infeasible. Every general purpose optimization algorithm must be able to handle the existence of these equality and inequality constraints. Based on [55], the following are the most common ways and their use makes sense depending on the nature of the algorithm (gradient based, non-gradient based, etc.): •Restauration method (also called feasible search ) : Designs that violate constraints are automatically restored to feasibility via the minimization of the active global constraint functions. •Rosen’s projection method: Provides search directions which guide the descent direction tangent to active constraint boundaries. •Penalty method: The fitness of the designs that violate constraints is artificially worsened, with the aim of forcing the optimization algorithm to abandon that search region of the solutions space. •Random design generation: When an infeasible design is detected, it is discarded and new random designs are generated within a (for instance, Gaussianshaped) probability density cloud about a desirable and feasible design until a new design is reached. The constraint handling techniques must be carefully chosen and implemented, because they have a direct effect on the convergence of the optimization algorithm. 1.4 Surrogate-based optimization One of the main concerns when running an optimization process is the computational cost. The evaluation of real-life objective functions can be very expensive, since it might involve solving non-linear systems, big meshes, etc. Taking into account that the optimization could need hundreds or thousands of objective function evaluations, the need for computational resources may seem scary. Therefore, for problems where objective function evaluations are already expensive and where the number of design variables is large (thus requiring many objective function evaluations), the only economically §1.4 Surrogate-based optimization 17 viable approach to optimization is to use an inexpensive and as accurate as possible surrogate model (a metamodel) instead of the actual high fidelity analysis method. Such surrogate models are known as response surfaces. They are very useful at the early stages of optimization, although progressively more complete physical models should be used as the global optimization process starts converging. Surrogate models are fitted through the available (often small) set of high fidelity values of the objective function. Once the response surface is created using an appropriate analytic formulation, it is very easy and fast to search such a surface for its optimum given a set of values of design variables supporting such a response surface. Some basic concepts related to the response surface generation methodology are presented in [6]. From the viewpoint of kernel interpolation/approximation techniques, many response surface methods are based on linear and non-linear regression and other variants of the least square technique. This group of mesh-free methods has been successfully applied to many practical, but difficult problems in engineering that are to be solved by the traditional mesh-based methods. The commercial optimization software IOSO (Indirect Optimization Based Upon Self-Organization) [56], a software known for its extraordinary speed and robustness, partly owes its success to the appropriate use of response surfaces. Due to the existence of several response surface techniques, their performance is to be evaluated according to some criterion. Here are some key aspects worth taking into account proposed by [57]: • Accuracy: The capability of predicting the system response over the design space of interest. • Robustness: The capability of achieving good accuracy for different problem types and sample sizes. • Efficiency: The computational effort required for constructing the metamodel and for predicting the response for a set of new points by metamodels. • Transparency: The capability of illustrating explicit relationships between input variables and responses. • Conceptual simplicity: Ease of implementation. Simple methods should require minimum user input and be easily adapted to each problem. A crucial aspect for the proper construction of a response surface by means of any algorithm is the location of the training points. If we are given freedom to choose 18 §1.5 Hybrid optimization methods the locations of the support points of a multi-dimensional response surface, a typical approach is to use Design of Experiments (DoE) for this purpose. For high dimensional problems, strategies such as the Latin Hypercube Sampling [58] or a variety of random number generators (e.g. the Sobol quasi-random sequences of numbers [59]) are used. However, when we do not have freedom to choose the number and the locations of the support points, all existing methods for generating response surfaces have serious problems with accuracy and robustness. This is mainly because arbitrary data sets provide inadequate uniformity of coverage of space of the design variables and clustering of the support points that leads to oscillations of the response surfaces. The most common multidimensional response surface fitting algorithms and their hybrids have been described in [60]: polynomial regression [61], kriging [62], radial basis functions [63], neural networks [64] and self-organizing algorithms [65]. Hybrid methods may also be constructed in order to overcome the shortfalls of single methods. The proposed hybrid methods are the Fittest Polynomial Radial Basis Function (FP-RBF) , Kriging Approximation with Fittest Polynomial Radial Basis Function (KRG-FP-RBF) , Hybrid Self-Organizing Model With RBF and the Genetic Algorithm Based Wavelet Neural Network (HYBWNN) . The article evaluates the performance of these algorithms on data sets containing either a scarce, small, medium or large number of points. 1.5 Hybrid optimization methods The “no free lunch theory” [4] has been introduced in a previous section. This theory says it is impossible to affirm that one single algorithm will outperform all others for all classes of optimization problems. Therefore, the usual way to proceed is the creation of hybrid methods, also called metaoptimization or hyperheuristics. Two are the key aspects to be discussed regarding hybridization: •Constitutive algorithms of a hybrid method may be combined sequentially, in parallel or in a mixed sequential/parallel way [60]. In sequential hybridization, a control algorithm performs automatic switching among the constituent algorithms at each stage during the optimization when the rate of convergence becomes unsatisfactory, the process tends towards a local optimum, or some other undesirable aspect of the iterative process appears [55]. In parallel hybridization, constitutive optimization algorithms run in parallel and contribute a portion of each new generation’s population. The portion that each search contributes to the new generation is dependent on the success of the algorithm to §1.5 Hybrid optimization methods 19 provide past useful solutions to the search. Finally, a sequential/parallel method is a mix of the two other hybridization strategies. Nevertheless, the resulting hybrid method is expected to be more robust and converge faster than its individual constituent optimization algorithms no matter the selected hybridization technique. In the context of this Doctoral Thesis, only sequential hybridization has been considered. •Combining both deterministic and stochastic methods allows achieving a right balance between exploration and exploitation. Deterministic methods are in general computationally faster (they require fewer objective function evaluations) than stochastic methods, although they can converge to a local minimum or maximum, instead of the global one. On the other hand, stochastic algorithms can ideally converge to a global maximum or minimum, although they are computationally slower than the deterministic ones. Indeed, stochastic algorithms can require thousands of evaluations of the objective functions and, in some cases, become non-practical. This is why stochastic methods are usually employed to find the region where the global optimum is located and deterministic methods to get the exact optimal point. The simplest sequential hybrid optimization algorithm is composed by a global (stochastic) search method and a local (deterministic) search method. Regarding the stochastic method, evolutionary algorithms hold a prominent place since their appearance in 1975, being well suited for both continuous and discrete optimization. Genetic algorithms are the most extended method in the evolutionary algorithms’ family. Regarding the deterministic method, it differs depending on the continuous or discrete nature of the optimization problem. Historically, single-objective continuous optimization problems have been thoroughly studied and gradient methods are the most widely used local search methods. However, there is less consensus in the case of continuous multi-objective optimization problems [5]. A survey of deterministic methods for continuous single-objective optimization is carried out in [6]. According to this study, Newton’s method converges more rapidly than the conjugate gradient method , but it has the drawback of long calculation times for the Hessian matrix coefficients. Therefore, it is preferred to approximate the Hessian based only on first order derivatives and avoiding the calculation of second order derivatives. This is done by means of the so called quasi-Newton methods, which have a slower convergence rate than the Newton’s method, but are overall computationally faster. Anyway, quasi-Newton methods converge more rapidly than the conjugate gradient 26 §1.7 Parallelization of genetic and other population-based optimization algorithms times, for example). Hence, the total time of the optimization process is to be reduced by means of parallelization in order to make the use of genetic algorithms viable. Genetic algorithms are naturally prone to parallelism since the operations on the individuals are relatively independent from each other. According to [23], parallelization techniques can be divided into software and hardware parallelization. In the case of genetic algorithms, when both classes of parallelization techniques are applied together an exceptional characteristic arises: the behavior of the parallel algorithm is better than the sum of the separate behaviors of its component sub-algorithms, i.e. the new algorithm is not just the parallel version of a sequential algorithm intended to provide speed gains. Software parallelization is accomplished by using a structured population, either in the form of a set of islands or a diffusion grid, and often leads to superior numerical performance even when the algorithms run on a single processor. Comparing with natural evolution, software parallelization is based on the fact that species form a large population distributed in a certain number of semi-isolated breeding subgroups. The local selection and reproduction rules allow the species evolve locally, and diversity is enhanced by migrations of individuals among the interconnected subgroups. Different search techniques may be utilized in each subgroup. Hardware parallelization is an additional way of speeding up the execution of the algorithm and consists on running a sequential algorithm in several processors. The results obtained in the parallel execution are the same as in the sequential execution, but a certain speedup is achieved in run time. What acceleration could be expected? According to [69], it is possible to have super-linear speedup for certain problems and parameterizations by using hardware parallelization exclusively, both in homogeneous and in heterogeneous parallel machines. The concept of super-linear speedup is understood as the fact that using m processors leads to an algorithm that runs more than m times faster than the sequential version. This assertion is compliant with [70], where Donaldson et al. showed that there is no theoretical upper limit for the speedup in heterogeneous systems. Another important concept related to parallelization is heterogeneity, where two different kinds may be distinguished again. Heterogeneity at software level refers to the use of different search techniques (coding, operators, parameters, etc.) in the population subgroups created by means of software parallelization. If the same search techniques are applied in all subgroups, the algorithm is considered homogeneous at software level. Heterogeneity at hardware level refers to the existence of processors with dif- §1.7 Parallelization of genetic and other population-based optimization algorithms 27 ferent characteristics (architecture, clock rate, etc.) when a genetic algorithm is run on a network of computers or in massively parallel computers. If the characteristics of all processors where the optimization is run are the same, the hardware is considered homogeneous. Finally the concept of scalability is introduced. Scalability measures the ability of a parallel machine and a parallel algorithm to use efficiently a larger number of processors, and depends on both the communication patterns of the algorithm and the infrastructure provided by the machine. Several scalability measures are described in [71]. 1.7.2 Main parallelization strategies According to [23], parallelization strategies of genetic algorithms are commonly divided into 3 types: global parallelization (type 1), coarse grain parallelization (type 2) and fine grain parallelization (type 3). Type 1 corresponds to hardware parallelization. In types 2 and 3, the classification coarse/fine grain relies on the computation/communication ratio. If this ratio is high, the Parallel Genetic Algorithm (PGA) is called a coarse grain algorithm, whereas if low it is called a fine grain PGA. Coarse grain PGAs are the most popular techniques and are also known as distributed or island GAs (dGAs), whereas fine grain PGAs are known as cellular (cGAs), diffusion or massively-parallel GAs. Hybrid algorithms have also been proposed, combining different parallel GAs at two levels in order to enhance the search in some way (see schemes d, e and f in Fig. 1.6). An interesting comparison of various parallel implementations may be found in [72]. Type 1: Global parallelization (master-worker parallelization) Also called explicit parallelization, this strategy implements hardware parallelization and distributes genetic operations and/or objective functions evaluations among several processors (a fraction of the population is assigned to each of the processors). This type presents a viable choice only for problems with a time-consuming function evaluation. Otherwise the communication overhead is higher than the benefits of the parallel execution. According to [73], the global parallel GA is called synchronous if it proceeds in the same way as a sequential GA, i.e. stopping and waiting to receive the fitness values for all the population before starting the next generation. This is the most usual strategy. The global parallel GA is called asynchronous when there is no clear division between 28 §1.7 Parallelization of genetic and other population-based optimization algorithms - 9In a parallel GA there exist many elementary GAs working on separate sub-populations Pi(t). Each sub-algorithm includes an additional phase of periodic communication with a set of neighboring sub-algorithms located on some topology. This communication usually consists in exchanging a set of individuals, although nothing prevents the sub-algorithms of exchanging other kind of information such as population statistics. All the sub-algorithms are thought to perform the same reproductive plan. Otherwise the PGA is heterogeneous [2], [34] (even the representation could differ among the islands, posing new challenges to the exchange of individuals [44]). In a distributed GA demes are loosely-coupled islands of strings (Figure 8b). A cellular GA (Figure 8c) defines a NEWS neighborhood (North-East-West-South in a toroidal grid) in which overlapping demes of 5 strings (4+1) execute the same reproductive plan. In every neighborhood (deme) the new string computed after selection, crossover, and mutation replaces the current one only if it is better (binary tournament), although many other variants are possible. This process is repeated for all the neighborhoods in the grid of the cellular GA (there are as many neighborhoods as strings). With regard to the classes of Figure 5 we suggest in Figure 8 several implementations of PGAs. Global parallelization consists in evaluating, and maybe crossing+mutating in parallel all the structures, while selection uses the whole population. Interesting considerations on global parallelization can be found in [16]. This model provides lower runtime only for slow objective functions; an additional limitation is that the search mechanism uses a single population. The automatic parallelization is rarelly found since the compiler must provide the parallelization of the algorithm automatically. The other hybrid models in Figure 8 combine different parallel GAs at two levels in order to enhance the search in some way. Interesting surveys on these and other parallel models can be found in [1], [16], [17]. Hierarchies of GAs are the most recurrent models found in the literature. In Figure 8d we can appreciate a distributed algorithm in which every island runs a cellular GA. In Figure 8e several algorithms using global parallelization are used to create a ring of islands. Finally, Figure 8f shows two levels of coarse grain PGAs, the inner level having a full-connected topology, and the outer level having a simple ring topology. ... Master Workers Master Workers Workers (e) (f)(a) (b) (c) (d) Figure 8. Different models of PGA: (a) global parallelization, (b) coarse grain, and (c) fine grain. Many hybrids have been defined by combining PGAs at two levels: (d) coarse and fine grain, (e) coarse grain and global parallelization, and (f) coarse grain plus coarse grain. We want to point out that coarse (cgPGA) and fine grain (fgPGA) PGAs are subclasses of the same kind of parallel GA consisting in a set of communicating sub-algorithms. We propose a change in the nomenclature to call them distributed and cellular GAs (dGA and cGA), since the grain is usually intended to refer to their computation/communication ratio, while actual differences can also be found in the way in which they both structure their population (see Figure 9). While a distributed GA has a large sub-population (>>1) a cGA has typically only one string in every sub-algorithm. For a dGA the sub-algorithms are loosely connected, while for a cGA they are tightly connected. In addition, in a dGA there exist only a few sub-algorithms, while in a cGA there is a large number of them. - 9In a parallel GA there exist many elementary GAs working on separate sub-populations Pi(t). Each sub-algorithm includes an additional phase of periodic communication with a set of neighboring sub-algorithms located on some topology. This communication usually consists in exchanging a set of individuals, although nothing prevents the sub-algorithms of exchanging other kind of information such as population statistics. All the sub-algorithms are thought to perform the same reproductive plan. Otherwise the PGA is heterogeneous [2], [34] (even the representation could differ among the islands, posing new challenges to the exchange of individuals [44]). In a distributed GA demes are loosely-coupled islands of strings (Figure 8b). A cellular GA (Figure 8c) defines a NEWS neighborhood (North-East-West-South in a toroidal grid) in which overlapping demes of 5 strings (4+1) execute the same reproductive plan. In every neighborhood (deme) the new string computed after selection, crossover, and mutation replaces the current one only if it is better (binary tournament), although many other variants are possible. This process is repeated for all the neighborhoods in the grid of the cellular GA (there are as many neighborhoods as strings). With regard to the classes of Figure 5 we suggest in Figure 8 several implementations of PGAs. Global parallelization consists in evaluating, and maybe crossing+mutating in parallel all the structures, while selection uses the whole population. Interesting considerations on global parallelization can be found in [16]. This model provides lower runtime only for slow objective functions; an additional limitation is that the search mechanism uses a single population. The automatic parallelization is rarelly found since the compiler must provide the parallelization of the algorithm automatically. The other hybrid models in Figure 8 combine different parallel GAs at two levels in order to enhance the search in some way. Interesting surveys on these and other parallel models can be found in [1], [16], [17]. Hierarchies of GAs are the most recurrent models found in the literature. In Figure 8d we can appreciate a distributed algorithm in which every island runs a cellular GA. In Figure 8e several algorithms using global parallelization are used to create a ring of islands. Finally, Figure 8f shows two levels of coarse grain PGAs, the inner level having a full-connected topology, and the outer level having a simple ring topology. ... Master Slaves Master Slaves Workers (a) (b) (c) (d) (e) (f) Figure 8. Different models of PGA: (a) global parallelization, (b) coarse grain, and (c) fine grain. Many hybrids have been defined by combining PGAs at two levels: (d) coarse and fine grain, (e) coarse grain and global parallelization, and (f) coarse grain plus coarse grain. We want to point out that coarse (cgPGA) and fine grain (fgPGA) PGAs are subclasses of the same kind of parallel GA consisting in a set of communicating sub-algorithms. We propose a change in the nomenclature to call them distributed and cellular GAs (dGA and cGA), since the grain is usually intended to refer to their computation/communication ratio, while actual differences can also be found in the way in which they both structure their population (see Figure 9). While a distributed GA has a large sub-population (>>1) a cGA has typically only one string in every sub-algorithm. For a dGA the sub-algorithms are loosely connected, while for a cGA they are tightly connected. In addition, in a dGA there exist only a few sub-algorithms, while in a cGA there is a large number of them. Figure 1.6: Different models of PGA (extracted from [23]): (a) global parallelization, (b) coarse grain parallelization, (c) fine grain parallelization, (d) coarse + fine grain hybrid, (e) coarse grain + global hybrid, and (f) coarse grain + coarse grain hybrid. generations, i.e. if any worker processor that finishes evaluating an individual returns it to the master and receives another individual. A high level of processor utilization is achieved with this strategy, despite the workers’ heterogeneous processor speeds. However, the obtained search results differ from those achieved by the sequential GA. The global parallelization uses a single population. The main advantage is the simplicity of implementation, because parallelization takes place only at the level of objective function calculation. The disadvantage is that no software parallelization is utilized. A schematic representation is provided in Fig. 1.6 (a). The method can be implemented efficiently on sharedand distributed-memory computers [73]. On a shared-memory multiprocessor, the population can be stored in shared memory and each processor could read a fraction of the population and write back the evaluation results without any conflicts. On a distributed-memory computer, the population is stored in one processor. This “master” processor is the responsible for sending the individuals to the other processors (the “workers”) for evaluation, collecting the results, and applying the genetic operators to produce the next generation. The difference with a shared-memory implementation is that the master has to send and receive messages explicitly. §1.7 Parallelization of genetic and other population-based optimization algorithms 29 Type 2: Coarse grain parallelization (island GA, distributed GA) As it was mentioned in a previous section, this strategy is based on distributing the whole population in a certain number of semi-isolated breeding subgroups called islands. An independent optimization takes place in each island, where a new population of individuals is created from the old one by applying genetic operators such as selection, crossover, mutation and replacement. From time to time some individuals are exchanged among the interconnected subgroups (see Fig. 1.6 (b)). Migration is the operator that guides this exchange of individuals among islands. The main concepts related to migration, which are the object of study of most publications on the island GA field [74], are mentioned hereafter. The migration gap is the number of steps (generations) in every sub-population between two successive exchanges. Migration may take place in every island periodically or by using a given probability of migration to decide in every step whether migration will take place or not. Island GAs are usually synchronous, i.e. the phases of migration and reception of external migrants are embedded in the same portion of the algorithm (migration takes place when all islands achieve a fixed number of generations). However, synchronous migration has the drawback of being slow for some problems [23]. In asynchronous island GAs, migrants are sent whenever it is needed and migrants are accepted whenever they arrive. This behavior is accomplished by implementing a chromosome buffer [74]. The migration rate is the parameter determining the number of individuals that undergo migration in every exchange. It is not clear which is the best value for the migration rate, but best results have been obtained for low percentages (ranging from 1% to 10% of the population size). Selection and replacement operators in the migration procedure of each island are commonly the same as the ones used in the independent optimization process taking place in that island. The topology is the map of interconnections between islands, being the general tendency to use static topologies that are set before launching the algorithm and remain unchanged [74]. It seems that the ring and hyper-cube are two of the best topologies for most problems, whereas full-connected and centralized topologies have shown problems in their parallelization and scalability due to the tight connectivity. Another approach is the use of dynamic topologies, giving rise to decentralization. In this kind of topology, all populations can be connected to all other populations directly. The main disadvantage of this approach is the complexity of implementation and, moreover, the efficiency and 30 §1.7 Parallelization of genetic and other population-based optimization algorithms usefulness of dynamic topologies has not been proved [74]. The performance of distributed GAs is often better than for the sequential GAs. The two main reasons are that subpopulations are run simultaneously using several processors (reducing the whole processing time) and that this kind of search maintains samples of very different promising zones of the search space, increasing the efficacy of the algorithm. The coarse-grain parallelization has several advantages with respect to fine-grain parallelization. On one hand, it is possible to adapt existing sequential algorithms for being used in each subpopulation (island). On the other hand, the number of available processors does not affect the result of the optimization. Thus, the strategy is suitable in the case of having limited resources. Moreover, it could be adapted more easily than the cellular GA to grid computing [75], where the number and performance of processors may vary during the optimization. The main disadvantage of a multiple-population GA might be that the critical path of a fine-grained algorithm seems to be shorter [76]. Type 3: Fine grain parallelization (cellular GA, diffusion GA, massively parallel GA) A cellular GA is implemented at large computer terminals [74] and consists of one spatially distributed population in which overlapping subpopulations execute the same reproductive plan. In every neighborhood the new individual computed after selection (usually by means of a binary tournament), crossover and mutation replaces the current one only if it is better. This process is repeated for all the neighborhoods in the grid of the cellular GA, where there are as many neighborhoods as individuals. The most popular algorithm uses a toroidal grid of individuals (toroidal topology) and defines a NEWS neighborhood (North-East-West-South) in which subpopulations of 5 individuals (4+1) execute the reproductive plan. Thus, “good” individuals are spread over the whole distributed population by means of migrations happening between neighboring subpopulations (4 subpopulations available per individual if a NEWS neighborhood is used). A schematic representation is provided in Fig. 1.6 (c). This strategy is very well suited for massively parallel computers, because it divides the population into a large number of parts. One individual is usually evaluated in each processor. Thus, the cellular GA is a completely decentralized model in which communications are equally distributed and where each cell has to wait only to a few other cells. Moreover, it is proven in [76] that the critical path of a fine-grained algorithm is shorter than that of a multiple-population GA. This means that if enough §1.7 Parallelization of genetic and other population-based optimization algorithms 31 processors were available, massively parallel GAs would need less time to finish their execution, regardless of the population size. Note, however, that considerations such as the communications bandwidth or memory requirements were not included in the mentioned theoretical study. The main disadvantage of this strategy is that the use of a small number of processors results in degeneration of the whole population, leading the genetic algorithm into a local optimum [74]. Therefore, this model should not be utilized unless a minimum number of processors is available. Another drawback may be that sequential algorithms cannot be adapted to fine-grain parallelization, being necessary the development of completely new algorithms. Hybrid parallelization models When two or more GA parallelization methods are combined they form a hierarchy [76], i.e. a hybrid parallelization strategy, and a better performance than with any of the constituent methods alone is achieved. Several distinct combinations may be proposed, as it is shown in Fig. 1.6 (d), (e) and (f). Note that these three models use a coarse-grained algorithm at the upper level, which is usual in hierarchical GAs [73]. A common hybrid strategy consists on combining a coarse-grained GA at the upper level and a global GA (master-worker model) at each of the subpopulations [76]. Thus, independent optimizations using the master-worker strategy are carried out in each island and some individuals migrate between islands from time to time. This approach is useful when objective functions that need a considerable amount of computation time are evaluated. 1.7.3 Shortcomings of standard parallelization techniques Description of shortcomings The integration of optimization methodologies (and particularly of genetic algorithms) with computational analyses/simulations has a profound impact on product design. Nowadays the main challenges related to the application of such methodologies arise from high-dimensionality of problems, computationally-expensive analysis/simulations and unknown function properties (i.e. black-box functions). An extensive review on these topics is available in [77]. The consequence is that the computational load of the optimization process may be considerable. This fact makes the development of appropriate parallelization techniques a critical issue in order to benefit from optimization strategies in real-world design problems. 32 §1.7 Parallelization of genetic and other population-based optimization algorithms The goal of every parallelization strategy is to maximize CPU usage by minimizing the communication overhead and by balancing the computational load correctly, thus avoiding idleness of allocated processors. In the case of traditional optimization algorithms (and particularly of genetic algorithms), the principal causes of computational load imbalance are listed below. •Inappropriate ratio (no. individuals / no. processor groups) Traditionally, when several individuals may be evaluated simultaneously, each individual is assigned to a group of processors. In standard algorithms, all groups are formed by the same number of processors and remain unchanged during optimization. Thus, the number of individuals assigned to a certain group of processors is calculated as the total number of individuals in the population divided by the number of processor groups. In case the residual of this quotient is not zero, some processor groups will evaluate one more individual than the other groups. Hence, the processor groups with fewer individuals to evaluate will remain idle while the rest of the groups finish their pending evaluations. For this explanation, all individuals are assumed to have homogeneous evaluation times and the number of allocated processor groups to be equal or smaller than the number of individuals in the population. Such a problem is reported for example in [78], where the optimum shape design of aerodynamic configurations is studied. The genetic algorithm evaluates an objective function that calls an unstructured grid-based CFD solver. The objective evaluation time represents between 80% and 90% of the total elapsed time. A two-level parallelization strategy is utilized, having a master-worker model at the upper level and a parallel evaluation of the objective function at the lower level. Two equally sized processor groups are responsible of concurrently evaluating two individuals. The evaluation of an odd number of airfoils may cause computational load imbalance, influencing negatively the parallel performance of the algorithm. •Heterogeneous parallel computer systems A computer system is a collection of n -processors interconnected by a communication network, specified by the pairwise latency and the bandwidth between processors. Ideally, computer systems would be homogeneous, i.e. all nodes would have identical architecture, clock rate and latency. Unfortunately, a certain degree of heterogeneity is hardly avoidable even between nodes sharing the same architecture, because additional factors such as the hardware’s age, refrigeration, etc. affect their performance. This hardware heterogeneity translates into §1.7 Parallelization of genetic and other population-based optimization algorithms 33 non-homogeneous objective evaluation times and the consequent loss of parallel performance in the case of optimization algorithms. In this regard, a general model to define, measure and predict the efficiency of applications running on heterogeneous parallel computer systems is presented in [79]. The effect of the imbalance caused by heterogeneous parallel systems is noticeable, for instance, if the parallel genetic algorithm implemented by Cantú-Paz [80] is used. This algorithm was designed for a homogeneous parallel computer system and distributes the same number of individuals per processor. Therefore, a considerable loss of parallel performance can be expected in heterogeneous environments. The most extreme situation regarding hardware heterogeneity occurs probably in grid computing. Grid computing consists of a geographically distributed infrastructure gathering computer resources around the world, and has emerged as an effective environment for the execution of parallel applications that require great computing power [75]. Grid computing environments provide an attractive infrastructure for implementing parallel metaheuristics. However, the fact that grid resources are distributed, heterogeneous and non-dedicated makes writing parallel grid-aware applications very challenging, because one has to address the issues of grid resource discovery and selection, grid job preparation, submission, monitoring and termination which differ from one middleware to another. •Heterogeneous objective function evaluation time Heterogeneity in the computational cost of evaluating the objective functions may cause an important load imbalance, provided that the simulation time of the genetic algorithm is dominated by the evaluation time of the objective functions. This phenomenon happens when the evaluation cost of individuals is dependent on the optimization variables, i.e. the input variables of the objective functions, and is not unusual in heat transfer and nonlinear mechanics applications. This kind of imbalance is reported for example in [78], a publication already mentioned in this document. As it was explained, two individuals are concurrently evaluated in two processor groups. Computational load imbalance is observed and caused by the fact that the required number of iterations for evaluating the objective function depends on the shape of the airfoil, i.e. depends on its optimization variables. 34 §1.7 Parallelization of genetic and other population-based optimization algorithms Current research Ongoing research regarding the previously introduced shortcomings is summarized hereafter. The main causes of computational load imbalance are listed again, together with some currently available solutions. •Inappropriate ratio (no. individuals / no. processor groups) The standard and simplest solution is to adjust the ratio so that the residual of the quotient becomes zero. This may be achieved by either modifying the number of processor groups or resizing the population managed by the genetic algorithm. Nevertheless, the coupling between these two terms is annoying and might involve some drawback when configuring an optimization algorithm to be used in massively parallel computers. As an example, a population with an even number of airfoils was created in [78] once the load imbalance was detected due to the use of an odd number of individuals. •Heterogeneous parallel computer systems Several studies have been found which try to tackle this source of imbalance. Some of them are cited hereafter, as a sample of the solutions proposed up to date. Genetic algorithms are discussed from an architectural perspective in [81], offering a general analysis of performance of GAs on multi-core CPUs and on many-core GPUs. As a conclusion, the authors propose the best parallel GA scheme for multi-core, multi-socket multi-core and many-core architectures. The use of information of processors’ heterogeneity to make the distribution of individuals nonhomogeneous is proposed in [79]. However, the approach is said to be highly undesirable because it requires that the program interacts with the resource management software, which contains the speeds of the processors allocated to the job. The implementation of such a strategy results in a non-portable program, due to the lack of standard interfaces for supplying this information. In the context of grid computing, some tools for managing hardware heterogeneity have been recently developed. WoBinGO [75], for instance, is a framework for solving optimization problems over heterogeneous resources, including HPC clusters and Globus-based grids. It uses a master-worker parallelization model, which is easily replaceable by a hierarchical parallel GA with master-worker demes or by a parallel cellular GA. §1.7 Parallelization of genetic and other population-based optimization algorithms 35 •Heterogeneous objective function evaluation time A very basic (and not advisable) solution to homogenize all objective function evaluation times is to limit their maximum allowed duration, e.g. by establishing a maximum number of iterations to their solver (this method is used for instance in [78]). Nonetheless, it must be borne in mind that this alternative affects the obtained results and is not suitable for general use. •Solutions to any kind of heterogeneity (hardwareor software-based) Although the source of hardwareand software-based heterogeneity is different, the consequence is common: both affect the time needed to evaluate each individual. This is the reason why some proposed solutions are valid for these two types of heterogeneity and have been summarized in this section. A master-worker model with a constant and homogeneous number of individuals assigned to each processor is described in [73]. However, the need of balancing the computational load among processors using a dynamic scheduling algorithm like guided self-scheduling is mentioned. The Adaptive Parallel Genetic Algorithm is proposed in [79]. This algorithm automatically changes the number of individuals to be evaluated on each processor depending on the evaluation time of each individual. The algorithm uses a server–client blocking message architecture, in which the server node is responsible for the distribution of work and the evolution of the genetic algorithm. The client nodes evaluate the fitness function for the individuals and return their values to the server node (see Fig. 1.7). Communication between processors is necessary for distributing and/or balancing the evaluation of a population over the nodes. Using this scheme a fast processor will request more work than a slow one and as a result, the algorithm will send more individuals to the faster processors adapting the algorithm to the heterogeneity of the system. The same will happen in the case of processors that receive individuals with short evaluation times, adapting the algorithm to heterogeneous objective function evaluation times. A similar alternative is described in [82], where an asynchronously global parallel genetic algorithm with 3-tournament elimination selection is introduced. The difference between the traditional master-worker GA and this algorithm is in the tasks performed by the master and the workers. In the traditional algorithm workers only evaluate individuals, whereas the master distributes individuals among workers as well as performing all genetic operations. In the new algorithm, the master only initializes the population, whereas workers perform the whole 42 §1.8 Test suites for optimization algorithms Schwefel Minimize f(〈x1,..., xn〉)= n X i=1−xisin³p|xi|´+418.9829n xi∈[−512.03,511.97] (1.4) where the minimum f(x)≈0 is located at x=(420.9687,420.9687,...). -2 -1 0 1 2 -2 -1 0 1 2 0 1000 2000 - 2 -1 0 1 2 0 1000 2000 (a) -5 -2.5 0 2.5 5 -5 -2.5 0 2.5 5 0 500 1000 1500 -5 -2.5 0 2.5 5 0 500 1000 1500 -400 -200 0 200 400 -400 -200 0 200 400 -500 0 500 - 400 -200 0 200 400 -500 0 500 -500 -250 0 250 500 -500 -250 0 250 500 0 50 100 150 - 500 -250 0 250 500 0 50 100 150 -40 -20 0 20 40 -40 -20 0 20 40 0 1 2 3 -40 -20 0 20 40 0 1 2 3 (b) (c) (d) (e) Figure 1.8: Test functions for single-objective optimization (extracted from [90]): (a) Rosenbrock, (b) Rastrigin, (c) Schwefel, (d) Griewank, and (e) Griewank (detail). §1.8 Test suites for optimization algorithms 43 1.8.2 Multi-objective optimization tests In the case of multi-objective optimization, typical test suites are gathered in [5] under the names of MOP, ZDT, DTLZ, OKA, WFG, etc. Two typical tests for unconstrained optimization have been selected, namely MOP4 (see Fig. 1.9) and ZDT6 (see Fig. 1.10): MOP4 (Kursawe’s function) Minimize the following objective functions f1(x)= 2 X i=1³−10exp ³−0.2qx2 i+x2 i+1´´ (1.5) f2(x)= 3 X i=1¡|xi|0.8+5sin(xi)3¢(1.6) where −5≤xi≤5i=1,...,3. ZDT6 (Zitzler-Deb-Thiele’s function N. 6) Minimize the following objective functions f1(x)=1−exp(−4x1)sin6(6πx1) (1.7) f2(x,g)=g(x)·µ1−µf1(x) g(x)¶2¶(1.8) where 0 ≤xi≤1i=1,...,10 and g(x)=1+9ÃP10 i=2xi 9!0.25 (1.9) 44 §1.8 Test suites for optimization algorithms -12 -8 -4 0 -20 -18 -16 -14 f2(x) f1(x) Figure 1.9: True Pareto Front of Kursawe’s function (MOP4). 0 0.5 1 1.5 0.5 0.75 1 f2(x) f1(x) Figure 1.10: True Pareto Front of Zitzler-Deb-Thiele’s function N. 6 (ZDT6). §1.10 Conclusions 45 1.9 Random number generators It has been said that evolutionary algorithms (among others) are stochastic methods, which means that they employ randomness to some degree. This implies that the quality of the produced results also relies on the quality of the selected random number generator. Consequently, it is necessary to carry out a brief state of the art research related to this topic. A pseudorandom number generator is an algorithm for generating a sequence of numbers whose properties approximate the properties of sequences of random numbers. However, the generated sequence is not truly random, being completely determined by a relatively small set of initial values. These values are called the generator’s seed. Distributions of most usual programming languages like C, C++ and Java include pseudorandom number generators. Nevertheless, their quality has been strongly questioned, as in [90], making the search for other alternatives highly advisable. The renowned mathematical library Intel ® Math Kernel Library (MKL) [91] provides various pseudorandom, quasi-random, and non-deterministic random number generators: Wichmann-Hill pseudorandom number generator [92], Mersenne Twister MT19937 pseudorandom number generator [93], SIMD-oriented Fast Mersenne Twister SFMT19937 pseudorandom number generator [94], Sobol quasi-random number generator [59], Niederreiter quasi-random number generator [95], etc. Studies of advantages and deficiencies of random number generators are also available, for instance in [96]. Regarding the consulted bibliography on evolutionary algorithms, references to the choice of an appropriate random number generator have been found in [60], where Sobol’s pseudorandom sequence generator [59] is employed, and in [97], where the Mersenne Twister MT19937 [93] is used. All these considerations have been taken into account in Chapter 2, where a random number generator is selected for implementation. 1.10 Conclusions The state of the art of optimization algorithms has been studied in this chapter. After a brief introduction to the goals of optimization, the main concepts have been defined and a classification of currently available search techniques for both single-objective and multi-objective optimization has been presented. All those algorithms must be used together with efficient constraint handling methods in order to guarantee the overall efficiency of the search. 46 References The evaluation of objective functions can be computationally expensive in real-life problems, making the cost of the whole optimization process prohibitive. Surrogatebased optimization has been introduced as an effective method of getting a lower fidelity model of the original function which can be simulated in a reasonable time, making possible the application of mathematical optimization techniques to such problems. Another difficulty faced by the optimization community is the lack of a single search technique which always outperforms all other available techniques. Thus, the best practical solution consists on creating hybrid optimization methods, which are able to select the best performing constituent algorithm at each moment. Evolutionary algorithms are the most widely used global optimization methods to date, among which genetic algorithms hold a prominent place. The main characteristics of such algorithms have been defined and an extensive review of the current state of the art regarding their parallelization has been carried out. Common test suites for comparing the performance of optimization algorithms have also been introduced and a selection of some tests has been made. Finally, the importance of choosing a good random number generator has been remarked in the case of using stochastic optimization algorithms. The most widespread generators have been mentioned, together with some software packages that include their implementation. Acknowledgments This work has been financially supported by a FPU Doctoral Grant (Formación de Profesorado Universitario) awarded by the Ministerio de Educación, Cultura y Deporte, Spain (FPU12/06265) and by Termo Fluids S.L. The author thankfully acknowledges these institutions. References [1] C. S. Beightler, D. T. Phillips, and D. J. Wilde. Foundations of optimization. Prentice Hall Inc., 1979. [2] R. L. Haupt and S. E. Haupt. Practical Genetic Algorithms. John Wiley & Sons, 2004. [3] D. E. Goldberg. Genetic algorithms in search, optimization, and machine learning. Addison-Wesley, 1989. References 47 [4] D. H. Wolpert and W. G. Macready. No free lunch theorems for optimization. IEEE Transactions on Evolutionary Computation, 1(1):67–82, 1997. [5] C. A. Coello Coello, G. B. Lamont, and D. A. Van Veldhuizen. Evolutionary algorithms for solving multi-objective problems. Springer Science & Business Media, 2007. [6] M. J. Colaço and G. S. Dulikravich. A survey of basic deterministic, heuristic and hybrid methods for single-objective optimization and response surface generation. In METTI IV - Thermal Measurements and Inverse Techniques, volume 1, Rio de Janeiro, Brazil, 2009. [7] G. B. Dantzig and R. Cottle. The basic George B. Dantzig. Stanford University Press, 2003. [8] J. A. Nelder and R. Mead. A simplex method for function minimization. The Computer Journal, 7(4):308–313, 1965. [9] M. J. Box. A comparison of several current optimization methods, and the use of transformations in constrained problems. The Computer Journal, 9(1):67–77, 1966. [10] H. P. Schwefel. Evolution and optimum seeking. Sixth-generation computer technology series, 1995. [11] M. J. D. Powell. An efficient method for finding the minimum of a function of several variables without calculating derivatives. The Computer Journal, 7(2):155–162, 1964. [12] C. G. Broyden. A class of methods for solving nonlinear simultaneous equations. Mathematics of Computation, 19(92):577–593, 1965. [13] R. Fletcher. Generalized inverses for nonlinear equations and optimization, Numerical Methods for Non-linear Algebraic Equations. Gordon & Breach (R. Rabinovitz, Ed.), London, UK, 1963. [14] D. Goldfarb and L. Lapidus. Conjugate gradient method for nonlinear programming problems with linear constraints. Industrial & Engineering Chemistry Fundamentals, 7(1):142–151, 1968. [15] D. F. Shanno. An accelerated gradient projection method for linearly constrained nonlinear estimation. SIAM Journal on Applied Mathematics, 18(2):322–334, 1970. 48 References [16] D. G. Luenberger and Y. Ye. Linear and nonlinear programming, volume 2. Springer, 1984. [17] J. L. Zhou and A. Tits. User’s guide for FFSQP version 3.7: A Fortran code for solving optimization programs, possibly minimax, with general inequality constraints and linear equality constraints, generating feasible iterates. Technical Report SRC-TR-92-107r5, Institute for Systems Research, University of Maryland, 1997. [18] J. Pearl. Heuristics: Intelligent search strategies for computer problem solving. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1984. [19] F. Glover and M. Laguna. Tabu search. Kluwer Academic Publishers, 1997. [20] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. Science, 220(4598):671–680, 1983. [21] K. E. Parsopoulos and M. N. Vrahatis. Recent approaches to global optimization problems through particle swarm optimization. Natural Computing, 1(2-3):235– 306, 2002. [22] M. Dorigo and L. M. Gambardella. Ant colony system: a cooperative learning approach to the traveling salesman problem. IEEE Transactions on Evolutionary Computation, 1(1):53–66, 1997. [23] E. Alba and J. M. Troya. A survey of parallel distributed genetic algorithms. Complexity, 4(4):31–52, 1999. [24] J. H. Holland. Adaptation in natural and artificial systems. Ann Arbor: University of Michigan Press, 1975. [25] T. Bäck, D. B. Fogel, and Z. Michalewicz. Handbook of evolutionary computation. IOP Publishing Ltd., Bristol, UK, 1st edition, 1997. [26] R. Storn and K. Price. Differential evolution: a simple and efficient adaptive scheme for global optimization over continuous spaces. Technical Report TR-95012, International Computer Science Institute, Berkeley, California, 1995. [27] L. J. Fogel. Artificial intelligence through simulated evolution. John Wiley, New York, 1966. References 49 [28] J. R. Koza. Genetic programming: on the programming of computers by means of natural selection, volume 1. MIT press, 1992. [29] R. G. Reynolds. An introduction to cultural algorithms. In A. V. Sebald and L. J. Fogel, editors, Evolutionary Programming - Proceedings of the Third Annual Conference, pages 131–139, San Diego, CA, USA, 24-26 February 1994. World Scientific Press. [30] A. Osyczka. Multicriteria optimization for engineering design. Design Optimization, 1:193–227, 1985. [31] D. A. Van Veldhuizen. Multiobjective evolutionary algorithms: classifications, analyses, and new innovations. PhD thesis, Wright Patterson AFB, OH, USA, 1999. AAI9928483. [32] A. Liefooghe, L. Jourdan, and E. G. Talbi. A unified model for evolutionary multiobjective optimization and its implementation in a general purpose software framework: ParadisEO-MOEO. Research Report RR-6906, INRIA, 2009. [33] J. D. Schaffer. Multiple objective optimization with vector evaluated genetic algorithms. In Proceedings of the 1st International Conference on Genetic Algorithms, pages 93–100. L. Erlbaum Associates Inc., 1985. [34] E. Zitzler and S. Künzli. Indicator-based selection in multiobjective search. In International Conference on Parallel Problem Solving from Nature (PPSN VIII), volume 3242, pages 832–842. Springer-Verlag, 2004. [35] N. Beume, B. Naujoks, and M. Emmerich. SMS-EMOA: Multiobjective selection based on dominated hypervolume. European Journal of Operational Research, 181(3):1653–1669, 2007. [36] M. Ehrgott. Multicriteria optimization. Springer Science & Business Media, 2005. [37] K. Deb, M. Mohan, and S. Mishra. Evaluating the ε -domination based multiobjective evolutionary algorithm for a quick computation of Pareto-optimal solutions. Evolutionary Computation, 13(4):501–525, 2005. [38] J. Molina, L. V. Santana, A. G. Hernández-Díaz, C. A. Coello Coello, and R. Caballero. g-dominance: Reference point based dominance for multiobjective metaheuristics. European Journal of Operational Research, 197(2):685–692, 2009. 50 References [39] S. Russell and P. Norvig. Artificial intelligence: A modern approach. Prentice-Hall, Upper Saddle River, New Jersey, 1995. [40] C. M. Fonseca and P. J. Fleming. Genetic algorithms for multiobjective optimization: formulation, discussion and generalization. In Proceedings of the Fifth International Conference on Genetic Algorithms, volume 93, pages 416–423, University of Illinois at Urbana-Champaign, San Mateo, California, 1993. Morgan Kaufmann Publishers. [41] H. Horn, N. Nafpliotis, and D. E. Goldberg. A niched Pareto genetic algorithm for multiobjective optimization. In Proceedings of the First IEEE Conference on Evolutionary Computation, IEEE World Congress on Computational Intelligence, volume 1, pages 82–87, Piscataway, New Jersey, 1994. [42] N. Srinivas and K. Deb. Multiobjective optimization using nondominated sorting in genetic algorithms. Evolutionary Computation, 2(3):221–248, 1994. [43] J. D. Knowles and D. W. Corne. Approximating the nondominated front using the Pareto archived evolution strategy. Evolutionary Computation, 8(2):149–172, 2000. [44] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2):182–197, 2002. [45] E. Zitzler and L. Thiele. Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach. IEEE Transactions on Evolutionary Computation, 3(4):257–271, 1999. [46] E. Zitzler, M. Laumanns, and L. Thiele. SPEA2: Improving the strength Pareto evolutionary algorithm. In Eurogen, volume 3242, pages 95–100, Athens, Greece, 2001. [47] D. W. Corne, J. D. Knowles, and M. J. Oates. The Pareto envelope-based selection algorithm for multiobjective optimization. In International Conference on Parallel Problem Solving from Nature (PPSN VI), volume 1917, pages 839–848. SpringerVerlag, 2000. [48] K. Deb, M. Mohan, and S. Mishra. Towards a quick computation of well-spread Pareto-optimal solutions. In Second International Conference on Evolutionary Multi-Criterion Optimization, volume 2632, pages 222–236, Faro, Portugal, 2003. Springer. References 51 [49] A. Liefooghe, M. Basseur, L. Jourdan, and E. G. Talbi. ParadisEO-MOEO: A framework for evolutionary multi-objective optimization. In International Conference on Evolutionary Multi-Criterion Optimization (EMO 2007), volume 4403, pages 386–400, Matsushima, Japan, 2007. Springer. [50] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and V. Grunert da Fonseca. Performance assessment of multiobjective optimizers: an analysis and review. IEEE Transactions on Evolutionary Computation, 7(2):117–132, 2003. [51] K. Deb. Multi-objective optimization using evolutionary algorithms, volume 16. John Wiley & Sons, Chichester, UK, 2001. [52] M. Basseur, F. Seynhaeve, and E. G. Talbi. Design of multi-objective evolutionary algorithms: Application to the flow-shop scheduling problem. In Proceedings of the IEEE Congress on Evolutionary Computation (CEC 2002), volume 2, pages 1151–1156, Piscataway, NJ, USA, 2002. [53] H. Meunier, E. G. Talbi, and P. Reininger. A multiobjective genetic algorithm for radio network optimization. In Proceedings of the IEEE Congress on Evolutionary Computation (CEC 2000), volume 1, pages 317–324, San Diego, USA, 2000. [54] C. Grosan, M. Oltean, and D. Dumitrescu. Performance metrics for multiobjective optimization evolutionary algorithms. In Proceedings of Conference on Applied and Industrial Mathematics (CAIM), Oradea, 2003. [55] G. S. Dulikravich, T. J. Martin, M. J. Colaço, and E. J. Inclan. Automatic switching algorithms in hybrid single-objective optimization. FME Transactions, 41(3):167– 179, 2013. [56] IOSO NM Version 1.0, User’s Guide. IOSO Technology Center, Moscow, Russia, 2003. [57] R. Jin, W. Chen, and T. W. Simpson. Comparative studies of metamodeling techniques under multiple modeling criteria. In Proceedings of the 8th AIAA / USAF / NASA / ISSMO Multidisciplinary Analysis & Optimization Symposium, Long Beach, CA, 2000. [58] M. D. McKay, W. Conover, and R. J. Beckman. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21(2):239–245, 1979. 58 §2.1 Introduction 2.1 Introduction 2.1.1 General design requirements of a library The aim of the work presented in this chapter is the creation of a new optimization framework, i.e. a set of classes that embody an abstract design for solutions to a family of related problems [1]. However, writing a new library is always a tough and time consuming task. Therefore, it is necessary to take a while for thinking of a good conceptual design and its appropriate implementation. This allows obtaining a wellstructured and extendable code, easing future maintenance and making the use of the library more comfortable for other researchers. A good review of key design aspects may be found in [1]. The main features expected from a library are summarized hereafter: •Code reusability: Reusability may be defined as the ability of software components to build many different applications. Code development is time consuming and error-prone, what makes development of new code from scratch every time a new problem arises highly undesirable. On the contrary, the use of thoroughly tested and well-documented libraries can save considerable time and headaches. It must be said that the object-oriented paradigm is particularly well-suited to develop reusable, flexible and extendable libraries by means of a hierarchical set of classes. •Conceptual separation between the solution method and the problem: The solution method must be abstract enough to be able to solve similar but distinct problems. •Code correctness and reliability: The developed code must be free of bugs, provide detailed error messages, avoid deadlocks at run time, have well tested algorithms, etc. •An adequate programming language: It is essential to choose the language that best suits the requirements of parallel and distributed programming, other software expected to be coupled to the framework, performance, etc. •Portability: The framework must be deployable on platforms with variable architectures (networks of PCs and workstations, massively parallel machines, etc.) and their associated operating systems. Therefore, it is important to code using a portable language and standard libraries. §2.1 Introduction 59 •Performance and efficient parallelization: The good performance of the code is a must due to the high computational cost of simulations. Moreover, an efficient and easy-to-handle parallelization is desirable in order to take advantage of High Performance Computing (HPC) resources. •Ease-of-use: The user-friendliness of the framework must guarantee access to full functionality with a minimum effort. This involves the implementation of a graphical user interface, simulation monitoring, documentation, etc. 2.1.2 Specific design requirements of the new optimization library Apart from the already mentioned general design requirements, each new development involves additional specific characteristics. In agreement with the objectives of this Doctoral Thesis, these specific features of the new library are listed hereafter: • The aim of the Doctoral Thesis is the development of a generic mathematical optimization tool, applicable in any field of science and engineering. Nevertheless, being this research activity hosted by the Heat and Mass Transfer Technological Center (CTTC), a special focus is to be put on the application of the library to the fields of expertise of the Center: Computational Fluid Dynamics and Heat Transfer (CFD & HT), multi-physics simulation, etc. • CTTC is developing its own software called TermoFluids [2] since several years ago. It is a must for the new library to be compatible with the software packages belonging to TermoFluids and to provide an appropriate coupling interface. The interaction between libraries may be motivated either because the new optimization framework needs to access basic in-house libraries or because an optimization problem is defined using in-house CFD & HT libraries. • CFD & HT is a field of engineering characterized by the need of solving huge nonlinear equations systems involving a high computational cost. Consequently, the use of High Performance Computing (HPC) infrastructures is a common practice. The new optimization library must be designed for maximizing its performance in such equipment, being portability a key aspect to ensure the usability of the code in worldwide supercomputers. Moreover, the irruption of the massively parallel computing paradigm is to be taken into account in the framework’s design. • The new library must be subject to the coding standards at CTTC. This involves 60 §2.1 Introduction using C++ as object-oriented programming language and MPI (Message Passing Interface) as the communication standard for parallel computing. • The development of a new library is always a tough task. The use of third-party libraries is allowed in case it is considered useful in terms of shortening the project’s duration or enhancing the quality of the final implementation. However, the following restrictions are to be taken into account: i) only open-source libraries protected by a GNU LGPL or a less restrictive license are accepted, i.e. CTTC must be free to independently decide the most appropriate license for the new optimization framework, ii) only libraries written in C or C++ are accepted, iii) the new library must hold the core algorithms in order to allow changing any feature of the behavior of the optimizer on demand. 2.1.3 Computing facilities at CTTC Although portability has already been highlighted as a key feature of the new code, its development and testing have been carried out using the supercomputing facilities available at CTTC. The information of both available computer clusters is attached hereafter (see also Fig. 2.1): •Cluster JFF2 (dating from 2009): This Beowulf HPC cluster called Joan Francesc Fernàndez 2nd Generation (JFF2) has 128 cluster nodes, each node counting 2 AMD Opteron Quad Core processors with 16 Gigabytes of RAM memory. The nodes are linked with an infiniband DDR 4X network interconnection with latencies of 2.6 microseconds and a 20 Gbits/s bandwidth. •Cluster JFF3 (dating from 2011): This Beowulf HPC cluster JFF third generation has 40 cluster nodes, each node counting 2 AMD Opteron with 16 Cores for each CPU linked with 64 Gigabytes of RAM memory and an infiniband QDR 4X network interconnection between nodes with latencies of 1.07 microseconds and a 40 Gbits/s bandwidth. The operating system in service is CentOS 6.5, and OpenMPI 1.8.5 is the used implementation of the MPI-3 standard for communications. §2.1 Introduction 61 196 APPENDIX C. COMPUTING RESOURCES Figure C.1: JFF supercomputer. Figure C.2: MareNostrum supercomputer. Figure 2.1: UPC’s JFF supercomputer consists of 168 computer nodes with 2024 cores and 4.6 TB of RAM in total. 2.1.4 Concluding remarks The development of the new optimization framework has to be carried out taking into account all requirements stated in the previous sections. Most of them are common features expected from a library. Regarding the specific requirements, they may be summarized by saying that the new library must perfectly fit the current optimization needs of CTTC and be fully compatible with the TermoFluids software. For finishing this section dedicated to design concepts, some concluding remarks are listed below: • According to the specific design requirements, the programming language to be used is C++ and the communications standard for parallel computing is MPI (Message Passing Interface). C++ is a widespread, portable and performant object-oriented language. Although Java, for instance, includes a lot of interesting concepts related to parallelism, its overall performance is worse. MPI is also a widely used and portable communication standard, which provides a convenient mechanism for modularizing parallelism through the use of “communicators” 62 §2.2 State of the art of optimization libraries [3]. Thus, both C++ and MPI seem to be good choices for the new optimization framework. • Fulfilling the portability requirement means that the code should run in different computer architectures and operating systems. For the development of this Doctoral Thesis, it has been considered enough that the code runs with the hardware and software stated in the previous section. • The coupling between the optimization framework and the models to be optimized has not been defined. As a first and most universal approach, it has been decided to use a direct coupling, i.e. the model to be optimized will behave as a black box for the optimizer, receiving values for the optimization variables from the optimizer and returning the values of the objective functions. • The user-friendliness of the framework is an important feature. However, it is easier to first focus on the algorithmic of the library and to take care of userfriendliness once an acceptable development level of the library has been reached. Thus, the creation of a graphical user interface, documentation, etc. have not been considered in the scope of this Doctoral Thesis. • The possibility of using third-party libraries has been considered crucial due to the limited time available for the completion of the Doctoral Thesis, the scope of the topic and the fact that optimization theory has been studied for several decades. • The optimization framework will include single-objective and multi-objective capabilities for real-valued problems. In the scope of this Doctoral Thesis, a global search method (a genetic algorithm) and a local search method (a gradient method) will be implemented. Discrete-valued problems are considered of less importance for the kind of CFD & HT studies carried out at CTTC, so operators for discrete optimization will not be implemented for the moment. 2.2 State of the art of optimization libraries 2.2.1 Description of the project’s needs The possibility to reuse already existing code is a good chance of accelerating the development of the new optimization framework. Plenty of scientific and mathematical §2.2 State of the art of optimization libraries 63 software has been released in the last years, easily accessible thanks to Internet. Part of this development has been supported by governmental initiatives through public research funding. However, the considerable amount of available software makes selection a hard task unless the project’s needs are well defined. The available software may be divided in proprietary software (usually commercial software) and open-source software (usually free software). The advantages of proprietary software are that the distributor offers a guarantee and technical support to the customer, and also that the implemented algorithms are usually a result of already settled knowledge. The main disadvantage is that the customer is unable to see the source code and, of course, a certain amount of money is to be paid for the license. On the other hand, open-source software is usually distributed with no guarantee or technical support, but the code is for free and may be accessed with no restriction. It was explained in the previous section that only open-source software with a GNU LGPL or equivalent license is accepted to be used in this project. Fortunately, many open-source initiatives for developing mathematical and scientific libraries have already been carried out and may be of great value. Proprietary software cannot be used for the development of the new optimization framework, but it can serve as a benchmark for validation tests and offers a good description of the interfaces and state of the art to which the research community and industry are used to. The needs expected to be fulfilled by means of a third-party open-source library, and consequently the selection criteria, are the following: • A reference library containing a genetic algorithm for both single-objective and multi-objective real-valued optimization is needed. • A reference library containing a gradient-based single-objective local search method is needed. • Additional features like other optimization algorithms or mathematical methods will be positively evaluated, provided that the quality of the required algorithms (a genetic algorithm and a gradient-based method) is not decreased. • All libraries must be written in C or C++ and compile under Linux operating system. • A good object-oriented structure, abstraction level and readability will be appreciated. 64 §2.2 State of the art of optimization libraries • Although no guarantee may be expected from open-source software, the developer of the selected library must be a trustworthy researcher or research institution. 2.2.2 Free open-source software References to several interesting open-source software packages have been found in the literature. In [4] (published in 2007) a comparison of 6 multi-objective optimization frameworks is provided according to the following criteria: available metaheuristic(s), framework type (black-box or white-box), available metrics, available hybrid algorithms, programming language and parallel features. [5] (published in 1999) gives an overview of sequential and parallel genetic algorithms that existed at the time. [1] (published in 2004) presents a list of other 6 existing optimization frameworks. Although no more articles will be cited here, additional references to optimization software packages may be surely found. However, the articles mentioned before were written some years ago and may not be representative of the last advances in the field of mathematical optimization. A search has been carried out in the Internet seeking for more up-to-date opensource software. A list with the most promising codes that were found, together with a short review, is attached hereafter. Evolving Objects (EO) [6] / Parallel and Distributed Evolving Objects (Paradiseo) [1,7] EO and its further evolution Paradiseo are template-based, ANSI-C++ evolutionary computation libraries which help write stochastic optimization algorithms very fast. The framework is the result of a European joint work and allows finding solutions to all kind of hard optimization problems, from continuous to combinatorial ones. Its main characteristics are a flexible object-oriented design, portability, availability of evolutionary (including genetic algorithms) and discrete local search methods and parallelization options. The framework’s development was active until 2013 and took place at the University of Granada first and at INRIA (Institut National de Recherche en Informatique et en Automatique) finally. Good documentation is available and the library has a GNU LGPL like license. It can be compiled under Linux. GAlib [8] It is a set of C++ genetic algorithm objects developed at the Massachusetts Institute of Technology (MIT). GAlib is known to be a mature code with very efficient operators and §2.2 State of the art of optimization libraries 65 good documentation. The library may be compiled under Linux and seems to have a nice graphical interface. The original source code copyright is owned by MIT, but allows modification and licensing fulfilling some restrictions. The last release of the code took place in 2007. GAUL [9] The Genetic Algorithm Utility Library is a flexible programming library that implements genetic algorithms, apart from other stochastic, evolutionary and local search methods. The code has been designed for Linux, written in C language and parallelized using several standards, including MPI. The last release of the library dates from 2009. GAUL is a trademark of Stewart Adcock and is distributed under the GNU GPL license. GENEVA [10] Geneva is a software library which enables users to solve large scale optimization problems in parallel on devices ranging from multi-processor machines over clusters to Grids and Cloud installations. It currently supports evolutionary algorithms (including genetic algorithms), swarm algorithms, gradient descents and a form of simulated annealing. Performance and extensibility are at the core of Geneva’s object-oriented design. The code is written in C++, may be compiled in Linux environments and is distributed under a GNU Affero GPL v3 license, although there are additional licensing options available. Good documentation is provided with the software, whose last release dates from 2015. Geneva was developed and is maintained by Gemfony scientific, a spin-off from Karlsruhe Institute of Technology (KIT). DAKOTA [11] The DAKOTA toolkit is a software framework for systems analysis, encompassing optimization, parameter estimation, uncertainty quantification, design of computer experiments, and sensitivity analysis. It interfaces with a variety of simulation codes from a range of engineering disciplines, and it manages the complexities of a broad suite of capabilities through the use of object-oriented abstraction, class hierarchies, and polymorphism. The library includes many algorithms (also genetic algorithms), is written in C++, implements the MPI standard and may be compiled under Linux. The code has been developed at Sandia National Laboratories and its distribution is restricted by the GNU LGPL license. The last available version was released in 2016. 66 §2.2 State of the art of optimization libraries TRILINOS [12] The Trilinos Project is an effort to facilitate the design, development, integration and ongoing support of mathematical software libraries. That effort is particularly focused on developing parallel solver algorithms and libraries within an object-oriented software framework for the solution of large-scale, complex multi-physics engineering and scientific applications. Among the many packages conforming Trilinos, attention has been put on MOOCHO (Multifunctional Object-Oriented arCHitecture for Optimization). It is designed to solve large-scale, equality and inequality nonlinearly constrained, nonconvex optimization problems (i.e. nonlinear programs) using reduced-space successive quadratic programming (SQP) methods. The library is written in C++ and may be compiled under Linux. The code has been developed at Sandia National Laboratories and its distribution is restricted by the GNU LGPL license. The last available version was released in 2016. 2.2.3 Proprietary commercial software The search of proprietary software has been less intensive because it does not provide a source code which may be reused for creating the new optimization framework, and also because its use as a benchmark is not foreseen in the scope of this Doctoral Thesis. Nevertheless, a few renowned commercial codes are introduced next. IOSO [13] Indirect Optimization Based Upon Self-Organization (IOSO) is a new generation multidimensional nonlinear optimization software based on the response surface technology. Its strategy differs significantly from other well-known approaches to optimization, apparently increasing the efficiency and robustness with respect to standard algorithms. The software has been designed for dealing with heavy tasks with up to 100 variables and 20 objectives, and provides full-automatic optimization algorithms which do not need to be tuned up by the user. It also provides coupling interfaces with the most widespread CAD, CFD and FEA codes. MATLAB [14] This is a well-known scientific software commercialized by Mathworks which includes plenty of mathematical tools. There are two toolboxes regarding optimization, namely the Optimization Toolbox and the Global Optimization Toolbox. The Optimization Toolbox provides functions for finding parameters that minimize or maximize objectives §2.2 State of the art of optimization libraries 67 while satisfying constraints, including solvers for linear programming, mixed-integer linear programming, quadratic programming, nonlinear optimization, and nonlinear least squares. The Global Optimization Toolbox provides methods that search for global solutions to problems that contain multiple maxima or minima. It includes global search, multi-start, pattern search, genetic algorithm, and simulated annealing solvers. LINGO [15] This is a tool designed for building and solving linear, nonlinear (convex & nonconvex/- global), quadratic, quadratically constrained, second order cone, semi-definite, stochastic, and integer optimization models. LINGO provides a completely integrated package that includes a powerful language for expressing optimization models, a full featured environment for building and editing problems, and a set of fast built-in solvers. Apparently, the optimization software developed by Lindo Systems Inc. is in use at over half the Fortune 500 companies in the US, including 23 of the top 25. 2.2.4 Concluding remarks After having reviewed the state of the art concerning optimization software, it is time for deciding which libraries will be taken as a reference for the development of the new framework. The main features of the 6 previously mentioned open-source libraries are compared in the following lines. GAUL and EO/Paradiseo are similar software packages, but GAUL has the drawback of being protected by a GNU GPL license and of being written in C, not in C++. Its website affirms that the code is used in several universities, although no list of such institutions is provided. Moreover, the last release of the code took place in 2009, four years before the last release of EO/Paradiseo. Thus, it is understood that EO/Paradiseo outperforms GAUL’s capacities and the latter is discarded. GAlib seems to be a robust library, is written in C++ and has a graphical user interface. However, its copyright is owned by MIT, the last release was in 2007 and only includes genetic algorithms. Since a wider scope is covered by other packages, which are also subject to less restrictive licenses, GAlib is discarded. GENEVA offers many interesting features: good documentation, several optimization methods, software parallelization, etc. The only drawback at first glance is its GNU GPL license. EO/Paradiseo holds many interesting characteristics as well: genetic algorithms, local search methods, parallelization, good documentation and tutorials, etc. It is also written in C++, seems to be well structured and is subject to the GNU LGPL 74 §2.3 Main features of Optimus Crossover Crossover consists on recombining the genetic material of n parent individuals, although no more than 2 parents are used in Optimus. Several crossover operators have been proposed in the literature, and it is not usually known which one performs best for a certain optimization problem. This is why the possibility of combining various methods is offered in Optimus, with the aim of using the strengths of each method. The user just needs to define the desired crossover operators and to assign them a probability to be selected. Every time the genetic algorithm needs to cross over the chromosomes of two parents, one among the available crossover operators will be called by means of a roulette wheel selection. 3 crossover operators are implemented in Optimus for standard real-valued individuals: •Hypercube crossover: Offspring are uniformly generated on the hypercube whose diagonal is the segment joining both parents, i.e. by doing linear combinations of each variable independently. The user provides an alpha parameter at the beginning of the optimization, and the crossover operator generates uniformly a random number in the range [alpha,1+alpha] for each variable. This random number is the coefficient used for defining the linear combination of the parents’ values of that variable. •Segment crossover: Offspring are uniformly generated on the segment joining both parents, i.e. the operator constructs two linear combinations of the parents with a single random number uniformly generated in the range [alpha,1+alpha]. Alpha is a parameter provided by the user at the beginning of the optimization, and the random number is the coefficient used for defining the linear combination of the parents’ values for all variables. •Uniform crossover: This operator simply exchanges values of variables between the 2 parents, creating new offspring. In the case of using self-adaptive individuals, two crossover operators are to be selected for recombining 2 individuals: the first operator defines how to recombine the optimization variables of the individuals, whereas the second operator defines how to recombine the self-adaption parameters. The crossover operators available are uniform crossover and hypercube crossover (with the alpha coefficient fixed to 0 value), being the use of segment crossover not enabled for self-adaptive individuals. The roulette-wheel §2.3 Main features of Optimus 75 selection between various crossover operators is neither enabled, so a single operator is to be chosen for each part of the individual. Two parents are needed to apply a standard crossover operator in Optimus, either when standard or self-adaptive individuals are used. However, an additional option called global crossover is available for self-adaptive individuals, which consists on randomly selecting two parents for each gene of the offspring that will be created. Mutation Mutation consists on altering a certain percentage of genes of each offspring created by the crossover operator. Since the best performing operator is not usually known (as it happens in the case of crossover operators), Optimus offers the possibility to combine various methods to use the strengths of each. The user activates the desired mutation operators and assigns to each of them a probability to be selected. Every time the genetic algorithm needs to mutate the chromosome of an offspring, one among the available operators is chosen by means of a roulette wheel selection. A single mutation operator is available in Optimus for self-adaptive individuals: •Self-adaptive mutation: A normal mutation is applied to each optimization variable. Each of these mutations is defined according to a normal distribution centered in the original variable’s value and whose standard deviation is taken from the corresponding self-adaptation parameter carried by the individual. A mutation is also applied to each self-adaptation parameter. 3 mutation operators are implemented in Optimus for standard real-valued individuals: •Uniform mutation: This operator modifies all variables by choosing new values uniformly on an interval centered on the old value and of width 2 · epsilon, being epsilon defined by the user for each optimization variable. •Deterministic-uniform mutation: Exactly k variables are modified uniformly by choosing a new value from an interval centered on the old value and of width 2 · epsilon, being epsilon defined by the user for each optimization variable. •Normal mutation: Also called Gaussian mutation, this operator acts on every optimization variable of an individual creating new values according to a normal distribution centered in the original value and with a fix standard deviation parameter defined by the user. 76 §2.3 Main features of Optimus Selection The selection step consists on choosing the individuals that will be used to generate the offspring population, being its size fixed at the beginning of the optimization. As it was already said, the general behavior is that the better (fitter) an individual, the higher its chance of being selected. The selection operators implemented in Optimus for single-objective optimization (either with standard or self-adaptive individuals) are the following: •Deterministic tournament: This operator returns the best of T uniformly chosen individuals in the population. The number of tournaments to be carried out is equal to the number of needed parents. •Stochastic tournament: The operator chooses uniformly two individuals from the population and returns the best one with probability R (tournament rate), being the real parameter R in the range [0 . 5 , 1 . 0]. The number of tournaments to be carried out is equal to the number of needed parents. Note that a stochastic tournament with rate 1.0 is strictly identical to a deterministic tournament of size 2. •Roulette wheel: This is the classical selection method used by Goldberg [21] in which each parent is selected according to a probability proportional to its fitness. •Ranking: This method starts by assigning a worth, i.e. a modified fitness value, to each individual of the population. Then selection is carried out by means of a roulette wheel algorithm based on the previously calculated worth values. Two parameters are needed for the worth calculation: the pressure (ranging in (0 , 1]) and the exponent (always greater than 0). Worth values are contained in [ m,M ], where m= 2 −pressure / populationSize and M=pressure / populationSize . Inside these bounds, the spacing between worth values depends on the exponent. •Ordered sequential selection: This operator sorts the population from best to worst and returns as many individuals as required following the list. If the population is exhausted and more individuals are needed, it loops back to the beginning of the list and continues returning individuals until the required number of parents is satisfied. If the number of required parents is smaller than the size of the source population, the best individuals are selected once. If the number required parents is N times that of the source size, all individuals are selected exactly Ntimes. §2.3 Main features of Optimus 77 •Unordered sequential selection: This operator shuffles the population and returns as many individuals as required. If the population is exhausted and more individuals are needed, it loops back to the beginning of the list and continues returning individuals until the required number of parents is satisfied. The three multi-objective optimization algorithms included in Optimus (NSGA-II, SPEA2 and IBEA) use a binary deterministic tournament selection, so other selection methods have not been made available for multi-objective optimization to date. Replacement The replacement operator is applied after the birth of all offspring and consists on selecting the survivors from the current and offspring populations in some arbitrary way. In Optimus the population size is always kept constant from one generation to the next one, being the possibility of having a variable population size not implemented. The available replacement operators are valid for either standard or self-adaptive individuals and may be classified according to the following schemes: merge-reduce operators and reduce-merge operators. The merge-reduce scheme has two major steps, first merging both populations of parents and offspring and then reducing that big population to the right size. In the reduce-merge scheme parents are first reduced of the exact number of offspring and then merged with the offspring population. In the latter case, it is implicitly assumed that few offspring have been generated, although this is not mandatory. 3 merge-reduce operators have been implemented: •Comma replacement: This operator, which is common in Evolution Strategies, selects the best offspring and discards all parents. Hence, at least as many offspring as the size of the population must be created. •Plus replacement: It first merges the offspring and the parents and finally the best individuals among them become the next generation. This operator is also common in Evolution Strategies. •EP tournament: This is a classical operator in Evolutionary Programming (EP). First the offspring and parents are merged and then a global tournament of size T begins. It works by assigning a score to all individuals in the population. Starting with a score of 0, each individual I is opposed T times to a uniformly chosen individual. I’s score is incremented by 1 every time it wins, and by 0.5 every 78 §2.3 Main features of Optimus time it draws. Once all tournaments are finished, the individuals for the next generation are selected deterministically based on their scores. 3 reduce-merge operators have been implemented: •Worst replacement: The worst parents are replaced by all offspring. •Deterministic tournament: Each parent to be replaced is chosen by an inverse deterministic tournament, i.e. the operator returns the worst of T uniformly chosen individuals in the population. •Stochastic tournament: Each parent to be replaced is selected by an inverse binary stochastic tournament, i.e. each time the operator chooses uniformly two individuals from the population and returns the worst one with probability R (tournament rate), being the real parameter Rin the range [0.5,1.0]. The possibility of activating weak elitism is also implemented for single-objective optimization, which means that if the best fitness in the new population is worse than the best fitness of the parent population, the worst individual of the new population is replaced by the best parent. This strategy ensures that the overall best fitness in the population will never decrease. 2.3.5 Hybrid methods One single hybrid algorithm has been implemented in Optimus, only available for single-objective optimizations. It combines two constitutive algorithms sequentially: the genetic algorithm as global search method and the Trilinos/Moocho package as local search method. Two approaches to sequential hybrid algorithms have been found in the literature [22]: i) the best solution found by the GA is taken as the starting point for the local search method, and ii) the gradient method can be incorporated in the GA as a new operator and can be applied either to the best individual or to the individuals corresponding to local optima. The first approach has been implemented in Optimus. The control algorithm in charge of switching automatically from one method to the other is represented in the flow diagram in Fig. 2.2. The user must define the stagnation criterion for the GA, i.e. the maximum allowed number of generations (let us refer to it as N ) to get some improvement of the best known solution so far. If no improvement is obtained after N generations, the control algorithm checks if the best individual so §2.3 Main features of Optimus 79 far has been previously used for starting a local search. If so, it means that the local search was unable to improve it and, since it is a deterministic method, exactly the same result would be reached again. Thus, the local search is skipped and the flow returns to the genetic algorithm. If no local search was started in the past with that individual, it is launched (i.e. the Trilinos/Moocho package is called) and the obtained result is compared with the original individual. If the local search result is better, then the best individual in the population is replaced by the newly found solution. After this step, the flow returns to the genetic algorithm and continuation criteria area checked, as usual. The Trilinos/Moocho package (Multifunctional Object-Oriented arCHitecture for Optimization) has been designed to solve large-scale optimization problems using reduced-space successive quadratic programming (SQP) methods, as it was already explained in a previous section. Moocho transforms the original problem into an equivalent quadratic problem, having both of them the same solution. The advantage obtained with the transformation is a higher convergence rate. The new problem could be solved by means of a BFGS algorithm. Nevertheless, the lBFGS algorithm is preferred instead in order to decrease the memory requirement. This issue is of special concern in the case of problems having many variables. The only drawback of lBFGS with respect to BFGS is its linear convergence rate, whereas higher convergence rates are achieved by BFGS. No hybrid method is implemented for multi-objective optimization to date, but an expansion of the framework in this direction will be considered in the future. 2.3.6 Continuation criteria The continuation criteria (also called stopping criteria) are in charge of controlling if the optimization process shall end or go on. Several criteria were implemented in Optimus, being possible to enable more than one in the same simulation. In that case, a single stop signal sent by any active criterion is enough to make the optimization finish. The continuation criteria available in Optimus are described hereafter: • Stop optimization if ... –Resource limitation criteria * ... maximum number of generations reached: When the genetic algorithm has carried out a certain number of generations, the best solution so far is returned as the global optimum. 80 §2.3 Main features of Optimus YES Optimization start Optimal solution found Definition of the optimization problem Selection of GA parameters Creation of an initial population Evaluation of the initial population Selection of parents for mating Application of variation operators: crossover and mutation Evaluation of new individuals Replacement strategy Is any stopping criterion fulfilled? NO Best individual improved in last N generations? YES Best individual already used by LS? Local search Best individual improved? Insert indiv. into population YES NO YES NO YES NO Figure 2.2: Flow chart of a sequential hybrid algorithm, composed by a genetic algorithm and a local search method. §2.3 Main features of Optimus 81 * ... maximum number of evaluations reached: When the optimizer has carried out a certain number of objective function evaluations, the best solution so far is returned as the global optimum. * ... maximum simulation time reached: When the optimization process has spent a certain amount of wall clock time, the best solution so far is returned as the global optimum. –Solution’s stagnation criteria * ... best individual did not improve in N generations: When the optimizer has not been able to improve the best solution so far after a certain number of generations, that solution is returned as the global optimum. –Target accomplishment criteria * ... target fitness value reached: When the optimizer has obtained a solution with a fitness value which is equal to or better than a target fitness value established at the beginning of the optimization, the best solution so far is returned. • Continue optimization until ... –Minimum usage of resources criteria * ... minimum number of generations reached: The optimization process cannot be stopped until the genetic algorithm has carried out a certain amount of generations, even though some stopping criterion is fulfilled. 2.3.7 Statistics It is important to provide some information of the optimization process in run time in order to allow the engineer to check regularly the correctness of the simulation. Due to the high number of objective function evaluations usually needed, an early detection of unexpected behaviors might save a considerable amount of computational time. Moreover, the provided information must be enough to evaluate if the obtained result is the best attainable solution for the simulation that has been run. Thus, it was decided to calculate the following data after the genetic algorithm completes each generation: •Number of generations: Total amount of generations carried out so far. 82 §2.3 Main features of Optimus •Number of evaluations: Total amount of objective function evaluations carried out so far. •Wall clock time: Total amount of wall clock time spent so far by the optimization process. •Best fitness value: Fitness value of the best individual available in the current population. In the case of multi-objective optimization, the best value available in the current population for each objective is identified. •Statistics of the population: Average and standard deviation of the fitness, taking into account every individual in the current population. In the case of multi-objective optimization, the average and standard deviation of the current population for each objective is calculated. Some additional metrics are available after the completion of each generation when multi-objective optimizations are carried out. These metrics are contribution and entropy, and are computed for both the archive and the current population. 2.3.8 Parallelization Task management strategy After the thorough state of the art study of parallelization techniques carried out in Chapter 1, it was decided to first implement a two-level model (see Fig. 2.3 (a)): a master-worker strategy for the optimization algorithm in the upper level, and a case dependent fine-grain parallelization of individuals in the lower level. The parallel task manager was implemented from scratch, although some interesting notes about the parallel algorithms used by Paradiseo were found in [1]. The adopted parallelization strategy is used by both the genetic algorithm and the local search method and relies exclusively on hardware parallelization, thus obtaining in a shorter time identical results as sequential algorithms. This may be very beneficial when applied to optimizations with a high computational cost, as those in the field of CFD & HT. A similar approach is described in [22], where the problem of getting an optimum shape design of aerodynamic configurations is studied. Once that satisfactory results have been achieved, a 3-level parallelization model composed by an upper coarsegrain level (island model) and an intermediate master-worker level will be implemented for the genetic algorithm, keeping the fine-grain parallelization of individuals in the lower level (see Fig. 2.3 (b)). §2.3 Main features of Optimus 83 - 9In a parallel GA there exist many elementary GAs working on separate sub-populations Pi(t). Each sub-algorithm includes an additional phase of periodic communication with a set of neighboring sub-algorithms located on some topology. This communication usually consists in exchanging a set of individuals, although nothing prevents the sub-algorithms of exchanging other kind of information such as population statistics. All the sub-algorithms are thought to perform the same reproductive plan. Otherwise the PGA is heterogeneous [2], [34] (even the representation could differ among the islands, posing new challenges to the exchange of individuals [44]). In a distributed GA demes are loosely-coupled islands of strings (Figure 8b). A cellular GA (Figure 8c) defines a NEWS neighborhood (North-East-West-South in a toroidal grid) in which overlapping demes of 5 strings (4+1) execute the same reproductive plan. In every neighborhood (deme) the new string computed after selection, crossover, and mutation replaces the current one only if it is better (binary tournament), although many other variants are possible. This process is repeated for all the neighborhoods in the grid of the cellular GA (there are as many neighborhoods as strings). With regard to the classes of Figure 5 we suggest in Figure 8 several implementations of PGAs. Global parallelization consists in evaluating, and maybe crossing+mutating in parallel all the structures, while selection uses the whole population. Interesting considerations on global parallelization can be found in [16]. This model provides lower runtime only for slow objective functions; an additional limitation is that the search mechanism uses a single population. The automatic parallelization is rarelly found since the compiler must provide the parallelization of the algorithm automatically. The other hybrid models in Figure 8 combine different parallel GAs at two levels in order to enhance the search in some way. Interesting surveys on these and other parallel models can be found in [1], [16], [17]. Hierarchies of GAs are the most recurrent models found in the literature. In Figure 8d we can appreciate a distributed algorithm in which every island runs a cellular GA. In Figure 8e several algorithms using global parallelization are used to create a ring of islands. Finally, Figure 8f shows two levels of coarse grain PGAs, the inner level having a full-connected topology, and the outer level having a simple ring topology. ... Master Workers (e) (f)(a) (b) (d) Figure 8. Different models of PGA: (a) global parallelization, (b) coarse grain, and (c) fine grain. Many hybrids have been defined by combining PGAs at two levels: (d) coarse and fine grain, (e) coarse grain and global parallelization, and (f) coarse grain plus coarse grain. We want to point out that coarse (cgPGA) and fine grain (fgPGA) PGAs are subclasses of the same kind of parallel GA consisting in a set of communicating sub-algorithms. We propose a change in the nomenclature to call them distributed and cellular GAs (dGA and cGA), since the grain is usually intended to refer to their computation/communication ratio, while actual differences can also be found in the way in which they both structure their population (see Figure 9). While a distributed GA has a large sub-population (>>1) a cGA has typically only one string in every sub-algorithm. For a dGA the sub-algorithms are loosely connected, while for a cGA they are tightly connected. In addition, in a dGA there exist only a few sub-algorithms, while in a cGA there is a large number of them. Figure 2.3: PGA models in Optimus (extracted from [5]): (a) global parallelization (current implementation), (b) coarse grain + global hybrid (future implementation). The implemented 2-level strategy is similar to the self-scheduling strategy described in [23], where the master processor manages a single processing queue and maintains a prescribed number of jobs active on each group of workers. Once a group of workers has completed a job and returned its results, the master assigns the next job to this group. Thus, the workers themselves determine the schedule through their job completion speed. Heterogeneous processor speeds and/or job lengths are naturally handled, provided there are sufficient instances scheduled to balance the variation. Individuals belonging to the genetic algorithm’s population compose the batch of individuals to be evaluated in parallel most of the times. Note that a synchronous genetic algorithm has been implemented, i.e. all individuals belonging to a generation are evaluated before the next generation starts. However, when the local search method is being run, a batch of individuals is evaluated every time a gradient needs to be computed. 3 parameters are expected from the user of the optimization library in order to configure the parallel task manager: the number of available processors, the number of processors that will be used to evaluate each individual, and the processor partitioning model [23]: “dedicated master” or “peer partition” approach. The master processor is dedicated exclusively to task scheduling operations in the first approach, whereas it also participates in the computation of individuals in the latter approach. Note that the peer partition approach shall be used cautiously in order to avoid inter-processor communication delays, as it is explained further on in this section. Communications strategy Hardware may use two types of memory: shared memory and distributed memory. On a shared memory multiprocessor, information (e.g. the population) is stored in 90 §2.3 Main features of Optimus –Parallelization: The parameters to be defined are the processor partitioning model (dedicated master or peer partition strategy) and the number of processors for solving each individual. –Recovery information: Since optimization may be a costly process, it is crucial to regularly store some recovery information from which the optimization may be restarted in case the simulation is interrupted unexpectedly. This recovery information consists of two files: a file containing all the parameters defining the optimization algorithm and another file containing the simulation results obtained so far. If the optimization is interrupted, it is possible to restart it using these two files and exactly the same result as in the uninterrupted simulation will be obtained. The user may choose the names of the two files, and also the frequency with which the file containing simulation results is saved. When a simulation is wanted to be restarted using the information of the two recovery files, the names of the files to be used are specified inside this parameter group, as well. –Continuation criteria: The user may activate the desired continuation criteria among the available ones, which were enumerated in a previous section. –Variation operators: The parameters contained in this group are the bounds for the optimization variables or genes, the types and application probability of crossover operators, and finally the types and application probability of mutation operators. •Local search’s parameter sheet This file contains every parameter controlling the behavior of the local search algorithm. The Trilinos/Moocho package is the only available local search algorithm to date, so the format of its original parameter file, called Moocho.opt by default, is used. The parameters the user may define are the following: maximum number of iterations, maximum run time, convergence tolerance, level of detail of the generated outputs and some mathematical options specifying the behavior of the local search algorithm. Additional information is available in the Trilinos Project’s documentation [12]. Data output interface Simulation information may be extracted either onto the screen or into the hard drive. The purpose of this design is to allow the user to follow the execution of the optimization §2.3 Main features of Optimus 91 algorithm by reading the information shown on screen and, once the simulation has finished, to have extended information saved in the hard drive for post processing and later use. The user may customize the output information, as it was explained in section Data input interface. The hierarchical data structure of the output files is the following: • Optimization case directory –OptimusFiles directory –Results directory The user is expected to create a directory (the optimization case directory) in which every file related to the optimization case is stored. All output files created by the Trilinos/Moocho package are stored here, together with the parameter file used for recovery purposes. The subdirectory OptimusFiles contains as many subdirectories as individuals that can be simulated simultaneously. The files generated for solving the objective function(s) of each individual are saved in those subdirectories. The subdirectory Results, although it may be renamed by the user, contains the files with general simulation results and also the recovery files. 2.3.10 Other features Some additional aspects of the Optimus library are mentioned in the following paragraphs. Regarding the random number generator, the widely used Mersenne Twister MT19937 [24] has been implemented, already mentioned in Chapter 1 and available in the Paradiseo package. The use of surrogate models for reducing the evaluation time of computationally expensive objective functions is very interesting. It has not been implemented due to lack of time, but it is planned in a future extension of the library. Finally, some popular mathematical functions used as benchmarks for optimizers have been implemented. On one hand, the Rosenbrock, Rastrigin, Schwefel and Griewank functions are available for testing the single-objective optimization algorithms. On the other hand, Kursawe’s function (MOP4) and Zitzler-Deb-Thiele’s function N. 6 (ZDT6) were included for testing multi-objective optimization algorithms. These functions are intended to be used every time modifications are inserted into the library in order to check the correctness of the new code. 92 §2.4 Validation tests 2.3.11 Optimus vs. Paradiseo The Optimus library is clearly based on the structure and algorithms available in Paradiseo. Therefore, which is the advantage of using Optimus? The summary of the main differences between both libraries to date is the following: • The local search methods available in Paradiseo are suitable for discrete optimization, but no alternative is provided for real-valued optimization. A gradient-based local search method was added in Optimus, making possible to build a hybrid method (genetic algorithm + gradient method) for real valued optimization problems. • The existence of at least two development branches is very noticeable in Paradiseo. An important consequence is that the use of self-adaptive individuals is possible in single-objective optimization, but not in multi-objective optimization. Since the use of this kind of more sophisticated individuals clearly improves the obtained results, self-adaptive individuals for multi-objective optimizations have been implemented in Optimus. • The parallelization strategy in Optimus and Paradiseo differs considerably. Paradiseo is designed focusing on combinatorics problems whose objective function evaluation times are usually low. The master-worker parallelization algorithm in Optimus has been designed for objective functions having long simulation times, thus offering the possibility of evaluating each function in parallel. The parallel code of Optimus has been created from scratch, so any potential similarity with Paradiseo has occurred by chance. 2.4 Validation tests Now that the first version of the optimization library has been created, it is necessary to test the correctness of the implemented algorithms. An exhaustive check of every feature was carried out by the author. Most of the tests were more related to programming issues and are not shown in this Doctoral Thesis. The validation tests related to the conceptual development of the library are presented in 2 steps: i) based on mathematical functions contained in common optimization test suites, and ii) based on real-world simulations belonging to the field of Computational Fluid Dynamics & Heat Transfer (CFD & HT). §2.4 Validation tests 93 Validation tests measure the accuracy and performance of the optimization library. Accuracy is calculated by comparing the known optimal solution of the benchmark case with the solution provided by the optimization algorithm. The performance may be obtained in several ways: measuring the number of objective function evaluations, measuring simulation time, etc. At this development stage, accuracy is considered important and little attention is paid to performance, i.e. no comparison of different methods and parametrizations of the optimization algorithm is carried out. The reason is that no effort was made so far to optimize the performance of the code, so some measurements are shown but just for information purposes. The unique goal is to provide a few results which prove the correctness of the library’s implementation. 2.4.1 Benchmark mathematical functions In Chapter 1 several reputed test suites were mentioned, together with a selection of some functions for this Doctoral Thesis. On one hand, the Rosenbrock, Rastrigin, Schwefel and Griewank functions were selected for testing the single-objective optimization algorithms. On the other hand, Kursawe’s function (MOP4) and Zitzler-Deb-Thiele’s function N. 6 (ZDT6) were the preferred functions for testing multi-objective optimization algorithms. Since the evaluation time of these mathematical functions is extremely short, all simulations were run sequentially in 1 processor of the supercomputer and hence the parallel features of Optimus were not used. Single-objective tests The selected test functions were simulated with 2, 4 and 6 optimization variables. The stopping criterion for the optimizer was to reach a target objective value of 1e-9, being the real optimal value equal to 0 for all functions. Each function was optimized twice: once using exclusively the genetic algorithm, and once using the hybrid algorithm composed by the genetic algorithm and the local search method in Trilinos/Moocho package. Results are shown in Table 2.2 and include the obtained optimal objective value, the number of function evaluations carried out and the required wall clock time for each optimization. Due to the decimal precision used in the implementation of Schwefel’s function, the optimizer was unable to reach the required accuracy. Nevertheless, the obtained solutions are considered to be satisfactory. Regarding the Griewank’s test function, note that the results obtained by the hybrid optimizer are not shown. The reason is that it was unable to find the global optimum of 94 §2.4 Validation tests the function with the required accuracy. This fact is explained by the noisy nature of the function’s fitness landscape. Note also that the optimization variables were bounded in the range [−50,50]. Test function Benchmark optimum Genetic algorithm Hybrid algorithm Optimum No. function evaluations Time (seconds) Optimum No. function evaluations Time (seconds) Rosenbrock_ 2vars 0 6.6658e-10 50048 6.7528 8.8794e-12 169 0.0172 Rosenbrock_ 4vars 0 9.9656e-10 59920 7.1141 2.0650e-11 1031 0.1171 Rosenbrock_ 6vars 0 9.8136e-10 95270 8.9573 2.5048e-11 4120 0.4481 Rastrigin_ 2vars 0 4.6314e-10 1400 0.2051 0 332 0.0550 Rastrigin_ 4vars 0 6.9260e-10 5360 0.6044 0 2254 0.2403 Rastrigin_ 6vars 0 7.8748e-10 11680 0.9840 0 6859 0.6079 Schwefel_ 2vars 0 2.5455e-05 1376 0.1951 2.5455e-05 215 0.0668 Schwefel_ 4vars 0 5.0910e-05 5450 0.6062 5.0910e-05 1823 0.1825 Schwefel_ 6vars 0 7.6365e-05 14860 1.1752 7.6365e-05 3519 0.3304 Griewank_ 2vars 0 2.2690e-10 1184 0.1669 - - - Griewank_ 4vars 0 7.9315e-10 21220 1.4272 - - - Griewank_ 6vars 0 9.0648e-10 31900 3.0057 - - - Table 2.2: Results of the single-objective test functions. Multi-objective tests Multi-objective test functions MOP4 and ZDT6 are shown respectively in Fig. 2.4 and Fig. 2.5. The true optimal Pareto front (named Benchmark) and the Pareto front obtained by the optimizer (named Optimus) are represented in both figures. 50100 objective function evaluations were needed in order to build each one of the Pareto fronts corresponding to MOP4 and ZDT6. Performing those evaluations took 94.88 seconds in the first case and 28.45 seconds in the latter. It is considered that 1 evaluation involves calculating Objective 1 and Objective 2 in both tests. It can be seen that the accuracy and diversity of the Pareto front found by Optimus is acceptable, although better results could be found by assigning either more time or computational resources. §2.4 Validation tests 95 -12 -8 -4 0 -20 -18 -16 -14 f2(x) f1(x) Benchmark Optimus Figure 2.4: Pareto Front obtained by Optimus for the Kursawe’s function (MOP4). 0 0.5 1 1.5 0.5 0.75 1 f2(x) f1(x) Benchmark Optimus Figure 2.5: Pareto Front obtained by Optimus for the Zitzler-Deb-Thiele’s function N.6 (ZDT6). 96 §2.4 Validation tests 2.4.2 CFD & HT tests Two different real-world cases from the CFD & HT field are presented next. The first case consists on the energy labelling of a fridge, which is based on a real industrial problem solved at CTTC [25]. The second case is a CFD simulation of an incompressible fluid circulating through a rectangular pipe, where the optimal geometry of the pipe is searched. Optimization of the energy efficiency index of a fridge The energy efficiency index of a fridge is calculated according to 2 main characteristics: the energy consumption of the fridge and its useful (internal) volume. The lower the efficiency index, the fridge’s energy labelling is better and its price in the market higher. However, the optimal relation between the energy consumption and the useful volume is not trivial. The reason is that having thinner walls, i.e. greater internal volume, increases the heat losses due to lack of thermal isolation. The task of the optimization algorithm consists on finding the optimal balance between these two characteristics for a set of 3 problems of industrial interest. The first step has been to create a simplified mathematical model of the fridge, which may be solved both analytically and using more advanced computational tools. In this model, the fridge is composed by N walls. Each wall has the following characteristics: Ai:area ki:thermal conductivity di:thickness ∆Ti:temperature gradient The total conduction heat loss through the walls is defined by the following expression: Q= N X i=1 kiAi∆Ti/di(2.1) Two volumes are defined: the external volume of the fridge ( Vo ) and the internal or useful volume ( V ). Vo has a constant fixed value, whereas the internal volume is defined as: V=Vo− N X i=1 Aidi(2.2) §2.4 Validation tests 97 The energy efficiency index ( Ei ) is defined according to the laws in force in 2013 (see Fig. 2.6): Ei=Ea/Est (2.3) where the term Ea is related to the energy consumption of the appliance and Est is related to the internal volume. On one side, Eais defined as follows: Ea=100·k1·Q(2.4) where k1=365·24 1000·COP (2.5) being COP the Coefficient of Performance. On the other side, Est is defined as in Eq. 2.6 for A+ and A++ energy labels and as in Eq. 2.7 otherwise: Est =M·(AV)+N+CH (2.6) Est =M·(AV)+N(2.7) where M and N depend on the appliance class and CH is a correction factor (see Fig. 2.6). AV is defined as in Eq. 2.8 for A+ and A++ energy labels and as in Eq. 2.9 otherwise: AV =XµVc (25−Tc) 20 ·FF ·CC ·BI¶(2.8) AV =VR+ΩVF(2.9) where Tc is the compartment temperature, Vc is the compartment storage volume, FF / CC / BI are several correction factors tabulated in Fig. 2.6, VR is the fresh food storage compartment volume, VF is the frozen food storage compartment value and Ω is dependent on the appliance class (see Fig. 2.6). Two implementations of the described mathematical model were carried out. The first one consists on just writing the formulation using C++ functions. The second one consists on building a multi-physics model of the fridge using CTTC’s in-house multi-physics software, called NEST [19]. The aim of the latter implementation is to 98 §2.4 Validation tests test the coupling between Optimus and NEST, which was one of the design criteria of Optimus. The multi-physics model has a main system, called Fridge, at the top level. This system is composed by 2 subsystems, namely the refrigerator compartment and the freezer compartment. The refrigerator compartment subsystem is composed by 5 walls (named as refrigerator right side, refrigerator left side, refrigerator rear, refrigerator top and refrigerator door), whereas the freezer compartment subsystem is composed by 9 walls (named as bottom, compressor zone right side, compressor zone left side, compressor zone vertical side, compressor zone horizontal size, freezer right side, freezer left side, freezer rear, freezer door). Thus, the overall number of walls is 14 and the total volume and heat losses of the fridge are the sum of the volumes and heat losses of both compartments. Figure 2.6: Tables for calculating fridges’ and freezers’ energy efficiency index ( Ei ) according to the laws in force in 2013. §2.4 Validation tests 99 3 different optimization problems are presented next. The ∆T and the area of each wall are known input parameters, as well as the coefficients tabulated in Fig. 2.6. The aim of the optimizations is always to find a set of wall thickness values corresponding to the optimum design according to an objective function. The objective function is different in all 3 problems. The benchmark results have been obtained by using the analytical method of Lagrange multipliers applied to the simplest implementation of the objective function, i.e. the formulation written using C++ functions. On the other hand, the multi-physics model has been the objective function evaluated by Optimus. Since the proposed optimization problems are extracted from a real world industrial project, real geometrical data, operation points and other confidential information is not provided. This is why every wall thickness value of each optimization problem has been normalized using the biggest thickness in the benchmark results of that problem. Problem 1: Find the set of thickness values which minimizes the energy efficiency index (Ei) of the fridge. min f (d)=Ei(d) s.t.Ei>0 and di>0f or i =1,...,14 where dis a vector containing the 14 diwall thickness values. The obtained optimal sets of thickness values are shown in Table 2.3, whereas the associated heat losses, volumes and energy efficiency indexes are shown in Table 2.4. It took 2464 objective function evaluations to the genetic algorithm to reach those results, whereas it took 584 evaluations to the hybrid algorithm. 106 §2.4 Validation tests 0 100 200 300 400 500 0 0.1 0.2 0.3 0.4 0.5 0.6 Qloss [w] Vinsulation [m3 ] Optimus 20 30 40 50 60 70 80 90 100 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 Qloss [w] Vinsulation [m3 ] Optimus Figure 2.8: Variation of the minimal attainable total heat loss ( Qloss ) as a function of the insulation volume ( Vinsul ). Upper graph: Vinsulation is represented in the whole range. Lower graph: Vinsulation is represented in a smaller range, being every point extracted from the upper graph. §2.4 Validation tests 107 Both graphs were obtained by solving 2-objective optimization problems in which both objectives had to be minimized by assigning appropriate values to the set of 14 thicknesses. The obtained optimal Pareto fronts are shown in the graphs, which represent the best attainable values of Ei and Qloss for each value given to Vinsul . Note that a few points dominated by other points (according to Pareto’s criterion) are shown. This happened because several simulations for different Vinsul ranges were first run and the obtained results were finally assembled giving rise to the Pareto fronts in Fig. 2.7 and Fig. 2.8, being the dominated points not filtered after assembling the initial results. For concluding this section focused on the energy optimization of a fridge, it can be said that the Optimus library solved successfully all presented industrial optimization problems. Moreover, it has been proved that the coupling between Optimus and CTTC’s multi-physics software NEST works fine and is comfortable to use. Surface morphing optimization problem Surface morphing is one of the applications for which an optimization algorithm may be used. The need to find the optimal shape of a component which reduces aerodynamic drag, for instance, is a common problem in the field of aerodynamics. Finding such geometry is a difficult task and carrying out a parametric study may involve a large amount of computational time. The finer the mesh used for modeling the geometry, the more computational time is required. The most interesting aspects of facing such a problem with Optimus are twofold: i) it needs to be coupled to the TermoFluids package, the in-house CFD software at CTTC, and ii) the parallelization algorithm of Optimus can be tested. A two dimensional surface morphing problem is proposed next, based on the case of a duct with variable indentation taken from [26]. The geometric parameters of the mathematical model are shown in Fig. 2.9 and in the following equation, y(x)=     h|x|∈[0,x1] 0.5h[1−tanh(a(|x|−x2) )] |x|∈[x1,x3] 0|x|>x3 (2.12) where a= 4 . 14, x1= 4 b , x3= 6 . 5 b and x2= 0 . 5( x1+x3 ). h is the indentation, bounded in [0 , 0 . 38 ·b ]. b is a parameter representing the height of the duct, and for the present test case it is defined to be b= 1 . 0. The flow at the initial state is assumed to be completely 108 §2.4 Validation tests developed with velocity v=(6y(1−y),0,0)T(2.13) The inlet flow is set to be constant according to Eq. 2.14 U=Zb 0 6y(1−y)dy =1 (2.14) A pressure-based condition is applied at the outlet boundary. The Reynolds number is fixed to 507 and the Prandtl number to 0.71. As the indentation is varied, the velocity contour maps are altered as it is shown in Fig. 2.10. (a) Initial mesh. (b) Moved mesh on the plane XY . (c) Moved mesh on the plane XZ. Figure 6: Three-dimensional deformation of a sphere. (a) Plane XY . (b) Plane XZ. Figure 7: Initial mesh (left), moved mesh (middle) and detail of the moved mesh (right). Figure 8: Geometry of the test case based on a duct with a moving indentation, with b=1. 3.2. Comparison with the spring analogy method and CFD validation 3.2.1. Description of the benchmark case For the following tests (sections 3.2.2 and 3.2.3), the benchmark problem is based on a duct with a moving indentation [35, 36, 37, 38, 39]. The geometric parameters and the moving law are Figure 2.9: Geometry of the test case based on a duct with a moving indentation, with b= 1 (extracted from [26]). Figure 9: Mesh quality evolution with an unstructured mesh of 72.300 CV. the indented wall and behind is higher than the constant inlet rate. Around t∗=0.2 the first bubble recirculation appear. As time goes on other vortices occur behind the first one. They arise at either top and bottom walls and they move with the flow. At the end of the period the flow becomes that of the fully developed flow. Notice that on Figure 10 (see subfigure corresponding with t∗=0.7) eddies A, B, C and D have been pointed out in order to identify them easily below. Figure 10: Velocity contour maps along a cycle of deformation. In Figure 11 the wall shear stress on the indented wall is depicted. It indicates the strength and position of the eddies. Owing to local errors, unstructured meshes implies oscillations (see Figure 11a). Figure 11b illustrates the variation of the wall shear stress obtained with a structured mesh during a cycle. The results are in qualitative agreement with [35, 36]. Moreover, in Figure 12 the predicted position of crests and troughs corresponding to eddies B, C and D are depicted and compared with those experimentally obtained by Pedley and Stephanoff [38] as well as with the numerical results of Ralph and Pedley [39]. According to the latter, the adimensionalized abscissa is defined as x∗= (x−x1)(10St)1/3/b(37) Figure 2.10: Velocity contour maps for different indentation values (extracted from [26]). The proposed optimization problem makes use of this mathematical model and is stated as follows: §2.4 Validation tests 109 Problem 4: Find the optimal indentation of the pipe which minimizes the mean horizontal velocity of the fluid at x=0. min f (h)=Umean(h)at x =0 where 0 ≤h≤0.38 The problem has one single optimization variable and the mean horizontal velocity at x=0 can be calculated as follows: Umean =1 b−hZb h 6y(1−y)dy (2.15) Due to the simplicity of the problem, it is known that the optimal indentation value is h=0. The implementation has been carried out in the following way. A mesh for the rectangular pipe with no indentation has been first done by means of ICEMCFD [27]. When a new mesh is needed for a certain indentation value, a parallel radial basis function method for unstructured dynamic meshes [26] is used in order to displace step by step the nodes of the base mesh until the desired indented geometry is reached. Once the required mesh is constructed, the flow is solved by means of an explicit CFD solver from the TermoFluids package. It is made sure that the simulation time is enough to reach a steady turbulent state before the mean horizontal velocity at x= 0 is calculated. The cross section of the pipe has been discretized in 100 equally sized control volumes for Umean’s calculation. Two optimizations have been carried out using the genetic algorithm and the hybrid algorithm respectively. The same continuation criterion was selected in both cases, i.e. to stop the optimization when the best individual was not improved after a certain number of generations. On one hand, it took 925 function evaluations to the genetic algorithm to reach an objective value of 0.99036. The indentation corresponding to this result was 1.79338e-06. On the other hand, it took 676 function evaluations to the hybrid algorithm to reach an objective value of 0.990356, being the indentation corresponding to this result 2.94003e-05. The theoretical optimal objective value is 1, which is obtained for an indentation value of h= 0. It can be seen that the obtained results are very close, being deviations attributed to the inaccurate discretization of the integral used for calculating Umean. 110 §2.5 Conclusions The obtained results are considered satisfactory, mainly because a successful OptimusTermoFluids coupling was achieved and because the correctness of the parallel task manager’s implementation was proved. Good results were also provided for the optimization problem, although this was a secondary concern due to its simplicity. 2.5 Conclusions A new optimization library named Optimus has been implemented in this second chapter. The first step has been the definition of general and specific requirements of the library. The general requirements are based on common design concepts of a framework, whereas the specific requirements are conditioned by the research field and the infrastructure available at the research institution (CTTC). Thus, a brief description of the available computing facilities has been included. A state of the art study on already available optimization libraries has been conducted next, considering both free open-source software and proprietary commercial software. After comparing several packages, it has been decided to base the development of Optimus in two of them: Paradiseo and Trilinos/Moocho, which are both free opensource packages. Paradiseo has been chosen because of its implementation of genetic algorithms, whereas the availability of gradient-methods suitable for local search has been the appealing feature of Trilinos/Moocho. The main features of Optimus are explained in the subsequent section. The development strategy for the new library is first introduced. The implemented algorithms are described next, grouped under the following sections: definition of the optimization problem, single-objective vs. multi-objective optimization, genetic operators, hybrid methods, continuation criteria, statistics, parallelization, user interface and other features. Finally, a brief comparison between Optimus and Paradiseo is carried out with the aim of remarking the improvements available in the new library. The last section of this chapter contains the results obtained in two types of validation tests: benchmark mathematical functions and specific tests from the field of Computational Fluid Dynamics and Heat Transfer (CFD & HT). On one hand, the simulated mathematical functions have been extracted from common optimization test suites designed for both single-objective and multi-objective cases. On the other hand, CTTC has a vast experience on CFD & HT simulations and two cases have been selected: the optimization of the energy efficiency of a fridge, and the optimization of the geometry of a pipe. Optimus has been able to successfully carry out all tests, which proves that the library is ready to be used for solving real-world optimization problems. References 111 Acknowledgments This work has been financially supported by a FPU Doctoral Grant (Formación de Profesorado Universitario) awarded by the Ministerio de Educación, Cultura y Deporte, Spain (FPU12/06265) and by Termo Fluids S.L. The author thankfully acknowledges these institutions. References [1] S. Cahon, N. Melab, and E. G. Talbi. Building with ParadisEO reusable parallel and distributed evolutionary algorithms. Parallel Computing, 30(5):677–697, 2004. [2] O. Lehmkuhl, C. D. Perez-Segarra, R. Borrell, M. Soria, and A. Oliva. TERMOFLUIDS: A new parallel unstructured CFD code for the simulation of turbulent industrial problems on low cost PC cluster. In Parallel Computational Fluid Dynamics 2007, pages 275–282. Springer, 2009. [3] Message Passing Interface Forum. MPI: A Message-Passing Interface standard, version 3.1. High Performance Computing Center Stuttgart (HLRS), 2015. [4] A. Liefooghe, M. Basseur, L. Jourdan, and E. G. Talbi. ParadisEO-MOEO: A framework for evolutionary multi-objective optimization. In International Conference on Evolutionary Multi-Criterion Optimization (EMO 2007), volume 4403, pages 386–400, Matsushima, Japan, 2007. Springer. [5] E. Alba and J. M. Troya. A survey of parallel distributed genetic algorithms. Complexity, 4(4):31–52, 1999. [6] M. Keijzer, J. J. Merelo, G. Romero, and M. Schoenauer. Evolving objects: A general purpose evolutionary computation library. In Proceedings of the 5th International Conference on Artificial Evolution (EA’01), pages 231–242, Le Creusot, France, 2001. Springer. [7] Paradiseo, a software framework for metaheuristics. http://paradiseo.gforge.inria.fr. [8] M. Wall. GAlib: A C++ library of genetic algorithm components. Mechanical Engineering Department, Massachusetts Institute of Technology, 87:54, 1996. 112 References [9] S. Adcock. GAUL, the Genetic Algorithm Utility Library. http://gaul.sourceforge.net, 2009. [10] R. Berlich, S. Gabriel, A. Garcia, and M. Kunze. Distributed parametric optimization with the Geneva library. In Data Driven e-Science, pages 303–314. Springer, 2011. [11] B. M. Adams, L. E. Bauman, W. J. Bohnhoff, K. R. Dalbey, M. S. Ebeida, J. P. Eddy, M. S. Eldred, P. D. Hough, K. T. Hu, J. D. Jakeman, J. A. Stephens, L. P. Swiler, D. M. Vigil, and T. M. Wildey. Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis: version 6.0 user’s manual. Sandia Technical Report SAND2014-4633, July 2014. Updated November 2015 (Version 6.3). [12] M. A. Heroux, R. A. Bartlett, V. E. Howle, R. J. Hoekstra, J. J. Hu, T. G. Kolda, R. B. Lehoucq, K. R. Long, R. P. Pawlowski, and E. T. Phipps. An overview of the Trilinos project. ACM Transactions on Mathematical Software (TOMS), 31(3):397–423, 2005. [13] IOSO NM Version 1.0, User’s Guide. IOSO Technology Center, Moscow, Russia, 2003. [14] MATLAB Release 2016a. The MathWorks, Inc., Natick, Massachusetts, United States, 2016. [15] LINGO 16.0. Lindo Systems Inc., Chicago, Illinois, United States. [16] H. P. Schwefel. Numerische Optimierung von Computer-Modellen mittels der Evolutionsstrategie, volume 1. Birkhäuser, Basel, Switzerland, 1977. [17] H. P. Schwefel. Internal Report of KFA Juelich, KFA-STE-IB-3/80, 1980. [18] G. Rudolph. Globale Optimierung mit parallelen Evolutionsstrategien. Diploma thesis, University of Dortmund, 1990. [19] R. M. Damle, O. Lehmkuhl, G. Colomer, and I. Rodríguez. Energy simulation of buildings with a modular object-oriented tool. In ISES Solar World Congress 2011, pages 1–11, 2011. [20] A. Liefooghe, L. Jourdan, and E. G. Talbi. A unified model for evolutionary multiobjective optimization and its implementation in a general purpose software framework: ParadisEO-MOEO. Research Report RR-6906, INRIA, 2009. References 113 [21] D. E. Goldberg. Genetic algorithms in search, optimization, and machine learning. Addison-Wesley, 1989. [22] N. Marco and S. Lanteri. A two-level parallelization strategy for genetic algorithms applied to optimum shape design. Parallel Computing, 26(4):377–397, 2000. [23] M. S. Eldred, W. E. Hart, B. D. Schimel, and B. G. van Bloemen Waanders. Multilevel parallelism for optimization on MP computers: Theory and experiment. In Proc. 8th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, number AIAA-2000-4818, Long Beach, CA, volume 292, pages 294–296, 2000. [24] M. Matsumoto and T. Nishimura. Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Transactions on Modeling and Computer Simulation (TOMACS), 8(1):3–30, 1998. [25] Modular domestic refrigeration systems with high energy efficiency (KERS). Research Project F00299, Ref. IPT-020000-2010-30, INNPACTO (Spanish Government), Company: Fagor Electrodomésticos, Funding: 439552 Euros, Period: 20102013. [26] O. Estruch, O. Lehmkuhl, R. Borrell, C. D. Pérez-Segarra, and A. Oliva. A parallel radial basis function interpolation method for unstructured dynamic meshes. Computers & Fluids, 80:44–54, 2013. [27] Inc. ANSYS. ANSYS®ICEM CFD, Release 16.2. 114 References 3 Load balancing methods for parallel optimization algorithms Abstract. The aim of this third chapter is to evaluate the impact of several load balancing methods on the parallel efficiency of optimization algorithms. After carrying out a state of the art study of the available techniques, new load balancing algorithms are proposed and tested by means of an exhaustive theoretical case study. Finally, they are applied to the optimization of a real-world engineering model consisting on the refrigeration of a power electronic device. 115 122 §3.2 Approach to the load balancing problem Time Processor ID 1 2 3 4 (a) Time Processor ID 1 2 3 4 (b) Time Processor ID 1 2 3 4 (c) 1 1 2 3 1 2 3 (a) (b) (c) Figure 3.2: Examples of load divisibility: (a) Indivisible load (b) Modularly divisible load (c) Modularly divisible load where the second module is also arbitrarily divisible. Time Processor ID 1 2 3 4 (a) Time Processor ID 1 2 3 4 (b) Time Processor ID 1 2 3 4 (c) 1 1 2 3 1 2 3 (a) (b) (c) Figure 3.3: Interaction graphs of the tasks represented in Figure 3.2. §3.2 Approach to the load balancing problem 123 1 2 3 4 5 6 1.1 1.2 1.3 1.4 1.5 1.6 2.1 2.2 2.3 2.4 2.5 2.6 3.1 3.2 3.3 3.4 3.5 3.6 Time step: t 1 t 2 t 3 Figure 3.4: Example of the task interaction graph of a mesh with 6 cells simulated in 6 processors during 3 time steps. 124 §3.2 Approach to the load balancing problem Processing a mesh is an arbitrarily divisible job whose flexibility is constrained by the solver. A rigid job would be that which can only be run using a fixed number of processors, 6 for instance. This is comfortable for the designer of the solver in charge of processing the mesh, because it is enough to partition the domain only once and this can be done before beginning any simulation. Nonetheless, the task scheduler of the optimizer would not be able to modify the amount of resources assigned to the task and a degradation of the parallel efficiency could happen when evaluating a batch of individuals. If a moldable job was provided, the task scheduler would be able to decide the amount of processors assigned to each task at the beginning of the evaluation of the batch of individuals. The decision could also be changed at any time before the execution of the task started. This would involve the task to have the ability to call a mesh partitioning module in order to fit the mesh topology to the number of processors assigned dynamically by the task scheduler. A malleable job owns the preferred and most complex degree of flexibility, thanks to the ability of being expanded or shrunk by the task manager in order to adapt to the available resources in run time. This degree of flexibility maximizes the overall parallel efficiency, but also requires a more complex implementation. Finally, an evolving job would be similar to the malleable one with the only difference that the job itself would be able to ask for resources to the task scheduler. 3.2.2 Factors affecting parallel performance The importance of maintaining a balanced load among the processors in order to achieve high parallel performance in large-scale multiprocessor systems has been already mentioned. The reason is that the total execution time of a set of tasks is determined by the last processor to finish its assigned work. The main factors affecting the parallel performance of applications, i.e. the crucial factors when designing a task schedule, have been identified and are listed hereafter. Factors related to the hardware: •Hardware heterogeneity: The same task has different completion times for different resources when considering a heterogeneous system [2]. Thus, information on the processors’ speed is to be taken into account by the scheduler in order to obtain a good balance of the computational load. A potential solution, similarly to that proposed in [8], consists on testing the hardware before the execution of the real set of tasks is started so that the scheduler owns benchmark results of all participating processor nodes indicating their performance. §3.2 Approach to the load balancing problem 125 •Network topology: A supercomputer is composed by several multiprocessor nodes connected to each other creating a network. Each node’s processors have a common memory space called shared-memory, and all the interconnected sharedmemory islands give rise to a distributed-memory network. On one hand, each processor is able to access almost instantaneously the shared-memory owned by its node. On the other hand, communication delays are considerably higher in a distributed-memory environment due to the fact that communications are carried out by message passing [3]. Indeed, communication on the network is becoming the bottleneck for the scaling of parallel applications [6]. Two conclusions may be extracted from these considerations. First, tasks using OpenMP, MPI or hybrid models (OpenMP+MPI) have to be carefully placed upon the machines so that hardware affinities are efficiently handled for optimal performance (processor-toprocessor affinity). Second, if the data needed by a task are distributed among the local memories of a distributed-memory network, that parallel task will have an affinity for a subset of the processors based on the locality of its memory references (task-to-processor affinity) [9]. Hence, the task scheduler should provide topology aware task placement techniques based on the mentioned processor-to-processor and task-to-processor affinities in order to improve the parallel performance. This feature is even more important in highly parallel applications like mesh-based CFD & HT simulations. •Amount of system’s memory: As it is mentioned in [9], the amount of memory available in the system can have an indirect impact on the total execution time by allowing the scheduler to modify and expand task-to-processor affinities and to reduce memory conflicts. For example, in a distributed-memory machine, heavilyshared data objects can be replicated in the memories of several processors. This replication expands the number of processors for which a particular task may have an affinity. Disadvantages of this replication are the cost of the additional memory, and perhaps more importantly, the time required to maintain coherence. •Processor partitioning model: In the “dedicated master” partitioning, a processor can be dedicated exclusively to scheduling operations. In the “peer partition” approach, the loss of a processor to scheduling is avoided and any processor can run the scheduler when necessary. Obviously, the dedicated master approach seems to be the most suitable model for performing sophisticated scheduling techniques since it exploits the full processing power of the designated processor for this aim [9]. 126 §3.2 Approach to the load balancing problem Factors related to the nature of tasks: •Load divisibility: A divisible computational load may be split into several tasks, what allows reaching a better load balance but inserts a scheduling overhead caused by the increase in the overall number of tasks. The ability of the scheduler for stopping and restarting a task is also useful as a load division technique. Nonetheless, [10] warns about the overhead of starting a task, which may be due to i) the time to transfer application input/output data to/from each computing resource, and ii) the potential latencies involved when initiating a computation or a communication. •Flexibility of jobs: As it has been previously explained, the internal design of arbitrarily divisible parallel jobs can ease or constrain their partitioning into smaller jobs, and this is the concept known as flexibility. The higher the flexibility of the jobs, the easier will be for the task scheduler to balance the overall computational load. Moreover, the use of malleable or evolving jobs allows implementing sophisticated dynamic scheduling algorithms able to adapt the load balance in run time. •Interdependencies between tasks: The precedence between tasks is another factor to be taken into account by the scheduler. Although tasks in the batches handled by optimization algorithms are independent and do not need any synchronization or inter-task communications, the use of load division techniques may generate precedence relationships between some of the newly generated tasks. •Heterogeneity in task run times depending on its input variables: The tasks contained in each batch handled by an optimization algorithm are of the same type, but their processing time may vary depending on the input variables (the optimization variables). Thus, the task scheduler should be aware of the processing time required by each job. [8] proposes to insert the parallel application into a class having a standardized interface (an enhancement of the FMI interface [11], for instance) that contains information of its performance and parallelizability and with which the task scheduler should be able to interact. Factors related to the scheduler’s configuration: •Limited scheduling time: As mentioned in [9], the time required to schedule the tasks can be in direct conflict with the desire to balance the computational §3.2 Approach to the load balancing problem 127 load. For instance, using small task sizes can lead to good load balancing since these small tasks can be used to fill in processor idle times. However, if some of the scheduling operations are performed at run time, the time required to perform the scheduling will be added directly to the execution time of each task. Since this scheduling time is generally independent of the execution time of the parallel tasks, the use of small task sizes can produce very high scheduling overheads, which can negate the advantages of using small tasks for load balancing. Hence, it is useful to limit the scheduling time in such a way that equilibrium between the computational load balance and the scheduling overhead is achieved. •Distributed scheduling: The scheduling operation can also be parallelized and distributed among several processors so as to get better task-to-processor maps by assigning greater computational power to the scheduler. However, it must be taken into account that the information required by the scheduler should be available in the memory of the processors in charge of running it, what may involve additional communication overhead [9]. 3.2.3 Applications using task schedulers The need of task scheduling algorithms is shared by several fields in the information technology. Operation of High Performance Computing (HPC) platforms may be one of the best known applications, but other emerging fields like cloud computing and grid computing have shown similar necessities. There is vast literature analyzing the algorithmic and last advances of Resource and Job Management Systems (RJMS), the software controlling the behavior of HPC platforms. Such system managers take into consideration several features for scheduling [6], like hierarchical resources, topology aware placement, energy consumption efficiency, quality of services, fairness between users or projects, etc. An important metric to evaluate the work of a RJMS on a platform is the observed system utilization. However, studies and logs of production platforms show that HPC systems generally suffer of significant underutilization rates. According to [6], less than 65% of the overall capacity of clusters is utilized throughout the year. Task scheduling is also a critical problem in cloud computing [2]. A cloud system could have plenty of users that may come from all over the world. Therefore, large-scale task scheduling happens frequently, seeking the cloud providers to reduce the total completion time of the tasks sent by the users. Grid computing shares computational power and data storage capacity over the 128 §3.2 Approach to the load balancing problem Internet and the goal of grid task schedulers is to achieve high system throughput and to match the applications’ needs with the available computing resources [12]. Scheduling on a grid is characterized by three main phases. Phase one is resource discovery, which generates a list of potential resources. Phase two involves gathering information about those resources and choosing the best set to match the applications’ requirements. In phase three the task is executed, which includes file staging and cleanup. Being the characteristics of these three applications (HPC platforms, cloud computing and grid computing) similar to those owned by the process of scheduling the evaluation of individuals from a batch created by an optimization algorithm, literature related to the mentioned fields has been considered a good starting point for knowing the state of the art of scheduling algorithms. 3.2.4 Methodology for developing load balancing algorithms An overview of the methodology that has been followed for developing and testing load balancing algorithms is explained in this section and can be summarized in four main points: 1. Some assumptions have been made regarding the parallel performance factors considered by the task scheduler in order to simplify the scheduling problem. Complexity will be increased by adding more factors gradually. 2. The load balancing problem has been split in 2 subproblems: tasks’ time estimation problem and task scheduling problem. 3. Metrics to evaluate the goodness of the results provided by the load balancing algorithm have been defined. 4. The way of modeling the execution of real tasks in a real High Performance Computing (HPC) system has been decided with the aim of ensuring the reproducibility of experiments as well as reducing their computational cost. Assumptions on the parallel performance factors Some assumptions on the parallel performance factors have been made with the aim of simplifying the scheduling problem. As scheduling algorithms that achieve satisfactory results are provided, the maximal complexity of the scheduling problem can be reached by modifying these assumptions gradually. §3.2 Approach to the load balancing problem 129 Assumptions related to the hardware: • Hardware heterogeneity: The hardware is considered to be homogeneous and scheduling tests are carried out in a single multiprocessor node composed by 32 processors. Thus, it is not necessary to provide any benchmark case for testing the performance of hardware components. • Network topology: Using a single multiprocessor node implies that all processorto-processor and task-to-processor affinities are equal. Since all processors have access to a common memory space (shared-memory), communications are quasiinstantaneous between any two of them without the need of message passing. • Amount of system’s memory: The memory of the utilized multiprocessor node is the only system’s memory available and may be accessed by any processor. Hence, it is not possible to increase any task-to-processor affinity by adding memory management algorithms. • Processor partitioning model: The “peer-partition” approach is used. Since a unique node of 32 processors will be available and initial scheduling algorithms are not expected to reach any considerable level of sophistication, it has been thought that sacrificing 1 processor for full-time scheduling is a greater waste of effective computational time than slowing down slightly the completion of some tasks in case a scheduler is run before their computation is started. Assumptions related to the nature of tasks: • Load divisibility: Every task is considered to be arbitrarily divisible by the scheduler, which can assign to the task a number of processors ranging between a minimum and a maximum defined by the user. The possibility of stopping and later restarting a task is not allowed, because the restart feature is not that common and/or it may involve a noticeable computational time overhead. Hence, the only way of dividing a task that has been considered is its parallel execution. • Flexibility of jobs: Every task is considered to be moldable, i.e. the scheduler assigns a number of processors to the task before its execution starts and this number remains constant until the completion of the task. • Interdependencies between tasks: All tasks belonging to a batch of individuals created by the optimizer are independent from each other, existing no precedence 130 §3.2 Approach to the load balancing problem relations among them. Since it has been assumed that tasks cannot be stopped and restarted, precedence relations will be neither generated in run time. • Heterogeneity in task run times depending on its input variables: It is possible to have tasks whose processing times vary according to the values of the input variables. Assumptions related to the scheduler’s configuration: • Limited scheduling time: It is necessary to limit the scheduling time in order to avoid an overall performance degradation caused by excessively long scheduling periods. • Distributed scheduling: The implementation of sequential task schedulers is foreseen. When satisfactory results are obtained, the possibility of parallelizing those schedulers will be considered. Splitting the load balancing problem in 2 subproblems Two steps are clearly differentiated in the resolution of the load balancing problem of a batch of individuals created by an optimizer and solved in a High Performance Computing (HPC) system: the estimation of the tasks’ evaluation times and the resolution of a combinatorial problem. The estimation of the tasks’ evaluation times consists on obtaining guess values of the amount of time required by each task if it is simulated in any subset of the available processors. As a result, a map of times is obtained for each task. The resolution of the combinatorial problem consists on designing a task schedule such that the makespan of the batch of individuals is minimized. The task scheduler uses the time maps calculated for every task in the previous step to estimate the overall simulation time of successive schedules and to try to optimize them by deciding where (in which processors) and when each task is to be processed. Definition of metrics for evaluating the goodness of a task schedule The definition of adequate metrics is important in order to evaluate the goodness of a method. In the case of Resource and Job Management Systems (RJMS) used in computing clusters, a typical metric to evaluate the performance of the scheduling algorithm is the observed system utilization. Nevertheless, there is a major difference between the queue of tasks managed by a RJMS and that managed by the task scheduler of an optimizer: the latter has a finite number of similar tasks, whereas the RJMS receives a §3.2 Approach to the load balancing problem 131 continuous flow of heterogeneous tasks coming from different users. Consequently, the makespan of the batch of tasks (i.e. the overall execution time taking place between the start time of the first processed task and the end time of the last task) is found to be a more representative metric than the observed system utilization in order to measure the quality of the obtained schedules for the batches of tasks to be evaluated by the optimizer. Modeling real tasks in an HPC system During the development of task schedulers, there are several reasons to avoid testing them by processing real tasks in a production HPC platform: considerable computational load of real tasks, long queues due to tasks sent by other users, variable environmental conditions in the HPC platform affecting negatively the reproducibility of the tests, the fact that specific experiments may need specialized software which is often hard to install on production platforms, etc. Thus, it is advisable to model the tasks to be processed as well as the aspects of the HPC platform to be taken into consideration. On one hand, and according to [6], tasks can be modeled with up to 3 fidelity degrees depending on the specific features under evaluation: •Sleep applications: They consist on using the Unix command “sleep” followed by a number which defines the amount of time in seconds that a processor will stay idle. The main advantage of sleep jobs is that they represent the simplest type of application with a predefined steady duration that can be performed in any number of processors. However, this kind of jobs is not influenced by CPU, bandwidth, or memory stress. •Synthetic applications: They are defined based on profiles of real applications and are commonly used with the goal of testing specific parts of the system like CPU, memory, network or I/O. They implicate real computation but do not capture the whole complexity of real applications. Some typical synthetic applications are the following: NAS NPB3.3 benchmarks, Linpack-HPL, pchksum benchmark, Mandelbrot set, etc. •Real applications: They hold real-life’s complexity and are therefore the most representative way of testing the designed algorithms. 234 References [6] Y. Georgiou. Contributions for resource and job management in high performance computing. PhD thesis, Université de Grenoble, France, 2010. [7] D. M. Valdivieso. Towards a virtual platform for aerodynamic design, performance assessment and optimization of Horizontal Axis Wind Turbines. PhD thesis, Universitat Politècnica de Catalunya, Spain, 2017. [8] P. Nordin, R. Braun, and P. Krus. Job-scheduling of distributed simulation-based optimization with support for multi-level parallelism. In Proceedings of the 56th Conference on Simulation and Modelling (SIMS 56), October, 7-9, 2015, Linköping University, Sweden, number 119, pages 187–197. Linköping University Electronic Press, 2015. [9] B. Hamidzadeh, L. Y. Kit, and D. J. Lilja. Dynamic task scheduling using online optimization. IEEE Trans. Parallel Distrib. Syst., 11(11):1151–1163, November 2000. [10] Y. Yang, K. van der Raadt, and H. Casanova. Multiround algorithms for scheduling divisible loads. IEEE Transactions on Parallel and Distributed Systems, 16(11):1092–1102, Nov 2005. [11] T. Blochwitz, M. Otter, J. Akesson, M. Arnold, C. Clauss, H. Elmqvist, M. Friedrich, A. Junghanns, J. Mauss, D. Neumerkel, et al. Functional mockup interface 2.0: The standard for tool independent exchange of simulation models. In Proceedings of the 9th International MODELICA Conference; September 3-5; 2012; Munich; Germany, number 076, pages 173–184. Linköping University Electronic Press, 2012. [12] W. Zhou and Y. P. Bu. An adaptive genetic algorithm for the grid scheduling problem. In 2012 24th Chinese Control and Decision Conference (CCDC), pages 730–734. IEEE, 2012. [13] E. Y. H Lin. A bibliographical survey on some well-known non-standard knapsack problems. INFOR: Information Systems and Operational Research, 36(4):274–317, 1998. [14] Jr. E. G. Coffman, M. R. Garey, and D. S. Johnson. Approximation algorithms for NP-hard problems. pages 46–93. PWS Publishing Co., Boston, MA, USA, 1997. [15] A. Lodi, S. Martello, and D. Vigo. Recent advances on two-dimensional bin packing problems. Discrete Appl. Math., 123(1-3):379–396, November 2002. References 235 [16] H. Izakian, A. Abraham, and V. Snasel. Comparison of heuristics for scheduling independent tasks on heterogeneous distributed environments. In Proceedings of the 2009 International Joint Conference on Computational Sciences and Optimization - Volume 01, CSO ’09, pages 8–12, Washington, DC, USA, 2009. IEEE Computer Society. [17] A. Ecer, Y. P. Chien, H. U. Akay, and J. D. Chen. Load balancing for multiple parallel jobs. In European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS; September 11-14; 2000, Barcelona, Spain, 2000. [18] F. A. Omara and M. M. Arafa. Genetic algorithms for task scheduling problem. Journal of Parallel and Distributed Computing, 70(1):13–22, 2010. [19] M. Tripathy and C. R. Tripathy. Dynamic load balancing with “work stealing” for distributed shared memory clusters. In 2010 International Conference on Industrial Electronics, Control & Robotics (IECR), pages 43–47. IEEE, 2010. [20] M. M. Najafabadi, M. Zali, S. Taheri, and F. Taghiyareh. Static task scheduling using genetic algorithm and reinforcement learning. In 2007 IEEE Symposium on Computational Intelligence in Scheduling, pages 226–230. IEEE, 2007. [21] C. Franke, J. Lepping, and U. Schwiegelshohn. Greedy scheduling with complex obejectives. In 2007 IEEE Symposium on Computational Intelligence in Scheduling, pages 113–120. IEEE, 2007. [22] A. Nissimov and D. G. Feitelson. Probabilistic backfilling. In E. Frachtenberg and U. Schwiegelshohn, editors, Job Scheduling Strategies for Parallel Processing: 13th International Workshop, JSSPP 2007, Seattle, WA, USA, June 17, 2007. Revised Papers, pages 102–115, Berlin, Heidelberg, 2008. Springer. [23] B. M. Adams, L. E. Bauman, W. J. Bohnhoff, K. R. Dalbey, M. S. Ebeida, J. P. Eddy, M. S. Eldred, P. D. Hough, K. T. Hu, J. D. Jakeman, J. A. Stephens, L. P. Swiler, D. M. Vigil, and T. M. Wildey. Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis: version 6.0 user’s manual. Sandia Technical Report SAND2014-4633, July 2014. Updated November 2015 (Version 6.3). [24] R. H. Myers, D. C. Montgomery, and C. M. Anderson-Cook. Response surface methodology: product and process optimization using designed experiments, 2009. 236 References [25] A. A. Giunta and L. T. Watson. A comparison of approximation modeling techniques: Polynomial versus interpolating models. In 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, number AIAA-984758, St. Louis, MO, 1998, pages 392–404. American Institute of Aeronautics and Astronautics, 1998. [26] D. C. Zimmerman. Genetic algorithms for navigating expensive and complex design spaces. Final Report for Sandia National Laboratories contract AO-7736 CA, 2:30332–0150, 1996. [27] J. H. Friedman. Multivariate adaptive regression splines. The Annals of Statistics, 19(1):1–141, 1991. [28] M. J. L. Orr. Introduction to radial basis function networks, 1996. [29] A. Nealen. An as-short-as-possible introduction to the least squares, weighted least squares and moving least squares methods for scattered data approximation and interpolation. Technical report, Discrete Geometric Modeling Group, Technishe Universitaet, Berlin, Germany, 2004, URL: http://www.nealen.com/projects. [30] V. E. Bazterra, M. Cuma, M. B. Ferraro, and J. C. Facelli. A general framework to understand parallel performance in heterogeneous clusters: analysis of a new adaptive parallel genetic algorithm. Journal of Parallel and Distributed Computing, 65(1):48–57, 2005. 4 Conclusions and future work 237 238 §4.1 Conclusions 4.1 Conclusions The objective of this Doctoral Thesis was the development of parallel optimization algorithms to be deployed in massively parallel computers. A generic mathematical optimization tool applicable in any field of science and engineering has been developed for this aim, although a special focus has been put on the application of the library to the fields of expertise of the institution hosting this research activity, i.e. the Heat and Mass Transfer Technological Center (CTTC). The main contribution has been the development and implementation of several load balancing strategies which have demonstrated to be able to reduce the simulation time required by the optimization tool. This was proved by an exhaustive theoretical case study. Finally, the time reduction attained in a real-world engineering application that consists on optimizing the refrigeration system of a power electronic device was presented as an illustrative example. A summary of the conclusions extracted in the previous chapters is provided in the subsequent paragraphs. A thorough state of the art study was conducted in the first chapter with the aim of obtaining a global scope of the mathematical optimization techniques available to date. Special emphasis was laid on the genetic algorithm. Then a description of the main concepts and shortcomings of the standard parallelization strategies available for such population-based optimization methods was included. Since the computational cost of a real-world design problem may be considerable, the development of appropriate parallelization techniques is a critical issue in order to benefit from optimization strategies. The goal of every parallelization strategy is to maximize the CPU usage by minimizing the communication overhead and by balancing the computational load correctly, thus avoiding idleness of allocated processors. In the case of traditional optimization algorithms, and particularly of genetic algorithms, the principal causes of computational load imbalance are the following: • Inappropriate ratio (no. individuals / no. processor groups): Usually, when several individuals may be evaluated simultaneously, each individual is assigned to a group of processors. In standard algorithms, all groups are formed by the same number of processors and remain unchanged during optimization. However, if the ratio (no. individuals / no. processor groups) is not well set, some processors may remain idle while others are immersed in computation. §4.1 Conclusions 239 • Heterogeneous parallel computer systems: Hardware heterogeneity translates into non-homogeneous objective evaluation times and the consequent loss of parallel performance of optimization algorithms. • Heterogeneous objective function evaluation time: Heterogeneity in the computational cost of evaluating the objective functions may cause an important load imbalance provided that the simulation time of the genetic algorithm is dominated by the evaluation time of the objective functions. This phenomenon happens when the evaluation cost of individuals is dependent on the optimization variables and is not unusual in heat transfer and nonlinear mechanics applications. The new optimization library named Optimus was implemented in the second chapter. After carrying out a state of the art study on available optimization libraries, it was decided to base the development of Optimus in two open-source packages: Paradiseo and Trilinos/Moocho. Paradiseo was chosen because of its implementation of genetic algorithms, whereas the availability of gradient-based local search methods was the appealing feature of Trilinos/Moocho. Optimus makes use of the best performing characteristics of both libraries, but additional functionality was also added. Finally, validation tests of the new library were carried out at the end of the chapter. They include the optimization of benchmark mathematical functions and two specific tests from the field of Computational Fluid Dynamics and Heat Transfer (CFD & HT), namely the optimization of the energy efficiency of a fridge and the optimization of the geometry of a pipe. Optimus succeeded in every test, proving the suitability of the library for solving real-world optimization problems. In the third chapter load balancing methods for parallel optimization algorithms were developed and studied. The approach to the load balancing problem was started by introducing the divisible load theory and the concept of job flexibility. After that, the main factors affecting the parallel performance of optimization algorithms were identified and divided into three groups: factors related to the hardware, factors related to the nature of tasks and factors related to the scheduler’s configuration. The next step consisted on defining the methodology for developing and testing new load balancing algorithms, and it was decided to split the load balancing problem into 2 subproblems: a tasks’ time estimation problem and a task scheduling combinatorial problem. Regarding the task scheduling problem, a state of the art study was carried out including the review of classical combinatorial problems and of other applications using similar scheduling algorithms, such as the Resource and Job Management Systems (RJMS) installed in High Performance Computing (HPC) platforms. A second state of 240 §4.1 Conclusions the art study was carried out in order to tackle the subproblem consisting on estimating the tasks’ evaluation time. Several data fitting methods were analyzed, as well as some suggestions made by other authors. Two different task managers were then proposed, namely the static task manager and the dynamic task manager, and three theoretical case studies were carried out in order to evaluate their parallel performance with respect to the traditional selfscheduling task manager. A first approach was obtained by carrying out simulations of short sleep jobs with linear scalability. The next step was the simulation of long sleep jobs with linear scalability, as their average computational time is more representative of the real target application of the developed load balancing algorithms, i.e. the optimization of CFD & HT models. Finally, long sleep jobs with non-linear scalability were simulated by modeling the degradation of parallel performance taking place as they are assigned an increasing number of processors. The starting point of the tests consisted on proving the existence of the loss of parallel efficiency caused by the computational load imbalance when the self-scheduling task manager is used. This is unavoidable in case of having tasks with heterogeneous computational times, even if the (no. individuals / no. processor groups) ratio is properly selected by the user. Indeed, the main problem when configuring the self-scheduling task manager consists on the impossibility of determining the optimal number of processors that should be assigned per task in order to minimize the makespan of the whole batch of tasks. The static and dynamic task managers clearly proved their ability to successfully balance the computational load of batches of tasks with heterogeneous completion times, being the main drawbacks the scheduling overhead and handling the task time estimation errors. Moreover, the tests confirmed the superiority of the dynamic task manager subject to an appropriate configuration and to the availability of accurate enough tasks’ time estimations. Some observations can be made regarding the scheduling overhead. On one hand, the longer the scheduling time is, the better task distributions may be obtained. On the other hand, the greater the average completion time of the simulated tasks, the lower the relative scheduling overhead with respect to the average makespan of the batch of tasks. Hence, the static and dynamic task managers maximize their performance the greater the average computational time of the tasks is. Regarding the task time estimation errors, their presence is unfortunately unavoidable and has a negative impact on the load balance achieved by the task managers. An acceptable upper bound for estimation errors was provided according to the obtained results. §4.2 Future work 241 Note as well that an optimization may comprise the execution of numerous batches of tasks. Consequently, the benefits of reducing the makespan of a batch are multiplied and a significant global impact may be achieved. Moreover, the computational time savings are maximized as an increasing number of processors is used for the optimization. Such a behavior makes the dynamic task manager highly advisable for exascale computing applications. Although the tasks’ time estimation subproblem was not properly tackled due to lack of time, 3 global data fitting methods were implemented for this aim: kriging interpolation, radial basis functions (RBF) and artificial neural networks (ANN). Moreover, some preliminary tests were performed using the Rosenbrock’s function in order to evaluate their accuracy and computational cost. Finally, an illustrative example consisting on the optimization of a power electronic device was presented in order to extrapolate the conclusions obtained in the theoretical case study to engineering cases involving real data processing. Despite the apparent simplicity of the studied thermal model, the load balancing strategies developed in this chapter outperformed the results obtained by the self-scheduling task manager. 4.2 Future work The inherent parallel nature of evolutionary computation is a promising factor in developing robust and scalable optimization algorithms which can help reach the goal of a successful integration of Optimization and High Performance Computing capabilities. Therefore, parallel evolutionary computation is expected to be an active research area in the near future. The contribution of this Doctoral Thesis has consisted on the proposal of some load balancing algorithms that can be integrated in evolutionary computation software packages to improve their parallel efficiency. Some suggestions to further develop the proposed algorithms are included hereafter. The tasks’ time estimation subproblem is a key aspect of the load balancing algorithms, but it was not studied in detail due to lack of time. Thus, it is crucial to carry out research in this direction so that the proposed load balancing algorithms can be comfortably used for solving production optimization problems. The research should include two aspects: i) selection or development of appropriate data fitting techniques, and ii) definition of a robust sampling method. In this sense, some criterion able to measure the accuracy of the subsequently created tasks’ time maps and to activate the load balancing algorithms provided that a certain accuracy level has been reached is to be implemented. The author’s experience shows that it is preferable to use the 242 §4.2 Future work self-scheduling task manager unless a minimal accuracy level has been reached. The reason behind is that too great time estimation errors are prone to create severe load imbalance, causing the inefficient use of computational resources and serious simulation delays. Note that the parallel execution of data fitting methods may be necessary as the number of points used for modeling time maps increases. Regarding the task scheduling combinatorial subproblem solved by both the static and dynamic task managers, a rigorous study of the task scheduling algorithm is needed in order to improve its robustness and convergence speed. The performance of different genetic operators could be compared, for instance, or even the use of an optimization method other than the genetic algorithm could be considered in case it is better suited to solve such combinatorial problems. As a result of this study, the scheduling algorithm may be able to achieve greater generations’ makespan reductions at a lower cost. The static and dynamic versions of the task management algorithm in which the aforementioned scheduler is contained have also some improvement possibilities and are outlined in the next three paragraphs. The first aspect to be improved is the method used for calculating the time limitation for the task scheduling algorithm. According to the author’s mind, a good approach could be the simultaneous use of 2 criteria, being the scheduling time limited in each case by the most restrictive criterion: i) a limitation depending on the previous generation’s makespan (as it was implemented for this Doctoral Thesis), and ii) a limitation depending on the size of the combinatorial problem to be solved, as it is thought that having the combinatorial optimization problems a finite solutions space, the time required by a certain search method to find a good solution could be standardized for different problem sizes. Such a combined criterion may also be useful for balancing hybrid optimizations, in which the size of the batches of individuals created by the global and the local methods may differ significantly. Moreover, it could also be used when the scheduler is run several times per generation by the dynamic task manager, as the number of individuals to be scheduled is progressively decreased. A second aspect to be improved is the CPU usage during the scheduling process. In the proposed algorithms the scheduler is always run in a single processor while other processors remain idle waiting for the scheduler to finish. Some alternatives to avoid wasting computational time in this manner could be to either parallelize the execution of the scheduler or to run several schedulers simultaneously and use the best among the obtained task distributions. In case a “dedicated master” processor partitioning model was used, the scheduler could be uninterruptedly run in the master processor with the most updated information and provide the best available schedule when required by §4.2 Future work 243 any other processor. A third aspect would consist on finding an automated method able to set the optimal configuration of the task management algorithm regarding the number of schedulers to be run per generation and the percentage defining the maximal time limitation for each scheduler run based on the previous generation’s makespan. Finally, note that hardware homogeneity was assumed in every case studied in the scope of this Doctoral Thesis. Nevertheless, hardware heterogeneity cannot be avoided in massively parallel applications and it involves variable inter-processor communication delays and different data processing speeds. These factors need to be taken into consideration by the task management algorithm to optimally map tasks to processors, being the complexity of the optimization problem to be solved substantially increased.