PATENT DOCUMENT

Abstract:
The underlying invention generally relates to the field of Estimation of Distribution Algorithm, especially to optimization problems, including single-objective optimization and Multi-Objective Optimization. 
     The proposed method for optimization comprises six steps. In a first step it provides an initial population or a data set with a plurality of members respectively represented by parameter sets. Then one or a plurality of fitness functions are applied to evaluate the quality of the members of the population. In a third step offspring of the population is generated by means of a stochastic model using information from all members of the population. One or a plurality of fitness functions are applied to evaluate the quality of the offspring with respect to the underlying problem of the optimization. In a fifth step offspring is selected. Lastly the method goes back to the third step until the quality reaches a threshold value.

Full Description:
CROSS-REFEFRENCE TO RELATED APPLICATIONS 
   This application is related to U.S. application Ser. No. 11/042,991 filed on Jan. 24, 2005. 
   FIELD OF THE INVENTION 
   The underlying invention generally relates to the field of Estimation of Distribution Algorithms (EDA), especially for optimization problems, including single-objective optimization (SOO) and multi-objective optimization (MOO). 
   BACKGROUND OF THE INVENTION 
   A Hybrid Representation proposed in Okabe, T., Jin, Y. and Sendhoff, B. Evolutionary Multi-Objective Optimisation with a Hybrid Representation, In Proceedings of Congress on Evolutionary Computation (CEC-2003), pages 2262-2269, 2003, which is incorporated by reference in its entirety, is exploiting the different dynamics of Genetic Algorithm (GA) and Evolution Strategy (ES). This difference is caused by the different offspring distribution. The new algorithm shows superior performance compared to the state-of-the-art MOO methods. A natural for the extension of the Hybrid Representation algorithm is to try to model the optimal offspring distribution directly. This new algorithm generates a stochastic model from the offspring distribution in the Parameter Space (PS). Based on this stochastic model, offspring will be generated. Neither crossover nor mutation are used. This type of algorithms is often called Estimation of Distribution Algorithm (EDA).  FIG. 1  shows the difference between EAs (a) and EDAs (b) as it is known from conventional systems. Besides the way how they generate offspring, both of them are basically the same. 
   Genetic Algorithms (GAs) are well known to be powerful tools to obtain optimal solutions for complex optimization problems. As an extension of GAs, recently Estimation of Distribution Algorithms (EDAs) have received considerable attention, see Tsutsui, et al, cited below. The first authors to use the term “EDA” for estimating the distribution of promising points were Mühlenbein and Paass. In the state of the art, several EDAs have been proposed. Since the research for EDAs was triggered by the problem of disruption in GA, a lot of methods are suitable for discrete problems only. However, recently, several papers for the continuous problem have been also published. 
   According to a recent survey of EDAs, Pelikan, M., Goldberg, D. E. and Lobo, F. A Survey of Optimization by Building and Using Probabilistic Models. Technical Report 99018, University of Illinois, Urbana-Champaign Urbana, Ill. 61801, 1999 which is incorporated by reference herein in its entirety, the proposed methods can be classified into three classes: (1) no interaction, (2) pairwise interaction, and (3) multivariable interaction. 
   The first class, i.e. no interaction, does not consider any epistasis at all. Thus, each locus is treated independently. The second class, i.e., pairwise interaction, considers pairwise dependency only. Finally, multivariable interaction takes into account any type of dependency between variables. 
   The classification of the state of the art EDA methods can be seen in the following Table 1-a for discrete domains and in Table 1-b for continuous domains. All of the references cited in these tables are incorporated by reference herein in their entirety. 
   
     
       
             
           
             
             
             
           
         
             
               TABLE 1-a 
             
             
                 
             
             
               EDA methods for discrete domains. 
             
             
                 
             
           
           
             
                 
             
           
        
         
             
               1) No 
               Univariate Marginal 
               Mühlenbein, H. and Mahnig, T. Evolutionary 
             
             
               Interaction 
               Distribution 
               Algorithms: From Recombination to Search 
             
             
                 
               Algorithm (UMDA) 
               Distributions. In Theoretical Aspects of Evolutionary 
             
             
                 
                 
               Computing, pages 135-173, 2000, and 
             
             
                 
                 
               Mühlenbein, H. and Mahnig, T. Evolutionary Synthesis 
             
             
                 
                 
               of Bayesian Networks for Optimization. Advances in 
             
             
                 
                 
               Evolutionary Synthesis of Intelligent Agent, pages 429-455, 
             
             
                 
                 
               2001. 
             
             
                 
               Population Based 
               Baluja, S. Population-Based Incremental Learning: A 
             
             
                 
               Incremental Learning 
               Method for Integrating Genetic Search Based Function 
             
             
                 
               (PBIL) 
               Optimization and Competitive Learning. Technical 
             
             
                 
                 
               Report CMU-CS-94-163, Carnegie Mellon University, 
             
             
                 
                 
               1994. 
             
             
                 
               Compact Genetic 
               Harik, G. R., Lobo, F. G. and Goldberg, D. E. The 
             
             
                 
               Algorithm (cGA) 
               Compact Genetic Algorithm. Technical Report 97006, 
             
             
                 
                 
               University of Illinois, Urbana-Champaign Urbana, IL 
             
             
                 
                 
               61801, 1997, and 
             
             
                 
                 
               Harik, G. R., Lobo, F. G. and Goldberg, D. E. The 
             
             
                 
                 
               Compact Genetic Algorithm. In Proceedings of 
             
             
                 
                 
               Congress on Evolutionary Computation (CEC-1998), 
             
             
                 
                 
               pages 523-528, 1998. 
             
             
                 
               Hill-Climing with 
               Kvasnicka, V., Pelikan, M. and Pospichal, J. Hill 
             
             
                 
               Learning (HCwL) 
               Climbing with Learning (An Abstraction of Genetic 
             
             
                 
                 
               Algorithm). Neural Network World, 6: 773-796, 1996. 
             
             
                 
               Incremental 
               Mühlenbein, H. The Equation for the Response to 
             
             
                 
               Univariate Marginal 
               Selection and Its Use for Prediction. Evolutionary 
             
             
                 
               Distribution 
               Computation, 5(3): 303-346, 1998. 
             
             
                 
               Algorithm (IUMDA) 
             
             
                 
               Reinforcement 
               Paul, T. K. and Iba, H. Reinforcement Learning 
             
             
                 
               Learning Estimation 
               Estimation of Distribution Algorithm. In Proceedings of 
             
             
                 
               of Distribution 
               Genetic and Evolutionary Computation Conference 
             
             
                 
               Algorithm (RELEDA) 
               (GECCO-2003), pages 1259-1270, 2003. 
             
             
               2) Pairwise 
               Mutual Information 
               de Bonet, J. S., Isbell, J., Charles, L. and Viola, P. 
             
             
               Interaction 
               Maximization for 
               MIMIC: Finding Optima by Estimating Probability 
             
             
                 
               Input Clustering 
               Densities. Advances in Neural Information Processing 
             
             
                 
               (MIMIC) 
               Systems, 9: 424-431, 1996. 
             
             
                 
               Combining 
               Baluja, S. and Davies, S. Combining Multiple 
             
             
                 
               Optimizers with 
               Optimization Runs with Optimal Dependency Trees. 
             
             
                 
               Mutual Information 
               Technical Report CMU-CS-97-157, Carnegie Mellon 
             
             
                 
               Trees (COMIT) 
               University, 1997. 
             
             
                 
               Bivariate Marginal 
               Mühlenbein, H. and Mahnig, T. Evolutionary Synthesis 
             
             
                 
               Distribution 
               of Bayesian Networks for Optimization. Advances in 
             
             
                 
               Algorithm (BMDA) 
               Evolutionary Synthesis of Intelligent Agent, pages 429-455, 
             
             
                 
                 
               2001, and 
             
             
                 
                 
               Pelikan, M. and Mühlenbein, H. Marginal Distributions 
             
             
                 
                 
               in Evolutionary Algorithms. In Proceedings of the Forth 
             
             
                 
                 
               International Conference on Genetic Algorithms 
             
             
                 
                 
               (Mendel-1998), pages 90-95, 1998. 
             
             
               3) Multi- 
               Extended Compact 
               Harik, G. R. Linkage Learning via Probabilistic 
             
             
               variable 
               Genetic Algorithm 
               Modeling in the ECGA. Technical Report 99010, 
             
             
               Interaction 
               (ECGA) 
               University of Illinois, Urbana-Champaign Urbana, IL 
             
             
                 
                 
               61801, 1999, and 
             
             
                 
                 
               Lobo, F. G. and Harik, G. R. Extended Compact Genetic 
             
             
                 
                 
               Algorithm in C++. Technical Report 99016, University 
             
             
                 
                 
               of Illinois, Urbana-Champaign Urbana, IL 61801, 1999. 
             
             
                 
               Factorized 
               Mühlenbein, H. and Mahnig, T. FDA - A Scalable 
             
             
                 
               Distribution 
               Evolutionary Algorithm for the Optimization of 
             
             
                 
               Algorithm (FDA) 
               Additively Decomposed Functions. Evolutioanry 
             
             
                 
                 
               Computation, 7(1): 45-68, 1999, and 
             
             
                 
                 
               Mühlenbein, H. and Mahnig, T. The Factorized 
             
             
                 
                 
               Distribution Algorithm for Additively Decomposed 
             
             
                 
                 
               Functions. In Proceedings of Congress on Evolutionary 
             
             
                 
                 
               Computation (CEC-1999), pages 752-759, 1999, and 
             
             
                 
                 
               Mühlenbein, H. and Mahnig, T. Evolutionary Synthesis 
             
             
                 
                 
               of Bayesian Networks for Optimization. Advances in 
             
             
                 
                 
               Evolutionary Synthesis of Intelligent Agent, pages 429-455, 
             
             
                 
                 
               2001. 
             
             
                 
               Polytree 
               Larrañaga, P. and Lozano, J. A., editor. Estimation of 
             
             
                 
               Approximatin of 
               Distribution Algorithms. A New Tool for Evolutionary 
             
             
                 
               Distribution 
               Computation. Kluwer Academic Publishers, 2002. 
             
             
                 
               Algorithm (PADA) 
             
             
                 
               Estimation of 
               Larrañaga, P., Etxeberria, R., Lozano, J. A. and Pena, J. M. 
             
             
                 
               Bayesian Networks 
               Optimization by Learning and Simulation of 
             
             
                 
               Algorithm (EBNA) 
               Bayesian and Gaussian Networks. Technical Report 
             
             
                 
                 
               EHU-KZAA-IK-4/99, Department of Computer Science 
             
             
                 
                 
               and Artificial Intelligence, University of the Basque 
             
             
                 
                 
               Country, 1999. 
             
             
                 
               Bayesian 
               Khan, N., Goldberg, D. E. and Pelikan, M. Multi- 
             
             
                 
               Optimization 
               objective Bayesian Optimization Algorithm. Technical 
             
             
                 
               Algorithm (BOA) 
               Report 2002009, Univeristy of Illinois, Uvbana- 
             
             
                 
                 
               Champaign, Urbana, IL 61801, 2002, and 
             
             
                 
                 
               Laumanns, M. and Ocenasek, J. Bayesian Optimization 
             
             
                 
                 
               Algorithms for Multi-objective Optimization. In 
             
             
                 
                 
               Proceedings of Parellel Problem Solving from Nature 
             
             
                 
                 
               VII (PPSN-VII), pages 298-307, 2002, and 
             
             
                 
                 
               Pelikan, M., Goldberg, D. E. and Cantu-Pax, E. BOA: 
             
             
                 
                 
               The Bayesian Optimization Algorithm. In Proceedings 
             
             
                 
                 
               of Genetic and Evolutionary Computation Conference 
             
             
                 
                 
               (GECCO-1999), pages 525-532, 1999, and 
             
             
                 
                 
               Pelikan, M., Goldberg, D. E. and Cantu-Paz, E. Linkage 
             
             
                 
                 
               Problem, Distribution Estimation and Bayesian 
             
             
                 
                 
               Networks. Technical Report 98013, University of 
             
             
                 
                 
               Illinois, Urbana-Champaign Urbana, IL 61801, 1998, 
             
             
                 
                 
               and 
             
             
                 
                 
               Pelikan, M., Goldberg, D. E. and Cantu-Paz, E. BOA: 
             
             
                 
                 
               The Bayesian Optimization Algorithm. Technical 
             
             
                 
                 
               Report 99003, University of Illinois, Urbana-Champaign 
             
             
                 
                 
               Urbana, IL 61801, 1999, and 
             
             
                 
                 
               Zhang, B.-T. A Bayesian Framework for Evolutionary 
             
             
                 
                 
               Computation. In Proceedings of Congress on 
             
             
                 
                 
               Evolutionary Computation (CEC-1999), pages 722-728, 
             
             
                 
                 
               1999. 
             
             
                 
               Learning Factorized 
               Mühlenbein, H. and Mahnig, T. Evolutionary Synthesis 
             
             
                 
               Distribution 
               of Bayesian Networks for Optimization. Advances in 
             
             
                 
               Algorithm (LFDA) 
               Evolutionary Synthesis of Intelligent Agent, pages 429-455, 
             
             
                 
                 
               2001. 
             
             
                 
               Factorized 
               Larrañaga, P. and Lozano, J. A., editor. Estimation of 
             
             
                 
               Distribution 
               Distribution Algorithms. A New Tool for Evolutionary 
             
             
                 
               Algorithm (FDA) 
               Computation. Kluwer Academic Publishers, 2002. 
             
             
                 
               (extended version) 
             
             
                 
             
           
        
       
     
   
   
     
       
             
           
             
             
             
           
         
             
               TABLE 1-b 
             
             
                 
             
             
               EDA methods for continuous domains. 
             
             
                 
             
           
           
             
                 
             
           
        
         
             
               1) No 
               Univariate Marginal 
               Larrañaga, P., Etxeberria, R., Lozano, J. A. and Peña, 
             
             
               Interaction 
               Distribution 
               J. M. Optimization by Learning and Simulation of 
             
             
                 
               Algorithm (UMDAC) 
               Bayesian and Gaussian Networks. Technical Report 
             
             
                 
                 
               EHU-KZAA-IK-4/99, Department of Computer Science 
             
             
                 
                 
               and Artificial Intelligence, University of the Basque 
             
             
                 
                 
               Country, 1999. 
             
             
                 
               Stochastic Hill 
               Rudlof, S. and Köppen, M. Stochastic Hill Climbing 
             
             
                 
               Climbing with 
               with Learning by Vectors of Normal Distributions. In 
             
             
                 
               Learning by Vectors 
               Proceedings of the First Online Workshop on Soft 
             
             
                 
               of Normal 
               Computing (WSC1), Nagoya, Japan, 1997. 
             
             
                 
               Distributions 
             
             
                 
               (SHCLVND) 
             
             
                 
               Population Based 
               Sebag, M. and Ducoulombier, A. Extending Population- 
             
             
                 
               Incremental Learning 
               Based Incremental Learning to Continuous Search 
             
             
                 
               (PBILC) 
               Spaces. In Proceedings of Parallel Problem Solving 
             
             
                 
                 
               from Nature V (PPSN-V), pages 418-427, 1998. 
             
             
               2) Pairwise 
               Mutual Information 
               Larrañaga, P., Etxeberria, R., Lozano, J. A. and Peña, 
             
             
               Interaction 
               Maximization for 
               J. M. Optimization by Learning and Simulation of 
             
             
                 
               Input Clustering 
               Bayesian and Gaussian Networks. Technical Report 
             
             
                 
               (MIMICC) 
               EHU-KZAA-IK-4/99, Department of Computer Science 
             
             
                 
                 
               and Artificial Intelligence, University of the Basque 
             
             
                 
                 
               Country, 1999. 
             
             
               3) Multivariable 
               Estimation of 
               Larrañaga, P. and Lozano, J. A., editor. Estimation of 
             
             
               Interaction 
               Multivariate Normal 
               Distribution Algorithms. A New Tool for Evolutionary 
             
             
                 
               Algorithm (EMNA) 
               Computation. Kluwer Academic Publishers, 2002. 
             
             
                 
               Estimation of 
               Larrañaga, P. and Lozano, J. A., editor. Estimation of 
             
             
                 
               Gaussian Networks 
               Distribution Algorithms. A New Tool for Evolutionary 
             
             
                 
               Algorithm (EGNA) 
               Computation. Kluwer Academic Publishers, 2002. 
             
             
                 
               Iterated Density 
               Bosman, P. A. N. and Thierens, D. An Algorithmic 
             
             
                 
               Estimation Algorithm 
               Framework for Density Estimation Based Evolutionary 
             
             
                 
               (IDEA) 
               Algorithms. Technical Report UU-CS-1999-46, 
             
             
                 
                 
               Department of Computer Science, Utrecht University, 
             
             
                 
                 
               1999, and 
             
             
                 
                 
               Bosman, P. A. N. and Thierens, D. Continuous Iterated 
             
             
                 
                 
               Density Estimation Evolutionary Algorithms within the 
             
             
                 
                 
               IDEA Framework. Technical Report UU-CS-2000-15, 
             
             
                 
                 
               Department of Computer Science, Utrecht University, 
             
             
                 
                 
               2000, and 
             
             
                 
                 
               Bosman, P. A. N. and Thierens, D. IDEAs Based on the 
             
             
                 
                 
               Normal Kernels Probability Density Function. Technical 
             
             
                 
                 
               Report UU-CS-2000-11, Department of Computer 
             
             
                 
                 
               Science, Utrecht University, 2000, and 
             
             
                 
                 
               Bosman, P. A. N. and Thierens, D. Mixed IDEAs. 
             
             
                 
                 
               Technical Report UU-CS-2000-45, Department of 
             
             
                 
                 
               Computer Science, Utrecht University, 2000, and 
             
             
                 
                 
               Bosman, P. A. N. and Thierens, D. Negative Log- 
             
             
                 
                 
               Likelihood and Statistical Hypothesis Testing as the 
             
             
                 
                 
               Basis of Model Selection in IDEAs. In Proceedings of 
             
             
                 
                 
               the Tenth Belgium-Netherlands Conference on Machine 
             
             
                 
                 
               Learning, pages 109-116, 2000, and 
             
             
                 
                 
               Bosman, P. A. N. and Thierens, D. Advancing 
             
             
                 
                 
               Continuous IDEAs with Mixture Distributions and 
             
             
                 
                 
               Factorization Selection Metrics. In Proceedings of the 
             
             
                 
                 
               Optimization by Building and using Probabilistic 
             
             
                 
                 
               Models OBUPM Workshop at the Genetic and 
             
             
                 
                 
               Evolutionary Computation Conference (GECCO-2001), 
             
             
                 
                 
               pages 208-212, 2001, and 
             
             
                 
                 
               Thierens, D. and Bosmann, P. A. N. Multi-Objective 
             
             
                 
                 
               Mixture-based Iterated Density Estimation Evolutionary 
             
             
                 
                 
               Algorithms. In Proceedings of Genetic and Evolutionary 
             
             
                 
                 
               Computation Conference (GECCO-2001), pages 663-670, 
             
             
                 
                 
               2001. 
             
             
                 
               Parzen Estimation of 
               Costa, M. and Minisci, E. MOPED: A Multi-objective 
             
             
                 
               Distribution 
               Parzen-based Estimation of Distribution Algorithm for 
             
             
                 
               Algorithm (PEDA) 
               Continuous Problems. In Proceedings of the Second 
             
             
                 
                 
               International Conference on Evolutionary Multi- 
             
             
                 
                 
               Criterion Optimization (EMO-2003), pages 282-294, 
             
             
                 
                 
               2003. 
             
             
                 
               Marginal Histgram 
               Tsutsui, S., Pelikan, M. and Goldberg, D. E. 
             
             
                 
               Model (MHM) 
               Probabilistic Model-building Genetic Algorithms Using 
             
             
                 
                 
               Marginal Histograms in Continuous Domain. In 
             
             
                 
                 
               Proceedings of the KES-2001, Knowledge-based 
             
             
                 
                 
               Intelligent Information Engineering Systems and Allied 
             
             
                 
                 
               Technologies, volume 1, pages 112-121, and 
             
             
                 
                 
               Tsutsui, S., Pelikan, M. and Goldberg, D. E. 
             
             
                 
                 
               Evolutionary Algorithm Using Marginal Histogram 
             
             
                 
                 
               Models in Continuous Domain. Technical Report 
             
             
                 
                 
               2001019, University of Illinois, Urbana-Champaign 
             
             
                 
                 
               Urbana, IL 61801, 2001. 
             
             
                 
             
           
        
       
     
   
   Since the method proposed by the present invention belongs to the class of multivariable interaction, some of the popular methods in this category will be explained next. 
   To learn the linkage among parameters and therefore the structure of the problem, Bayesian networks are used. With the Bayesian Networks, the conditional probability is approximated. Each node and connection in the Bayesian networks correspond to the parameters and the conditional probability, respectively. Finally, the factorized probability is used to generate offspring. Recently, this method has been applied to MOO problems and has gathered much attention as per Khan, N., Goldberg, D. E. and Pelikan, M. Multi-objective Bayesian Optimization Algorithm. Technical Report 2002009, Univeristy of Illinois, Urbana-Champaign, Urbana, Ill. 61801, 2002 and Laumanns, M. and Ocenasek, J. Bayesian Optimization Algorithms for Multi-objective Optimization. In Proceedings of Parellel Problem Solving from Nature VII (PPSN-VII), pages 298-307, 2002 which are incorporated by reference herein in their entirety. Iterated Density Estimation Evolutionary Algorithm (IDEA) 
   Bosman and Thierens have proposed four types of EDAs that all belong to the class of IDEA. The first one is for the discrete domain where the conditional probability is used to build up the stochastic model. The others are for the continuous domain. A normalized Gaussian, a histogram method and a kernel method are used to generate the stochastic model. The kernel based method has been also applied to MOO, called Mixture-based IDEA (MIDEA). 
   Parzen-Based Estimation of Distribution Algorithm (PEDA) 
   To generate the stochastic model, a Parzen estimator is used to approximate the probability density of solutions. Based on the stochastic model, new offsprings will be generated. This method has been used for MOO problems. 
   Marginal Histogram Model (MHM) 
   For each parameter, the search space is divided into small bins. The ratio of the number of individuals in each bin to the whole number is assigned as the selection probability. With this probability, a bin is selected randomly. In the selected bin, an offspring is generated uniformly. 
   In view of these references, it is the object of the present invention to propose a more efficient algorithm for optimization. 
   SUMMARY OF THE INVENTION 
   A new algorithm, which can be called Voronoi-based Estimation of Distribution Algorithm (VEDA), is using a Voronoi mesh to represent the stochastic model. The Voronoi mesh is usually used in another context, i.e. finite element methods (FEMs) for generating the mesh. The mesh covers the parameter space. Each cell of the mesh is represented by its center. The rank of the center is determined by the fitness of the individuals belonging to its cell. 
   One of the characteristics of VEDA is to exploit all available information. In the canonical EAs and EDAs, only selected (promising) individuals contribute to the generation of offspring. However, even solutions of poor quality can serve as “bad examples”, therefore, they are not ignored but exploited by the technique according to the present invention. In particular, in the context of complex optimization problems with computationally expensive evaluations, it is of importance to use essentially all available information. In the context of EDAs poor solutions are helpful as they induce a low probability of the stochastic model in their neighborhood in the parameter space. 
   To generate offspring efficiently, clustering and Principal Component Analysis (PCA) are used in VEDA. In particular, in higher dimensional cases, it becomes more time consuming to generate the stochastic model. To reduce the computational cost, both techniques are used. 
   The basic idea behind the Voronoi-based Estimation of Distribution Algorithm (VEDA) stems from the observations on the intrinsic interplay between the distributions in the parameter space and in the fitness space. On the one hand, it is sensible to directly estimate the most appropriate search distribution for multi-objective optimization problems. Therefore, EDAs is a sensible starting point. On the other hand, both in canonical evolutionary algorithms as well as in EDAs, individuals which are not selected are completely discarded, see  FIG. 2 . The information that they carry, i.e. “do not search close to this solution again”, is wasted. In particular, in the context of distribution based algorithms, this does not seem sensible. If a human operator would choose the next generation, he/she would certainly take this “negative” information into account. Thus, the second target of the VEDA is to exploit not just the “good” individuals in the selection process but also low-performance individuals as “bad examples”. 
   One aspect of the invention relates to the offspring distribution in the Parameter Space (PS) with the information given in the Fitness Space (FS). Based on the offspring, a stochastic model is generated from which the next generation will be determined. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     Further advantages of the underlying invention result from the subordinate claims as well as from the following description of two preferred embodiments of the invention which are depicted in the following drawings. 
       FIG. 1  shows the difference between known Evolutionary Algorithms (EAs) and known Estimation of Distribution Algorithms (EDAs) by depicting the basic flow of EAs (a) and EDAs (b). 
       FIG. 2  shows the basic idea of Voronoi-based Estimation of Distribution Algorithm (VEDA) according to one aspect of the present invention, wherein the offspring distribution in the PS with the information given in the FS are used to generate the promising offspring. 
       FIG. 3  shows the basic flow of VEDA. 
       FIG. 4  shows the detailed flow in VEDA and particularly the flow of generating a stochastic model (a) and the flow of generating offspring (b). 
       FIG. 5  shows three clusters determined by k-means clustering whereby it is obvious that it is easier to identify stochastic models for each of the clusters separately instead of for the whole data set. 
       FIG. 6  shows sample data sets for data without epistasis (a) and data with epistasis (b). 
       FIG. 7  shows the results of Unimodal Normal Distribution Crossover (UNDX) on SCH1 (n=50). 
       FIG. 8  shows that the search power for the orthogonal direction will be lost, if all solutions are on one line. 
       FIG. 9  shows perturbation added in the orthogonal direction and extension in the principal direction, β and γ being the differences between the maximum value and the minimum value in the principal and the orthogonal direction, δ being the perturbation in the orthogonal direction, and α being the extension for the principal direction. 
       FIG. 10  shows a discrete Voronoi mesh. 
       FIG. 11  shows the sample flow of generation of the Voronoi mesh, wherein the thicker lines outline the boundary of the Voronoi mesh. 
       FIG. 12  shows a geometry distribution with PG=0.2, 0.5, 0.8. 
       FIG. 13  shows the generated stochastic model for OKA4 in the parameter space, the Pareto front being equal to x 2 =x 1 ±4√{square root over (x 1 )}+4. 
       FIG. 14  shows the results of VEDA and NSGA-II. 
       FIG. 15  shows the solutions in the PS on OKA4 by VEDA and NSGA-II. 
       FIG. 16  shows a comparison of VEDA concept according to the present invention with the conventional systems. 
       FIG. 17  shows a graphical output of VEDA, wherein the solid rectangle is the output of VEDA, and the dotted line is the true Pareto front in the PS. 
   

   DETAILED DESCRIPTION OF THE INVENTION 
   Note that all references to prior are documents in brackets are to be understood as explicitly referring to these documents and incorporating the disclosure thereof by reference. According to the invention, a new algorithm for tackling optimization problems, including Single-Objective Optimization (SOO) and Multi-Objective Optimization (MOO), is proposed, namely the Voronoi-based Estimation of Distribution Algorithm (VEDA). 
   Based on the offspring distribution in the PS with the information given in the FS, a stochastic model is generated using a Voronoi mesh. With the model, the offspring are generated stochastically. 
   The basic flow of VEDA is shown in  FIG. 3 . Since some issues of EDAs as described by Pelikan, M., Goldberg, D. E. and Lobo, F. A Survey of Optimization by Building and Using Probabilistic Models. Technical Report 99018, University of Illinois, Urbana-Champaign Urbana, Ill. 61801, 1999, which is incorporated by reference herein in its entirety include (1) how to make a stochastic model, and (2) how to generate new offspring with the stochastic model, the detailed flows for these steps are shown in  FIGS. 4(   a ) and ( b ), respectively. 
     FIG. 3  explains the overall flow of VEDA concept according to the present invention. First, the stochastic model will be generated from a database (S 3   a ). If the database does not exist yet, initial data will be generated (S 1   a ) randomly and evaluated (S 1   b ). Based on the generated model, new promising individuals will be generated (S 3   b ) and evaluated (S 3   c ). With the fast ranking method, see Deb, K., Pratap, A., Agarwal, S. and Meyarivan, T. A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2):182-197, 2002, which is incorporated by reference herein in its entirety, the rank of the individuals will be calculated (S 3   d ). Using the crowded tournament selection, individuals will be selected (S 3   e ) and stored (S 3   f ) in the database. The non-selected individuals will be stored in a different database as “bad examples”. By adding new data in the database, the rank information will become incorrect. Since the rank of the former data was calculated without the latest data, the rank of them should be updated (S 3   g ). If a given terminal condition is met, VEDA will stop, otherwise the same flow will be repeated (S 3   h ). In the database, design parameters, fitness values and rank are stored. 
     FIG. 4  explains the details of the approach to build up the model and to generate offspring. To generate the stochastic model, a clustering is used (S 4   a ), see Kaufman, L. and Rousseeuw, P. J. Finding Groups in Data—An Introduction to cluster Analysis. Wiley-Interscience, 1990 which is incorporated by reference herein in its entirety. In each cluster, Principal Component Analysis (PCA) is carried out (S 4   b ) to get the principal axis and the orthogonal axis of the distribution. See Hyvärinen, A. Independent Component Analysis. Wiley-Interscience, 2001 and Jolliffe, I. T. Principal Component Analysis. Springer Verlag, 2002 which are incorporated by reference herein in their entirety. The data will be projected (S 4   c ) to the new coordinate system determined by PCA. Some additional measures like adding a perturbation (S 4   d ) to the model have to be taken to avoid local minima, which is explained below in the section “Perturbation of Offspring Distribution”. The minimum and maximum value for each axis is then calculated (S 4   e ). Since PCA was carried out, epistasis among the design parameters should be minimal. In the new coordinate system, a Voronoi mesh is generated (S 4   f ) as the stochastic model. Based on the assigned rank, the probability will be calculated for each cell (S 4   g ). 
   To generate a new individual, a cell is selected (S 4   h ) based on the assigned probability. After choosing a cell, a new individual is generated (S 4   i ) in the selected cell uniformly. Finally, the new individual is projected (S 4   j ) to the real coordinate system. In a further embodiment of the present invention, a system for optimization problems, including Single-Objective Optimization (SOO) and Multi-Objective Optimization (MOO), comprising means for carrying out the steps of the presented optimization method is proposed. The use of this method can be for optimizing aerodynamic or hydrodynamic designs. 
   The following items distinguish the VEDA concept according to the present invention over a histogram method with a Voronoi mesh:
     1. The shape of a mesh in VEDA will change automatically. Since VEDA uses a Voronoi mesh, the shape of a mesh is not fixed.   2. The number of data in a cell is not counted. Although a histogram method needs a huge number of data, VEDA needs only few data.   3. VEDA uses a rank instead of frequency of data in a cell. Although a histogram method is based on frequency of data in a cell, VEDA does not use the frequency. Thus, VEDA can reduce the number of data to generate a stochastic model.   4. VEDA generates a different stochastic model from a histogram method. Since a histogram method is based on frequency of data in a cell, the cell without data has the lowest probability. However, the probability in a cell without data is approximated using the nearest data in VEDA. Thus, the probability in the cell without data is not always the lowest probability.   

   Since the way to make the model and generate offspring are issues in EDAs, each step is explained in more detail in the following sections. 
   Clustering 
   In the VEDA proposed in the present invention, data clusters are determined in the first step. The reason is that it is easier and more efficient to determine several local stochastic models, i.e. one for each cluster, instead of one global model. 
   Furthermore,  FIG. 5  illustrates that also fewer empty cells will occur when using a clustering method.  FIG. 5  shows three clusters determined by k-means clustering. It is obvious that it will be easier to identify stochastic models for each of the clusters separately instead of for the whole data set. 
   In the context of the present invention, e.g. k-means clustering proposed by MacQueen in Kaufman, L. and Rousseeuw, P. J. Finding Groups in Data—An Introduction to cluster Analysis. Wiley-Interscience, 1990 can be used, which is incorporated by reference herein in its entirety. The calculation of k-means clustering proceeds via the following steps:
     1. An initial partition of the data into k non-empty subsets is randomly generated.   2. Seed points as the centroids of the clusters of the current partition are computed.   3. Each datum is assigned to the cluster with the nearest seed point. If data change from one cluster to another, the seed points are recalculated.   4. If no change in the assignment of data has occurred, stop. If there has been changed, go to step 2.   

   As pointed out in Kaufman et al, the drawback of k-means clustering is the dependency on the order of the data. 
   To use the k-means clustering, one has to determine the value of k, i.e. the number of clusters. In this work, the value of k will be determined at random within a range of [1, 10]. Better performance can be reached when k is set using some a priori knowledge on the Pareto front in the parameter space, e.g., the number of cluster can be determined by the number of disconnected pieces of the Pareto front. This also indicates why k=1, i.e. no clustering, can be viable choice, e.g., for SCH1. 
   Principal Component Analysis 
   To reduce the dimensionality and to generate the model efficiently, Principal Component Analysis (PCA) is used.  FIG. 6  shows two different data sets. If there is epistasis between the variables, see  FIG. 6(   b ), it is reasonable to map them into a coordinate system minimizing the linear correlations, see  FIG. 6(   a ). Thereby, we can reduce the dimensionality and generate the model efficiently. 
   The PCA can be explained as follows (see Takahashi, M. and Kita, H. A Crossover Operator Using Independent Component Analysis for Real-coded Genetic Algorithms. In Proceedings of Congress on Evolutionary Computation (CEC-2001), pages 643-649, 2001 which is incorporated by reference herein in its entirety for additional details). Suppose there exist m sets of data of n variables, which are represented by an n×m matrix
 
 X={xij }( i =1 , . . . , n, j =1 , . . . , m ).
 
The variance-covariance matrix S={sij} of X is given by:
 
                   S   ij     =       1     m   -   1       ·       ∑     k   =   1     m     ⁢           ⁢       x   ik     ⁢     x   jk                   (     Eq   .           ⁢   1     )               
where data are normalized in advance so that the mean of each variable becomes zero. Transforming the data X with an n×m matrix A to Y=AX, the variance-covariance matrix
 
               S   ′     =     {       s   ′     ⁢   ij     }       ,         s   ′     ⁢   ij     =         1     m   -   1       ·       ∑     k   =   1     m     ⁢           ⁢       y   ik     ⁢     y   jk     ⁢           ⁢   of   ⁢           ⁢   Y         =     {   yij   }               
is given by:
 S′=ASA T    (Eq.2) 
   The target of the principal component analysis (PCA) is to find a matrix A that makes S′ the identity matrix. That is, to find a matrix that eliminates the correlation among variables of Y. Since S is a real symmetric matrix, there exists an orthogonal matrix P and a diagonal matrix Λ such that P T  SP=Λ. A diagonal matrix Λ whose diagonal elements are eigen values λ 1 , . . . ,λn of S, and a matrix P whose columns are corresponding eigenvectors are a pair of such matrices. Hence, one obtains:
 
 S ′=( AP )Λ( AP ) T    (Eq. 3)
 
   To make S′ the identity matrix, a transformation matrix A can be obtained as follows: 
   
     
       
         
           
             
               
                 A 
                 = 
                 
                   
                     diag 
                     ⁡ 
                     
                       ( 
                       
                         
                           1 
                           
                             
                               λ 
                               1 
                             
                           
                         
                         , 
                         
                           1 
                           
                             
                               λ 
                               2 
                             
                           
                         
                         , 
                         … 
                         ⁢ 
                         
                             
                         
                         , 
                         
                           1 
                           
                             
                               λ 
                               n 
                             
                           
                         
                       
                       ) 
                     
                   
                   ⁢ 
                   
                     P 
                     T 
                   
                 
               
             
             
               
                 ( 
                 
                   Eq 
                   . 
                   
                       
                   
                   ⁢ 
                   4 
                 
                 ) 
               
             
           
         
       
     
   
   Assuming all the eigenvalues are positive, we can obtain the inverse transformation A−1 as follows:
 
 A   −1   =P  diag(√{square root over (λ 1 )}, √{square root over (λ 2 )}, . . . , √{square root over (λ n )})  (Eq. 5)
 
   Perturbation of Offspring Distribution 
   In the research for Real-Coded Genetic Algorithms (RCGAs), the PCA or Gram-Schmidt orthogonalization are also often used. See Deb, K. A Population-based Algorithm-generator for Real-parameter Optimization. Technical Report 2003003, Indian Institute of Technology, Kanpur: Kanpur Genetic Algorithms Laboratory (KanGAL), Kanpur, PIN 208 016, India, 2003, and Deb, K., Anand, A. and Joshi, D. A Computationally Efficient Evolutionary Algorithms for Real-parameter Optimisation. Evolutionary Computation, 10(4):371-395, 2002, and Deb, K. and Agrawal, R. B. Simulated Binary Crossover for Continuous Search Space. Complex Systems, 9:115-148, 1995, and Deb, K. and Beyer, H. -G. Self-adaptation in Real-parameter Genetic Algorithms with Simulated Binary Crossover. In Proceedings of Genetic and Evolutionary Computation Conference (GECCO-1999), pages 172-179, 1999, and Eshelman, L. J. and Schaffer, J. D. Real-coded Genetic Algorithms and Interval-schemata. In Proceedings of Foundations of Genetic Algorithms 2 (FOGA-2), pages 187-202, 1993, and Eshelman, L. J., Mathias, K. E. and Schaffer, J. D. Crossover Operator Biases: Exploiting the Population Distribution. In Proceedings of the Seventh International Conference on Genetic Algorithms, pages 354-361, 1997, and Kita, H., Ono, I. and Kobayashi, S. Multi-parental Extension of the Unimodal Normal Distribution Crossover for Real-coded Genetic Algorithms. In Proceedings of Congress on Evolutionary Computation (CEC-1999), volume 2, pages 1581-1587, 1999, and Ono, I. Genetic Algorithms for Optimization Taking Account of Characteristics Preservation (Japanese). PhD thesis, Tokyo Institute of Technology, 1997, and Ono, I. and Kobayashi, S. A Real-coded Genetic Algorithm for Function Optimization Using Unimodal Normal Distribution Crossover. In Thomas Back, editor, Proceedings of the Seventh International Conference on Genetic Algorithms, pages 246-253, 1997, and Ono, I., Kita, H. and Kobayashi, S. A Robust Real-Coded Genetic Algorithm using Unimodal Normal Distribution Crossover Augmented by Uniform Crossover: Effects of Self-Adaptation of Crossover Probabilities. In Proceedings of Genetic and Evolutionary Computation Conference (GECCO-1999), volume 1, pages 496-503, 1999, and Ono, I., Satoh, H. and Kobayashi, S. A Real-coded Genetic Algorithm for Function Optimization Using the Unimodal Normal Distribution Crossover (Japanese). Journal of Aritificial Intelligence, 14(6):1146-1155, 1999, and Takahashi, M. and Kita, H. A Crossover Operator Using Independent Component Analysis for Real-coded Genetic Algorithms. In Proceedings of Congress on Evolutionary Computation (CEC-2001), pages 643-649, 2001, and Tsutsui, S. Sampling Bias and Search Space Boundary Extension in Real Coded Genetic Algorithms. In Proceedings of Genetic and Evolutionary Computation Conference (GECCO-2000), pages 211-218, 2000, and Tsutsui, S. and Goldberg, D. E. Search Space Boundary Extension Method in Real-coded Genetic Algorithms. Information Sciences, 133(3-4):229-247, 2001, and Tsutsui, S. and Goldberg, D. E. Simplex Crossover and Linkage Identification: Single stage Evolution vs. Multi-stage Evolution. In Proceedings of Congress on Evolutionary Computation (CEC-2002), pages 974-979, 2002, and Tsutsui, S., Goldberg, D. E. and Sastry, K. Simplex Crossover and Linkage Learning in Real-coded GAs. In Proceedings of Genetic and Evolutionary Computation Conference (GECCO-2001), page 785, 2001, and Tsutsui, S., Yamamura, M. and Higuchi, T. Multi-parent Recombination with Simplex Crossover in Real Coded Genetic Algorithms. In Proceedings of Genetic and Evolutionary Computation Conference (GECCO-1999), pages 657-664, 1999, and Ueda, T., Koga, N. and Okamoto, M. Efficient Numerical Optimization Technique Based on Real-coded Genetic Algorithm. Genome Informatics, 12:451-453, 2001. All of these references are incorporated by reference herein in their entirety. 
   Unimodal Normal Distribution Crossover (UNDX), Bimodal Normal Distribution Crossover (BNDX), Parent-Centric Recombination (PCX), Modified PCX (mPCX) are examples of the usage of orthogonalization. See Ono, I. and Kobayashi, S. A Real-coded Genetic Algorithm for Function Optimization Using Unimodal Normal Distribution Crossover. In Thomas Bäck, editor, Proceedings of the Seventh International Conference on Genetic Algorithms, pages 246-253, 1997, and Ono, I., Kita, H. and Kobayashi, S. A Robust Real-Coded Genetic Algorithm using Unimodal Normal Distribution Crossover Augmented by Uniform Crossover: Effects of Self-Adaptation of Crossover Probabilities. In Proceedings of Genetic and Evolutionary Computation Conference (GECCO-1999), volume 1, pages 496-503, 1999, and Takahashi, M. and Kita, H. A Crossover Operator Using Independent Component Analysis for Real-coded Genetic Algorithms. In Proceedings of Congress on Evolutionary Computation (CEC-2001), pages 643-649, 2001, and Ueda, T., Koga, N. and Okamoto, M. Efficient Numerical Optimization Technique Based on Real-coded Genetic Algorithm. Genome Informatics, 12:451-453, 2001, and Deb, K. A Population-based Algorithm-generator for Real-parameter Optimization. Technical Report 2003003, Indian Institute of Technology, Kanpur: Kanpur Genetic Algorithms Laboratory (KanGAL), Kanpur, PIN 208 016, India, 2003 which are all incorporated by reference herein in their entirety. However, the usage of orthogonalization in MOO has problems, see  FIG. 7 , where the result of UNDX on SCH1 with 50 dimension is shown. 
   Clearly, the solution sets are not on the Pareto front. However, they seem to be on similar curves. Taking a closer look, we can find out the reason for this problem. If all solutions are on one line in the parameter space, the search power of the orthogonalization will be lost. A rough image is drawn in  FIG. 8 , where it can be seen that the search power for the orthogonal direction is lost if all solutions are on one line. Since all solutions lie on one line which is not the Pareto front, all offspring will be generated on this line. Thus, there is no way to find the Pareto front. The curves that can be seen in  FIG. 7  correspond to the situation described above where the orthogonalization fails. The effect can be regarded as a kind of premature convergence. 
   To avoid this undesirable convergence, perturbation is introduced.  FIG. 9 . shows perturbation added in the orthogonal direction and extension in the principal direction, β and γ being the differences between the maximum value and the minimum value in the principal and the orthogonal direction, δ being the perturbation in the orthogonal direction, and α being the extension for the principal direction. 
   In  FIG. 9 , the “obtained area” is given by the maximum and the minimum value in each direction. The width of the obtained area are β and γ in the principal direction and the orthogonal direction, respectively. If γ becomes zero, γ is set to 0.01 β. Basically, offspring will be generated in this obtained area. However, if offspring are generated in this area only, the above problem can occur. Thus, perturbation is introduced. The “obtained area” is shifted in the orthogonal direction by an offset δ given by δ=±0.25 β. The direction, i.e., “+” or “−”, is chosen randomly. Furthermore, the width β is enlarged: α=1.25 β. The new “generative area” is given by α times γ shifted by δ. Although several parameters were introduced, the performance of VEDA is robust against these parameters. Thus, we fix these parameters through this work. 
   Voronoi-based Stochastic Model 
   Voronoi mesh is defined as follows: 
   Given a set S of m data points in            , Voronoi mesh is the partition of           into m polyhedral cells, v(p)(p ∈ S). Each cell v(p), namely Voronoi mesh of p, is defined as the set of points in           which are closer to p than to any other points in S, or more precisely,
 
 v ( p )={ x ∈             ; dist(   x, p )≦dist( x, q )∀ q ∈S\p}   (Eq. 6)
 
where dist( ) is the Euclidean distance.

   In the generated area in  FIG. 9 , a Voronoi mesh is calculated in the following way:
     1. The rank for all data in the generated area is determined.   2. The probability of each cell is determined by the nearest data.   

   To simplify the generation of the Voronoi mesh, a discrete Voronoi mesh is used, see  FIG. 10 . The space is divided into small cells. The Voronoi mesh will be generated with these small cells.  FIG. 11  shows the sample flow for the generation of the Voronoi mesh. The thicker lines outline the boundary of the Voronoi mesh. 
   First of all, a space is divided into small cells, see  FIG. 11 . The filled circles in  FIG. 11  are data. As the first step, neighbors of each datum are grouped. This step is indicated by 1 in  FIG. 11 . As the second step, from the neighbors indicated by 1, the same procedure is carried out, which is indicated by 2. Until all cells are grouped, this procedure is repeated. If a cell is grouped from different data simultaneously, this cell is grouped to the data with the worse rank. 
   The number of cells in each direction, Ds, is determined as:
 
 D   s   =ε×|N   C |  (Eq. 7)
 
wherein ε and |N C | are respectively a predefined parameter and the number of data in a cluster.
 
Generate Offspring
 
   In the above section describing the Voronoi-based Stochastic Model, all cells in the generative area were assigned a rank. With the assigned rank, the selection probability is calculated. To calculate the selection probability, the geometry distribution is used. Mathematically, the geometry distribution is calculated as follows:
 
 P=P   G (1 −P   G ) r−1   (Eq. 8)
 
wherein PG and r are the given value in the range [0.0, 1.0] and the rank, respectively. For PG=0.2, 0.5, 0.8, the geometry distributions are drawn in  FIG. 12 . In this work, PG=0.8 is used.
 
To generate offspring, a cell is selected according to the above selection probability first. In the selected cell, an offspring is uniformly generated. This flow will be repeated until all necessary offspring are generated.
 
Comparison of VEDA with NSGA-II
 
   The proposed method of the present invention is tested on the test functions SCH1 (n=2, 5), FON2 (n=2, 5) and OKA4 (n=2). The used parameters are shown in Table 2-a, wherein ε corresponds to the accuracy of a model. 
   
     
       
             
           
             
             
             
           
         
             
               TABLE 2-a 
             
             
                 
             
             
               Parameters for VEDA. 
             
             
                 
             
           
           
             
                 
             
           
        
         
             
                 
               Number of Data 
               100 
             
             
                 
               Number of Offspring 
               100 
             
             
                 
               ε 
               1.0 for n = 2 
             
             
                 
                 
               0.5 for n = 5 
             
             
                 
               Maximum Iterations 
                10 
             
             
                 
                 
             
           
        
       
     
   
   Since one of the merits of VEDA is less fitness evaluations, only 1000 fitness evaluations are allowed. To reduce computational cost in n=5, the accuracy is reduced from ε=1.0 to 0.5. 
   To compare the results, NSGA-II (non-dominated Sorting Genetic Algorithm II), see Deb, K., Pratap, A., Agarwal, S. and Meyarivan, T. A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2):182-197, 2002, is also tested with only 1000 fitness evaluations. The parameters for NSGA-II are shown in Table 2-b. 
   
     
       
             
           
             
             
             
           
         
             
               TABLE 2-b 
             
             
                 
             
             
               Parameters for NSGA-II. 
             
             
                 
             
           
           
             
                 
             
           
        
         
             
                 
               Number of Individuals 
               100 
             
             
                 
               Crossover 
               One-Point Crossover 
             
             
                 
               Crossover Rate 
               0.9 
             
             
                 
               Mutation Rate 
               0.01 
             
             
                 
               Number of Bits per a design 
               20 
             
             
                 
               parameter 
             
             
                 
               Maximum Iterations 
               10 
             
             
                 
                 
             
           
        
       
     
   
     FIG. 13  shows the generated stochastic model for OKA4 in the parameter space, the Pareto front being equal to x 2 =x 1 ±4√{square root over (x 1 )}+4. The gray-level corresponds to the rank. Clearly, the generated model has a high probability near the Pareto front in the PS. 
   The simulation results of SCH1, FON2 and OKA4 for VEDA and NSGA-II are shown in  FIG. 14 . Since OKA4 is designed with a non-linear Pareto front in the PS, the solutions in the PS on OKA4 by VEDA and NSGA-II are shown in  FIG. 15 . In the figures, all parents are plotted. Since many parents in NSGA-II don&#39;t converge to the shown area, the number of individuals seems to be less than VEDA. 
   On the SCH1 and FON2 with 2 dimensional cases, the differences are small. However, VEDA is better than NSGA-II. On the SCH1 and FON2 with 5 dimensional cases, the performance of both is not sufficient. However, clearly VEDA shows a better result than NSGA-II. On the OKA4, the difference is clear. VEDA nearly converges to the Pareto front but NSGA-II does not. 
   From these results, VEDA shows better performance than NSGA-II in the early generation. This indicates that a real-world optimization problem may be solved without many fitness evaluation by the proposed VEDA. However, VEDA has also drawbacks, e.g. the increase of the computational complexity in higher dimensional cases. Like other EDAs, the current capacity of VEDA is about 10 dimension. Since many real-world optimization problems have a large number of design parameters, the current capacity of EDAs is not sufficient. 
   Comparison of VEDA with Other EDAs for MOO 
   Although EDAs are mainly studied on single objective optimization (SOO) problems, recently some methods have applied EDAs to MOO. Thierens and Bosman have applied Mixture-based Iterated Density Estimation Evolutionary Algorithm (MIDEA) to MOO and Costa and Minisci have applied Parzen-based EDA (PEDA) to MOO. Since both of them have used the same test functions, we use these test functions for the comparison. Since MIDEA did not show a sufficient quality on test function ZDT4 (n=10), ZDT4 is not used in this comparison. The maximum number of iterations is set as Table 3. The numbers of data and offspring are 100 and the value of ε is 1.0 for test functions FON2 and KUR1 and 0.1 for test function DEB4. See Thierens, D. and Bosmann, P. A. N. Multi-Objective Mixture-based Iterated Density Estimation Evolutionary Algorithms. In Proceedings of Genetic and Evolutionary Computation Conference (GECCO-2001), pages 663-670, 2001 for MIDEA, which is incorporated by reference herein in its entirety and Costa et al. for PEDA. 
   
     
       
             
           
             
             
             
             
             
           
         
             
               TABLE 3 
             
           
           
             
                 
             
             
               Maximum number of iterations. 
             
           
        
         
             
                 
               Method 
               FON2 (n = 3) 
               KUR1 (n = 3) 
               DEB4 (n = 10) 
             
             
                 
                 
             
             
                 
               MIDEA 
               3754 
               10762 
               8426 
             
             
                 
               PEDA 
               3100 
               11000 
               8300 
             
             
                 
               VEDA 
               3100 
               10700 
               8300 
             
             
                 
                 
             
           
        
       
     
   
   The results and a comparison of VEDA with state of the art methods are shown in  FIG. 16 . By comparison with the results from Table 3, VEDA is superior to MIDEA and PEDA on FON2 and KUR1 but VEDA is inferior to them on DEB4. Taking a closer look, the reason for DEB4 is caused by the generation of many offspring in an infeasible region. Since the current VEDA has no restriction for perturbation, the design parameters can easily become less than zero, thus in the case of design infeasible. By adding additional constraints in VEDA, the performance may be improved. 
   Piece-wise Linear Approximation of Pareto Front with VEDA 
   In the state of the art, nearly all MOO methods output a set of solutions to represent the Pareto front. There is no MOO method to output a mathematical description of the Pareto front in the parameter space. Jin and Sendhoff have proposed a way to represent the Pareto front by piece-wise linear functions in the PS. Since many test functions have a linear Pareto front in the PS, they have exploited this characteristic to get more accurate solutions. See Jin, Y. and Sendhoff, B. Connectedness, Regularity and the Success of Local Search in Evolutionary Multi-Objective Optimization. In Proceedings of Congress on Evolutionary Computation (CEC-2003), volume 3, pages 1910-1917, 2003 which is incorporated by reference herein in its entirety. 
   One of the implicit characteristics of VEDA is the availability of a mathematical description instead of just a solution set. Since VEDA uses the PCA and the maximum and minimum values in each axis, VEDA can output mathematical description. As an example, the output of VEDA on SCH1 (n=2) is shown in Table 4. 
   
     
       
             
           
             
             
             
           
         
             
               TABLE 4 
             
             
                 
             
             
               Mathematical output of VEDA. 
             
             
                 
             
           
           
             
                 
             
           
        
         
             
                 
               Principal Axis 
               {0.701764, 0.712410} 
             
             
                 
               Orthogonal Axis 
               {0.712410, −0.701764} 
             
             
                 
               Minimum Value in Principal Axis 
               −0.042017   
             
             
                 
               Maximum Value in Principal Axis 
               2.841821 
             
             
                 
               Minimum Value in Orthogonal Axis 
               −0.166802   
             
             
                 
               Maximum Value in Orthogonal Axis 
               0.208333 
             
             
                 
                 
             
           
        
       
     
   
   The graphical output is shown in  FIG. 17 . In  FIG. 17 , the solid rectangle is the output of VEDA and the dotted line is the Pareto front in the PS. The principal axis shows the gradient of the solution set and the minimum and maximum values indicate the range. By the information of axis and its range, the mathematical output seems to be possible. 
   Summary of the Advantages of the Present Invention 
   The invention proposes the Voronoi-based Estimation of Distribution Algorithm (VEDA) for tackling MOO problems. In VEDA, the offspring distribution with the ranking information is directly used to generate promising solutions. Using the concept of Voronoi mesh, the stochastic model is made and the promising solutions are generated according to the generated stochastic model. 
   On several test functions, the performance of VEDA has been compared with a state-of-the-art MOO algorithm called NSGA-II proposed in Deb, K., Pratap, A., Agarwal, S. and Meyarivan, T. A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2):182-197, 2002 which is incorporated by reference herein in its entirety. The performance of VEDA is better than NSGA-II. With a small number of fitness evaluations, VEDA can identify the Pareto front in most cases. This indicates the possibility to reduce the computational cost in real-world problems using VEDA. Although VEDA shows good performance on several test functions, the remaining problem is the increase of the computational complexity in higher dimensional cases. Additionally, VEDA should be modified to consider additional constraints. 
   In the VEDA, several parameters are introduced. Since the performance of the VEDA seems to be stable against the parameters, the fixed values are used in this invention. However, the parameter dependency should be taken into account. Although most of MOO methods in the state of the art output only solution set, VEDA can output a mathematical description.