Postadress: Besöksadress: Telefon:
Box 1026 Gjuterigatan 5 036-10 10 00 (vx) 551 11 Jönköping
Multi-Objective Optimization
using Genetic Algorithms
Kaveh Amouzgar
THESIS WORK 2012
PRODUCT DEVELOPMENT AND MATERIALS
ENGINEERING
Postadress: Besöksadress: Telefon:
Box 1026 Gjuterigatan 5 036-10 10 00 (vx) 551 11 Jönköping
Multi-Objective Optimization
using Genetic Algorithms
Kaveh Amouzgar
This thesis work has been carried out at the School of Engineering in
Jönköping in the subject area Product Development and Materials Engineering.
The work is a part of the master’s degree.
The authors take full responsibility for opinions, conclusions and findings
presented.
Supervisor: Niclas Strömberg
Scope: 30 ECTS credits
Date: 2012-05-30
Abstract
In this thesis, the basic principles and concepts of single and multi-objective Ge-netic Algorithms (GA) are reviewed. Two algorithms, one for single objective and the other for multi-objective problems, which are believed to be more efficient, are described in details. The algorithms are coded with MATLAB and applied on several test functions. The results are compared with the existing solutions in literatures and shows promising results. Obtained pareto-fronts are exactly similar to the true pareto-fronts with a good spread of solution throughout the optimal region. Constraint handling techniques are studied and applied in the two algorithms. Constrained benchmarks are optimized and the outcomes show the ability of algorithm in maintaining solutions in the entire pareto-optimal region. In the end, a hybrid method based on the combination of the two algorithms is introduced and the performance is discussed. It is concluded that no significant strength is observed within the approach and more research is required on this topic. For further investigation on the performance of the proposed techniques, implementation on real-world engineering applications are recommended.
Keywords
Single Objective Optimization, Multi-objective Optimization, Constraint Han-dling, Hybrid Optimization, Evolutionary Algorithm, Genetic Algorithm, Pareto-Front, Domination.
Contents
1 Introduction 1
1.1 Background . . . 1
1.2 Purpose and aims . . . 1
1.3 Delimitations . . . 1 1.4 Outline . . . 2 2 Theoretical background 3 2.1 What is Optimization? . . . 3 2.2 Single-Objective Optimization . . . 4 2.2.1 Evolutionary Method . . . 4
2.2.2 Genetic Algorithm Concept . . . 4
2.2.3 Genetic Algorithm Principles . . . 5
2.2.4 Real Parameter GA . . . 6
2.2.5 Generalization Generation Gap Algorithm (G3) . . . 7
2.2.6 Parent-Centric Recombination Operator (PCX) . . . 8
2.2.7 Constraint Handling . . . 9
2.3 Multi-objective Optimization . . . 10
2.3.1 Multi-Objective Optimization Formulation . . . 10
2.3.2 Multi-Objective Optimization Definitions . . . 11
2.3.3 Approaches Towards Non-Dominated Set . . . 14
2.3.4 Approaches Towards Multi-Objective Optimization . . . 14
2.3.5 Evolutionary Algorithms . . . 15
2.3.6 MOEA Techniques . . . 16
2.3.7 Comparison of MOEAs . . . 17
2.3.8 SPEA2: Improved Strength Pareto Evolutionary Algorithm 18 2.3.9 Overall SPEA2 Algorithm . . . 18
2.3.10 Constraint Handling . . . 21
3 Implementation 24
3.1 Single Objective . . . 24
3.1.1 Unconstrained Test Functions . . . 24
3.1.2 Constrained Test Functions . . . 26
3.2 Multi objective . . . 29
3.2.1 Unconstrained Test Functions . . . 29
3.2.2 Constrained Test Functions . . . 32
3.3 Hybrid Approach . . . 35 4 Test Results 36 4.1 Single Objective . . . 36 4.1.1 Unconstrained Functions . . . 36 4.1.2 Constrained Functions . . . 38 4.2 Multi-Objective . . . 40 4.2.1 Unconstrained Functions . . . 40 4.2.2 Constrained Functions . . . 51 4.3 Hybrid Approach . . . 52 5 Conclusion 54 6 Bibliography 55 A Hand Calculation of G3 Algorithm with Constraints 59 B Hand Calculation of SPEA2 Algorithm 63 B.1 Constraint Handling Method of SPEA2 Algorithm . . . 70
List of Figures
1 Trade-off curve . . . 10
2 Min-Min pareto-front . . . 13
3 The welded beam problem. . . 28
4 Tension/compression string problem. . . 29
5 Convergence of Schwefel’s function . . . 37
6 Convergence of Rosenbrock function . . . 37
7 Convergence of Test Function 1 . . . 38
8 Convergence of welded beam problem . . . 38
9 Convergence of Tension/Compression Spring . . . 38
10 Pareto-front of Exercise 14, single objective . . . 40
11 Pareto-front of Exercise 14, multi objective . . . 40
12 Pareto-front of Kursawe test function . . . 41
13 Pareto-front of ZDT1 test function . . . 41
14 Pareto-front of ZDT2 test function . . . 41
15 Pareto-front of ZDT3 test function . . . 41
16 Pareto-front of ZDT4 test function . . . 42
17 Pareto-front of ZDT6 test function . . . 42
18 Scatter-plot matrix of Kursawe test function . . . 45
19 Scatter-plot matrix of ZDT1 test function . . . 46
20 Scatter-plot matrix of ZDT2 test function . . . 47
21 Scatter-plot matrix of ZDT3 test function . . . 48
22 Scatter-plot matrix of ZDT4 test function . . . 49
23 Scatter-plot matrix of ZDT6 test function . . . 50
25 Pareto-front of OSY test function . . . 51
26 Pareto-front of SRN test function . . . 51
27 Pareto-front of TNK test function . . . 51
28 ZDT1 Hybrid and random archive population . . . 52
29 ZDT3 Hybrid and random archive population . . . 52
30 ZDT6 Hybrid and random archive population . . . 52
31 Min-Ex pareto-front and initial solutions . . . 63
32 Constrained Min-Ex pareto-front, feasible region and initial solutions 70
List of Tables
1 Results of unconstrained test functions, single objective. . . 362 Comparison of welded beam results . . . 38
3 Comparison of tension/compression spring results . . . 39
4 Pre-defined parameters of unconstrained SPEA2 . . . 41
5 Pre-defined parameters of constrained SPEA2 . . . 51
6 Random initial solutions for G3 algorithm hand calculation example 60 7 Current and external initial random population of SPEA2 . . . 63
8 Fitness assignment procedure of SPEA2 . . . 65
1
Introduction
Although substantial amount of search in optimization is conducted with regards to single objective problems, optimization problems with multi conflicting objec-tives are inevitable in many topics specially engineering applications. Two main methods have been proposed by scientist for solving multi-objective optimization problems: 1) Classical method, 2) Evolutionary algorithms. Classical methods are able to reach one optimal solution at each run, while evolutionary algorithms are based on a population of solutions which will hopefully lead to a number of opti-mal solutions at every generation. The evolutionary algorithm method which had shown benefits over the classical approach can be categorized in several categories. Genetic Algorithm is one of the methods that mimic the evolution of genes and chromosomes.
1.1
Background
Previous works by Beasley and Bull (1993); Coello (2007); Deb (1995, 2001, 2002, 2004); Deb et al. (2001); Fonseca and Fleming (1993); Haupt et al. (2004); Kim et al. (2004); Kukkonen (2006); Man et al. (1996); Zitzler and Thiele (1998); Zitzler et al. (2001) on the theory, concepts and algorithms of single and multi-objective optimization using evolutionary algorithms.
Previous studies by Deb (2000); Deb et al. (2002); Jimnez et al. (1999); Kuri-Morales and A.Gutierrez-Garcia (2002); Mezura-Montes et al. (2003); T. Ray (2001) on constraint handling methods.
Test functions and their comparison has been studied by (Binh and Korn, 1997); Deb (1991, 1999); Gamot and Mesa (2008); Kuri-Morales and A.Gutierrez-Garcia (2002); Kursawe (1991); Osyczka and Kundu (1995); Srinivas and Deb (1994); Tanaka and Watanabe (1995) and Zitzler et al. (2000).
1.2
Purpose and aims
The aim is to develop a fast and efficient multi-objective optimization technique by using GA (Genetic Algorithm) method, in order to solve multi-objective opti-mization problems with constraints. Several benchmarks are to be optimized with the developed algorithm. Also a hybrid approach will be suggested and further studied. MATLAB shall be used to code the algorithm.
1.3
Delimitations
Genetic Algorithm is the only method used in developing the technique. Other evolutionary methods like Evolution strategies, Evolutionary programming and Genetic Programming are not considered in the thesis.
1.4
Outline
The organization of the thesis is in five different sections:
Section 1: Introductory, background and purpose is described.
Section 2: General theory of optimization, single and multi-objective optimiza-tion is explained. Evoluoptimiza-tionary algorithms especially Genetic Algorithms are discussed in details and two algorithms are suggested. Constraint handling techniques and a hybrid method are theoretically defined.
Section 3: A short description on implementation of algorithms and a number of benchmarks will be presented in this section.
Section 4: The results obtained from the benchmarks are illustrated and com-pared with references.
2
Theoretical background
2.1
What is Optimization?
Optimization is a process of making things better. Life is full of optimization problems which all of us are solving many of them each day in our life. Which route is closer to school? Which bread is better to buy having the lowest price while giving the required energy? Optimization is fine-tuning the inputs of a process, function or device to find the maximum or minimum output(s). The inputs are the variables, the process or function is called objective function, cost function or Fitness value (function) and the output(s) is fitness or cost (Haupt et al., 2004).
In the thesis minimization of cost is tackled, in functions which maximum of cost is required, by slapping a minus in front of objective function, the output will be minimized. Therefore all the problems and functions in the thesis are addressed as minimization problem.
When only one objective function involves in the problem, it is called single-objective optimization, however in most real world problems more than one objec-tive function is required to be optimized, and therefore these problems are named multi-objective optimization.
Deb (2001) classified optimization solving methods into following two major cat-egories:
• Classical methods
• Evolutionary methods
The classical methods commonly use a single random solution, updated in every iteration by a deterministic procedure to find the optimal solution. These methods are classified into two distinct groups: direct methods where only the objective function and the constraints value are used to find the optimum and gradient-based methods whereas the first and second derivative of objective function and/or constraints are applied to find the search direction and optimal solution (Deb, 2001).
2.2
Single-Objective Optimization
2.2.1 Evolutionary Method
This method was inspired by the evolutionary process of human being and the interests for imitating living being is increasing since 19600s. Evolutionary method mimics the evolution principle of nature which results in a stochastic search and optimization algorithm. It also can out pace the classical method in many ways (Gen and Cheng, 1997).
Evolutionary method (algorithm) uses an initial population of random solutions in each iteration, instead of using a single solution as in classical method. This initial population is updated in each generation to finally converge to a single optimal solution. Having a population of optimum solution in a single simulation run, is a unique characteristic of the method in solving multi-objective optimization problems (Deb, 2001).
Gen and Cheng (1997) divides the method into three main types: genetic al-gorithm, evolutionary programming and evolutionary strategy while Deb (2001) describes an additional type to the three above-mentioned: genetic programming.
2.2.2 Genetic Algorithm Concept
Many real-world optimization problems are extremely difficult and complex in terms of number of variables, nature of the objective function, many local optimal, continuous or discrete search space, required computation time and resources, etc. in various domains including service, commerce and engineering.
Genetic algorithm was first introduced by John Holland (1975) in 19700s, whereas one of his students David Goldberg had an important contribution in popularizing this method in his dissertation by solving a complex problem (Haupt et al., 2004).
Genetic algorithm is an inspiration of the selection process of nature, where in a competition the stronger individuals will survive (Man et al., 1996). In nature each member of a population competes for food, water and territory, also strive for attracting a mate is another aspect of nature. It is obvious the stronger in-dividuals have a better chance for reproduction and creating offspring, while the poor performers have less offspring or even non. Consequently the gen of the strong or fit individuals will increase in the population. Offspring created by two
fit individual (parents) has a potential to have a better fitness compared to both parents called super-fit offspring. By this principle the initial population evolves to a better suited population to their environment in each generation (Beasley and Bull, 1993).
2.2.3 Genetic Algorithm Principles
As mentioned, in genetic algorithm unlike other classical methods, a population of random solution is selected. Each solution of the problem is represented as a set of parameters which are known as genes. Joining genes create a binary bit string of values, denoting each member of population referred as chromosome. A chromosome evolves through iterations, called generation (Gen and Cheng, 1997).
After representation a fitness or objective function is required. Also during the run a selection of parents for reproduction and recombination for creating offspring is essential. These aspects are called GA0s operators (Beasley and Bull, 1993). A selection or reproduction operator during reproduction phase of GA selects par-ents from population which they create offspring by recombination comprising next generation. The main objective of selection operator is to keep and du-plicate the fit solutions and eliminate the poor chromosomes, while keeping the size of population constant. Deb (2001) describes some schemes for achieving the above objective: tournament selection, proportionate selection, ranking selection, roulette wheel selection (RWS) and stochastic universal selection (SUS). It is ob-vious this operator cannot create new chromosomes to the initial population, it only make copies of good solutions.
In reproduction phase the two parents nominated by selection operator recombine to create one or more offspring with crossover or mutation operators.
They are number of different crossover operators in literature but the main concept is selecting two strings of solution (chromosomes) from the mating pool of selection operator and exchanging some portion of these two strings from a random selected point(s). Single point cross over is one basic type of this operator for binary GA (Deb, 2001).
A mutation operator is applied to individual solutions after cross over operator which a gene(s) is randomly changed in a string with a small probability to create a new chromosome. The aim of this operator is to maintain the diversity of the
population and increase the possibility of not losing any potential solution and find the global optimal, while cross over operator is a technique of rapid exploration of search space (Beasley and Bull, 1993).
To sum up, the selection operator selects and maintains the good solutions; while crossover recombines the fit solutions to create a fitter offspring and mutation operator randomly alter a gene or genes in a string to hopefully find a better string (Deb, 2001).
2.2.4 Real Parameter GA
There are some difficulties in binary-coded GAs, including inability to solve the problems where the values of variables have continuous search space or when the required precision is high. According to Deb (2001) hamming cliffs related to certain strings (01111 or 11110) is one of the difficulties where altering to a near neighbour string requires changes in many genes. He also claims necessity of large strings (chromosomes with many genes) in order to fulfil a necessary precision which in result increases the size of population, as another struggle for binary GAs. Therefore using floating point numbers to represent the variables in most problems is more logical which requires less storage than binary coded strings. In addition, since there is no need for decoding the chromosomes before evaluation of objective function in selection phase the real parameter GA (in some literature called continuous GA) is inherently faster than binary GA (Haupt et al., 2004).
Since the real value of parameters are directly used to find the fitness value in selection operator and there is no decoding to a string in real parameter GAs, this operator does not alter with binary GA selection operators and the same operators can be used in real parameter GAs. On the other hand, since the cross over and mutation operators used in binary GAs are based on strings and alteration in genes (bits), new cross over and mutation operators shall be defined for this type of GA.
Deb (2001) outlines some real parameter crossover operators such as linear crossover, naive crossover, blend crossover (BLX), simulated binary crossover (SBX), fuzzy recombination operator, unimodal normally distributed crossover (UNDX), sim-plex crossover (SPX), fuzzy connectives based crossover and unfair average crossover. Other cross over operators including parent centric crossover (PCX), modified PCX (mPCX) are recommended in literatures (Deb and Joshi, 2001).
Since in real parameter crossover operator two or more parents directly recombine to create on or more offspring and it has the same concept as mutation operator, a question comes up: Is there a good reason for using a mutation operator along with crossover operator? The debate still remains, however Deb (2001) argues, the different between these two operators is in the number of parent solutions selected for perturbation. He claims if offspring is created from one parent the operator is mutation while offspring created from more than one parent is crossover. He also mentions some common mutation operators in his book: Random mutation, non-uniform mutation, normally distributed mutation and polynomial mutation.
2.2.5 Generalization Generation Gap Algorithm (G3)
Deb (2002) proposes a population-based, four steps, real-parameter optimization algorithm-generator called Generalization Generation Gap (G3) model. In the same paper performance of the G3 algorithm is studied on three commonly used test problems and is compared with a number of evolutionary and classical opti-mization algorithms, also Deb (2004) performs a systematic parametric study on G3 model, both of these studies concludes to out performance of the algorithm to a number of existing classical and evolutionary algorithms. The G3 algorithm is coded in MATLAB for solving single-objective optimization problems in the thesis.
Generalization generation gap algorithm is modified steady-state GA to make it computationally faster, in which in every iteration only two new solutions are updated in the GA population. This model preserves elite solutions from previ-ous generation (Deb, 2002). Four plans are used which are Selection Plan (SP), Generation Plan (GP), Replacement Plan (RP) and Update Plan (UP).
The steps in algorithm are as follows:
Step 1 (SP): From the population B (set B) of size N, the best parent and µ − 1 other parents are randomly selected. These µ solutions create Set P.
Step 2 (GP): λ offspring are created from µ parents in set P, with using any recombination operator, which creates Set C.
Step 3 (RP): r random parents are chosen from set B, which creates Set R.
Step 4 (UP): r random parent chosen in step 3 (RP) are replaced with r best solutions from the combined set C ∪ R (set RC), in set B.
Several parametric studies such as Deb (2001, 2002, 2004); Kita (2001), compare the performance of G3 model to other evolutionary algorithms, and in all of the studies G3 model has shown a better performance and robustness. Also using different recombination operators has been examined and the overall result shows faster computation time and lower number of evaluation required to meet a desired accuracy of a parent centric recombination operator (PCX) proposed by Deb et al. (2001), which will be briefly described.
2.2.6 Parent-Centric Recombination Operator (PCX)
Deb et al. (2001) suggests a variation operator (combination of crossover and mu-tation operator) for this algorithm, called parent centric recombination operator (PCX). A parent centric operator ensures identically of population mean of the total offspring population to that of the parent population while mean centric operators preserve the mean between the participating parents and resulting off-spring. The paper states the benefit of parent centric recombination operators over mean centric operator, as the parents are selected from the fittest solution in selection plan and in most real parameter optimization problems it is assumed that the solutions near the parents can be the potential good solutions. Therefore, creating new solutions close to parents as how it is in PCX is a steady and reliable search technique.
The mean vector ~g of the chosen parents is computed. For each offspring, one parent ~x(p) is chosen with equal probability. The direction vector ~d = ~g − ~x(p) is
calculated. Thereafter from each of the µ − 1 parents perpendicular distance Di to
the line ~d(p) are computed and their average ¯D is found. The offspring is created as follows: ~ y = ~x(p)+ wζd~(p)+ µ X i=1,i6=p wηD~¯e(i)
where ~e(i) are the µ − 1 orthonormal bases that span the subspace perpendicular
to ~d(p). Parameters w
ζ and wη are zero-mean normally distributed variables with
variance w2
2.2.7 Constraint Handling
Most existing constraint handling methods in literatures are classified in five cat-egories which Deb (2001) describes them briefly:
• Method based on preserving feasibility of solutions.
• Method based on penalty functions.
• Methods biasing feasible over infeasible solutions.
• Methods based on decoders.
• Hybrid methods.
In the thesis, the method based on penalty function is used for single-objective optimization.
Penalty function method transforms a constrained optimization problem to an unconstrained problem usually by using an additive penalty term or penalty mul-tiplier. Penalty method can also be categorized in seven different type:
• Death Penalty • Static Penalties • Dynamic Penalties • Annealing Penalties • Adaptive Penalties • Segregated GA
In the Static Penalty method which is implemented in this section, the penalty parameters do not change within generations and is only applied to infeasible solutions.
There are number of approaches in this method suggested by authors but Morales et al. (1997) penalizes the objective function of infeasible solutions by using the information on the number of violated constraints. His approach is formulated as follows:
F (x) = ( f (x), if xisf easible, K −Ps i=1 K m , otherwise.
where s is the number of non-violated constraints and m is the total number of constraints. K is a large positive constant. Morales et al. (1997) uses 1 × 109 for this constant which should be large enough to assign a bigger fitness to infeasible solution compared to feasible individual. A simple single-objective constrained optimization problem is solved for one generation by using a step by step hand calculation of G3 algorithm in appendix A.
2.3
Multi-objective Optimization
In real world applications, most of the optimization problems involve more than one objective to be optimized. The objectives in most of engineering problems are often conflicting, i.e., maximize performance, minimize cost, maximize reliability, etc. In the case, one extreme solution would not satisfy both objective functions and the optimal solution of one objective will not necessary be the best solution for other objective(s). Therefore different solutions will produce trade-offs between different objectives and a set of solutions is required to represent the optimal solutions of all objectives.
Figure 1 shows the trade-off curve of decision making involved in buying a house problem.
The trade-off curve reveals that considering the extreme optimal of one objective (price) requires a compromise in other objective (house area). However there exists number of trade-off solutions between the two extreme optimal, that each are better with regards to one objective.
2.3.1 Multi-Objective Optimization Formulation
Basically a multi-objective optimization problem has more than one objective function, in engineering problems usually two objectives, to be optimized. In the thesis, minimization problems with only two objectives investigated, while maximization problems are transformed to minimizing optimization types.
Figure 1: Trade-off solutions illustrated for a house-buying decision-making
including inequality, equality and/or variable bounds to be satisfied. However in real engineering applications usually more than one constraint is involved in the problem. A general formulation of a multi-objective optimization problem is defined as follows: Minimize/Maximize fm(x), m = 1, 2, ..., M ; Subject to gj(x) ≤ 0, j = 1, 2, ..., J ; hk(x) = 0, k = 1, 2, ..., K; x(L)i ≤ xi ≤ x (U ) i , i = 1, 2, ..., n.
2.3.2 Multi-Objective Optimization Definitions
In order to fully understand multi-objective optimization problems (MOOP), al-gorithms and concepts some definitions must be clarified.
• Decision variable and objective space: The variable bounds of an opti-mization problem restrict each decision variable to a lower and upper limit which institutes a space called decision variable space.
In multi-objective optimization values of objective functions create a mutli-dimensional space called objective space. Each decision variable on variable
space corresponds to a point in objective space.
• Feasible and infeasible solutions: A solution that satisfies all the con-straints (inequality and equality) and variable bounds is referred to as a feasible solution. On the other hand, a solution that does not satisfy all constraints and variable bounds is called an infeasible solution.
• Ideal objective vector: If x∗(i) is a vector of variables that optimizes (minimize or maximize) the ith objective in a multi-objective optimization problem with M conflicting objectives:
∃ x∗(i)∈ Ω, x∗(i)=hx∗(i)1 , x∗(i)2 , ..., x∗(i)M iT : fi(x∗(i)) = OP T fi(x).
Then, the vector
z∗ = f∗ = [f1∗, f2∗, ..., fM∗ ]T
where fM∗ is the optimum of the Mth objective function, is ideal for a multi-objective optimization problem and the point in <n which determines this
vector is the ideal solution, therefore called the ideal objective vector. Gen-erally, ideal objective vector is a solution that does not exist. The reason is that the optimal solution of each objective in a MOP is not necessary the same solution. However, if the optimal solutions of all objectives are identical the ideal vector is feasible.
• Utopian objective vector: A vector that all of its components are marginally smaller (in case of minimization MOP) than that of the ideal objective vector is called utopian objective vector. In other words:
∀ i = 1, 2, ..., M : zi∗∗= zi∗− i, i > 0.
• Linear and non- linear MOOP: If all objectives and constraints are linear the problem is named a linear optimization problem (MOLP). In contrast, if one or more of the objectives and/or constrains are non-linear the problem in non-linear MOOP (Deb, 2001).
• Convex and Non-convex MOOP: The problem is convex if all objective functions and feasible region are convex. Therefore a MOLP problem is convex (Deb, 2001).
Convexity is an important issue in MOOPs, where in non-convex problems the solutions obtained from a preference-based approach will not cover the
non-convex part of the trade-off curve. Moreover many of the existing algo-rithms can only be used for convex problems. Convexity can be defined on both of spaces (objective and decision variable space). A problem can have a convex objective space while the decision variable space is non-convex.
• Domination (dominated, dominating and non-dominated): Most of real world applications consist of conflicting objectives. Optimizing a solution with respect to one objective will not result in an optimal solution regarding the other objective(s). For a M objective MOP, the operator / between two solutions i and j as i / j is translated as solution i is better than solution j on a particular objective. Also, i . j means that solution i is worse than solution j on this objective. Therefore, if the MOP is a minimization case, the operator / denotes < and vice versa. Now a general definition of domination for both minimization and maximization MOP can be made: A feasible solution x(1) is said to dominate another feasible solution x(2) (or
mathematically x(1) x(2)), if and only if:
1. The solution x(1) is no worse than x(2) with respect to all objectives
value, or fj(x(1)) 7 fj(x(2)) for all j = 1, 2, ..., M .
2. The solution x(1) is strictly better than x(2) in at least one objective
value, or f¯j(x(1)) C f¯j(x(2)) for at least one ¯j ∈ {1, 2, ..., M }.
Therefore solution x(1)dominates solution x(2), solution x(1)is non-dominated by solution x(2) or solution x(2) is dominated by solution x(1).
• Pareto- optimal set (non-dominated set): A solution is pareto-optimal if it is not dominated by any other solution in decision variable space. The pareto-optimal is the best known (optimal) solution with respect to all objec-tives and cannot be improved in any objective without worsening in another objective. The set of all feasible solutions that are non-dominated by any other solution is called the pareto-optimal or non-dominated set. If the non dominated set is within the entire feasible search space, it is called globally pareto-optimal set. In other words, for a given MOP, the pareto-optimal set, P∗, is defined as:
P∗ = {x ∈ Ω | ¬∃ x0 ∈ Ω F (x0) F (x)}.
• Pareto-front: The values of objective functions related to each solution of a pareto-optimal set in objective space is called pareto-front. In other words,
Figure 2: Pareto-front of a Min–Min problem
for a given MOP, F (x), and pareto-optimal set, P∗, the pareto-front, P F∗ is given by:
P F∗ := {u = F (x) | x ∈ P∗}.
Figure 2 illustrates a typical pareto-front of a two objective minimizing type optimization problem in objective space.
Since the concept of domination enables comparison of solutions with respect to multi-objectives, most of multi-objective optimization algorithms practice this concept to obtain the non-dominated set of solutions, consequently the pareto-front.
2.3.3 Approaches Towards Non-Dominated Set
They are several methods and algorithms towards finding the non-dominated set of solutions from a given population in an optimization problem. Deb (2001) describes three of the most common methods in his book from a naive and slow to an efficient and fast approach.
Approach 1: Naive and slow
Approach 2: Continuously updated
Approach 3: Kung et al.s efficient method
Approach 3 has the least computational complexity among the three and according to Kung and Luccio (1975) is the most efficient method. In all methods the concept of domination is used to compare the solution with respect to different objective functions.
2.3.4 Approaches Towards Multi-Objective Optimization
Extensive studies have been conducted in multi-objective optimization algorithms. But most of the researches avoid the complexity in the true multi-objective opti-mization problem by transforming the problem into single- objective optiopti-mization with the use of some user defined parameters.
Deb classifies the approaches towards solving multi-objective optimization in two groups.
• Ideal multi-objective optimization, where a set of solutions in form of a trade-off curve is obtained and the desired solution is selected according to some higher level information of problem.
• Preference based multi-objective optimization, which by using the higher level information a preference vector transforms the multi-objective prob-lem to a single-objective optimization. The optimal solution is obtained by solving the single-objective problem.
The ideal approach is less subjective compared to preference based approach, where analysis of non-technical, qualitative and experimental information is re-quired to find the preference vector. Therefore the second approach will not be further discussed in the thesis.
In absence of higher level information in an optimization problem within ideal approach none of the pareto-optimal solutions is preferred over others. Therefore in the ideal approach the main objective is to converge to a set of solution as close as possible to true pareto-optimal set, which is the common objective of all optimization tasks. However, diversity in the obtained pareto-optimal set is the second objective specific to multi-objective problems. With a more divers set of solutions that covers all parts of pareto-front in objective space, the decision making process at the next level using the higher level information is easier. Since two spaces are involved in MOOP, diversity of solutions in both decision and objective space is defined. Solutions with a large Euclidean distance in variable and objective space are referred as divers set of solutions in variable and objective space, respectively. The diversity in the two spaces are often Symmetric, however in complex and non-linear problems this property may not be true. Hence, Deb (2001) assumes that there are two goals in multi-objective optimization:
1. To find a set of non-dominated solutions with the least distance to pareto-optimal set.
2. To have maximum diversity in the non-dominated set of solutions.
Recall from section 2.1, that classifies optimization solving methods into; classi-cal and evolutionary method , the classification is also valid for multi-objective optimization problems.
In the classical method objectives are transformed to one objective function by means of different techniques. The easiest and probably most common is the weighted sum method which the objectives are scalarized to one objective by multiplying the sum of objectives to a weight vector (Deb, 2001). Other techniques are such as considering all objectives except one as constraints and limiting them by a user defined value ( − constraint) (Haimes and A., 1971). Deb (2001) very well presents some of the most important classical methods in one chapter of the book.
2.3.5 Evolutionary Algorithms
The characteristic of evolutionary methods which use a population of solutions that evolve in each generation is well suited for multi-objective optimization prob-lems. Since one of the main goals of MOOP solvers is to find a set of non-dominated solutions with the minimum distance to pareto-front, evolutionary algorithms can generate a set of non-dominated solutions in each generation.
Requirement of little prior knowledge from the problem , less vulnerability to shape and continuity of pareto-front, easy implementation, robustness and the ability to be carried out in parallel are some of the advantages of evolutionary algorithms listed in Goldberg (2005).
The first goal in multi-objective optimization is achieved by a proper fitness as-signment strategy and a careful reproduction operator. Diversity in the pareto-set can be obtained by designing a suitable selection operator. Preserving the elitism during generations shall be carefully considered in evolutionary algorithms. Elite-preserving operators, as Deb (2001) names them, are introduced to directly carry over the elit solutions to the next generation.
Coello (2007) presents the basic concepts and approaches of multi-objective opti-mization evolutionary algorithms. The book further explores some hybrid methods
and introduces the test functions and there analysis. Various applications of multi-objective evolutionary algorithms (MOEA) are also discussed in the book. Deb (2001) is another comprehensive source of different MOEAs. The book divides the evolutionary algorithms into non-elitist and elitist algorithms.
2.3.6 MOEA Techniques
All researchers are agreed upon that the invention of first MOEA is devoted to David Schaffer with his Vector Evaluation Genetic Algorithm (VEGA) in the mid-1980s, aimed at solving optimization problems in machine learning.
Deb (2001) and Coello (2007) both name various MOEAs which shows the differ-ence in the frame work and their operators as follows:
• Vector Evaluated GA (VEGA)
• Vector Optimized Evolution Strategy (VOES)
• Weight Based GA (WBGA)
• Multiple Objective GA (MOGA)
• Niched Pareto GA (NPGA, NPGA2)
• Non-dominated Sorting GA (NSGA,NSGA-II)
• Distance-Based Pareto GA (DPGA)
• Thermodynamical GA (TDGA)
• Strength Pareto Evolutionary Algorithm (SPEA, SPEA2)
• Multi-Objective Messy GA (MOMGA-I, II, III)
• Pareto Archived Evolution Strategy (PAES)
• Pareto Enveloped Based Selection Algorithm (PESA, PESA II)
• Micro GA-MOEA (µGA, µGA2)
Coello (2007) describes the concept of each EA along with an illustration of algo-rithm and short notes on advantages and disadvantages. At the end he summarizes all EAs in a table. While Deb (2001) devotes two complete chapter of the book
to fully define the concept and principle of each EA by step-by step description of algorithm, hand calculation, discussion on advantages and short comings, calcu-lating the computational complexity and simucalcu-lating an identical test problem for all algorithms.
2.3.7 Comparison of MOEAs
Since there exist several MOEAs, a question of which algorithm has the best performance is a common question among scientist and researchers. In order to settle to an answer several test problems has been designed and various amount of researches is carried out. In Deb’s book, a few significant studies on comparison of EAs are discussed. (Deb, 2001)
Konak et al. (2006) demonstrates the advantages and disadvantages of most well-known EAs in a table.
However the most representative, discussed and compared evolutionary algorithms are Non-dominated Sorting GA (NSGA-II) (Deb et al., 2002), Strength Pareto Evolutionary Algorithm (SPEA, SPEA2) (Zitzler and Thiele, 1998; Zitzler et al., 2001), Pareto archived Evolution Strategy (PAES)(Knowles, 1999, 2000) , and Pareto Enveloped Based Selection Algorithm (PESA, PESA II) (Corne and Knowles, 2000; Corne et al., 2001).
Extensive comparison studies and numerical simulation on various test problems shows a better overall behavior of NSGA-II and SPEA2 compared to other al-gorithms. In cases where more than two objectives are present SPEA2 seems to indicate some advantages over NSGA-II. Strength Pareto Evolutionary Algorithm (SPEA2) is comprehensively described in next section. Also SPEA2 is coded and implemented on a number of test functions.
2.3.8 SPEA2: Improved Strength Pareto Evolutionary Algorithm
Zitzler et al. (2001) improves the original SPEA (Zitzler and Thiele, 1998), and addresses some potential weaknesses of SPEA.
SPEA2 uses an initial population and an archive (external set). At the start, random initial and archive population with fixed sizes are generated. The fitness value of each individual in the initial population and archive is calculated per iteration. Next, all non-dominated solutions of initial and external population are
copied to the external set of next iteration (new archive). With the environmental selection procedure the size of the archive is set to a predefined limit. After wards, mating pool is filled with the solutions resulted from performing binary tournament selection on the new archive set. Finally, cross-over and mutation operators are applied to the mating pool and the new initial population is generated. If any of the stopping criteria is satisfied the non-dominated individuals in the new archive forms the pareto-optimal set.
Kim et al. (2004) adds two new mechanisms to SPEA2 in order to improve the searching ability of algorithm. The SPEA2 + algorithm, as it is named, uses a more effective crossover (Neighborhood Crossover) and new archive mechanism to diversify the solutions in both objective and variable spaces.
Kukkonen (2006) introduces a pruning method, which can be used to improve the performance of SPEA2. The idea of pruning is to reduce the size of a set of non-dominated solution to a pre-defined limit, while the maximum possible diversity is encountered.
2.3.9 Overall SPEA2 Algorithm
The overall algorithm of SPEA2 is as follows:(Zitzler et al., 2001)
Algorithm (SPEA2 Main Loop)
Input : N (population size) N (archive size)
T (maximum number of generations) Output: A (non-dominated set)
Step 1: Initialization: an initial population P0 and archive (external set) P0 is
generated. Set t = 0.
Step 2: Fitness assignment : Fitness values of individuals in Pt and Pt are
calcu-lated. (Fitness Assignment section)
Step 3: Environmental selection: All non-dominated individuals in Pt and Pt
shall be copied to Pt+1. If size of Pt+1 exceeds N , reduction of Pt+1 is
dominated individuals in Pt and Pt, if size of Pt+1 is less than N .
(Environ-mental Selection)
Step 4: Termination: If t > T or another stopping criterion is satisfied then, the non-dominated individuals in Pt+1 creates the output set A.
Step 5: Mating Selection: Binary tournament selection with replacement is per-formed on Pt+1 in order to fill the mating pool.
Step 6: Variation: Recombination and mutation operators shall be applied to the mating pool and the resulting population is set to Pt+1 . Increment
generation counter (t = t + 1) and go to Step 2.
Fitness Assignment
Each individual i in the archive Pt and the population Pt is assigned a strength
value S(i), representing the number of solutions it dominates:
S(i) = |{j|j ∈ Pt+ Pt∧ i j}||,
the raw fitness R(i) of an individual i is calculated:
R(i) = X
j∈Pt+Pt,j i
S(j).
The density estimation technique is adopted from the k-th nearest neighbor method (Silverman 1986), where the density at any point is a (decreasing) function of the distance to the k-th nearest data point. In SPEA2 the inverse of the distance to the k-th nearest neighbor is considered as the density measurement. The density D(i) corresponding to i is defined by:
D(i) = 1 σk i + 2 , where, k =pN + N .
and σik is the distance of solution i to the k-th nearest neighbour. Finally, the fitness of an individual i is calculated by adding D(i) to the raw fitness value R(i):
F (i) = R(i) + D(i)
Environmental Selection
The first step is to copy all non-dominated individuals, i.e., the ones with fitness value lower than one, from archive and population to the external set of the next generation:
Pt+1 = {i|i ∈ Pt+ Pt∧ F (i) < 1}.
If the size of non-dominated solutions is exactly the same as archive size (|Pt+1| =
N ) the environmental selection step is completed. Otherwise, there can be two situations:
• The archive is too small (|Pt+1| < N ): The best N − |Pt+1| dominated
individuals in the previous external set and population are copied to the new archive. This can be achieved by sorting the multi-set Pt+ Pt+1 from
lowest to highest fitness values and copy the first N − |Pt+1| individuals i
with F (i) > 1 from the sorted list to Pt+1 .
• The archive is too large (|Pt+1| > N ): In this case, an archive truncation
procedure is invoked which iteratively removes individuals from Pt+1 until
(|Pt+1| = N ) . Here, at each iteration the solution i is chosen for removal
for which i ≤dj for all j ∈ Pt+1 with:
i ≤d j :⇔ ∀ 0 < k < |Pt+1| : σjk = σ k i ∨ ∃ 0 < k < |Pt+1| : ∀ 0 < l < k : σl i = σ l j ∧ σ k i < σ k j , where σk
i denotes the distance of i to a user-predefined (k-th) nearest
neigh-bor in Pt+1 . In other words, at each stage removed solution will be the one
with the least distance to the k-th neighbor; if there is more than one solu-tion with the same distance the judgement will be upon the second smallest distance and so forth.
In appendix B, hand calculation and step by step simulation of a simple example minimization problem is fully described. This will help on better understanding of algorithm and the working principle of each step.
2.3.10 Constraint Handling
Handling constraints within MOEAs is an essential issue which must be consid-ered carefully, especially when dealing with real world engineering applications where constraints are always involved. Constraints can be in form of equality or inequality. Another classification of constraints are hard and soft constraints. A hard constraint is a must to be satisfied, while on the other hand, a soft one can be relaxed in order to accept a solution (Coello, 2007; Deb, 2000). Normally only inequality constraints are handled in MOEAs, however equality constraints can be easily transformed to inequality using:
|h(x)| − 6 0
where h(x) = 0 is the equality constraint and is very small value.
Constraints divide the decision space into two separate parts: feasible and infea-sible regions. A solution in the feainfea-sible region of search space satisfies all the constraints and it is called a feasible solution, otherwise the solution is infeasible.
The most popular and common way of handling constraints is the penalty function method. However sensitivity of penalty method to the penalty parameter is a drawback in this method.
In addition to penalty method, Jimnez et al. (1999) proposed a systematic con-straint handling procedure. Two other method which are more credited and elab-orated are the Ray-Tai- Seows constraint handling approach (T. Ray, 2001) and the Deb et al. (2002) proposed constraint handling method, which is implemented in NSGA II algorithm.
In Debs method a binary tournament selection operator is used for any two so-lutions selected from the population. Therefore in presence of constraints three scenarios will occur: 1) Both solutions are feasible; 2) One is feasible and the other is infeasible; 3) Both solutions are infeasible. In the method for each scenario fol-lowing rule is applied:
• Scenario 1) the solution with better objective function is selected (Crowded comparison).
• Scenario 2) the feasible solution will win the tournament.
• Scenario 3) the solution with less constraint violation is selected.
Deb modifies the definition of domination as solution i constraint dominates solu-tion j, if any of the following condisolu-tions is true:
1. Solution i is feasible and j not.
2. Solution i and j are both infeasible, but solution i has a smaller overall constrained violation.
3. Solution i and j are both feasible and solution i dominates solution j.
In the SPEA2 algorithm proposed for the thesis, binary tournament is applied to the archive population (Pt+1), which holds the non-dominated individuals to create
the mating pool. Afterwards the genetic operators are used to generate the child from the mating pool which is the initial population of next generation (Pt+1). In
constrained problems the modified definition of domination is implemented and the non-dominated solutions are selected according to constraint domination concept. Appendix B.1 simulates the principle and procedures of constraint domination concept on a simple problem.
2.4
Hybrid Multi-Objective Optimization Approach
In real world engineering problems there is no prior knowledge on the true global pareto-front. Although Evolutionary algorithms have shown a good convergence in benchmarks, hybrid methods have been proposed to ensure the convergence of an algorithm to the true pareto-front. Several hybridization techniques (combining an MOEA with other methods) are discussed in literatures.
Coello (2007) comprehensively deliberate the use of local search and co-evolutionary techniques as a hybrid method in a complete chapter of his book. He specifies local search decision space approaches such as depth-first search (hill-climbing), simu-lated annealing and Tabu search for consideration in hybridization.
Deb (2001) also argues the use of local search techniques with an MOEA. Accord-ing to Goldberg, the best way to achieve convergence to the exact pareto-front is implementing the local search techniques on the solutions obtained from an EA. However Deb proposes two other ways to use local search techniques; 1) during EA generations, 2) at the end of an EA run.
Here a new method of hybridization is introduced and tested on benchmarks to investigate the performance of the technique. A combination of single and multi-objective optimization evolutionary algorithms discussed in previous subsections are applied to obtain the global optimal solutions. The archive population in SPEA2, which holds the non-dominated solutions of each generation, is created using the single-objective genetic algorithm optimization method introduced in earlier sections called G3 algorithm.
First, the objectives are transformed to a single objective function by using the weighted sum method. A number of random weights equal to size of population are multiplied to each objective to scalarize the objective function. Then, every scalarized function is optimized with G3 single objective GA. After finding the optimal solution for each weighted function, the required initial population is obtained. Finally, the multi-objective algorithm (SPEA2) is used to optimize the function. Therefore, the hybridizing technique is applied before EA generations to create the required initial archive population.
3
Implementation
All the algorithms in the thesis are coded with MATLAB. Several benchmarks are encompassed and solved with the coded algorithms to ensure the accuracy and efficiency of algorithm.
3.1
Single Objective
3.1.1 Unconstrained Test Functions
Sphere Function fSphere(x) = n X i=1 x2i
has a global minimum of 0 at x∗ = (0, 0)T.
Ellipsoidal Function
The behaviour of the algorithms for a poorly scaled objective function is discussed using the following objective function:
fEllipsoid(x) = ax21+ n
X
i=1
x2i
where a is a positive parameter. If a = 1, the function has a valley structure. In this experiment, a = 0.01 is used. The global minimum is 0 at x∗i = 0, i = 1, ..., n.
Schwefel’s Function fSchwef el(x) = n X i=1 i X j=1 xi !2
Goldstein-Price Function
The Goldstein-Price function is given by:
fGoldstein(x1, x2) = (1 + (x1 + x2+ 1)2)(19 − 14x1+ 3x21− 14x2+ 6x1x2+ 3x22))
×(30 + (2x1− 3x2)2(18 − 32x1+ 12x12+ 48x2− 36x1x2+ 27x22))
has a global minimum of 3 at x∗ = (0, −1)T. The typical search range is −2 ≤ xi ≤ 2, i = 1, 2.
Rosenbrock Function
This function is used to discuss the behaviour of the algorithms for functions having complex non-separable structure, such as a curved, deep valley, given by
fRosenbrock(x) = n
X
i=2
(100(x21 − xi)2+ (1 − xi)2),
and has a global minimum of 0 at x∗i = 1. The typical search range is −5.12 ≤ xi ≤ 5.12, i = 1, ..., n.
Colville Function
The Colville function is defined as
fColville(x1, x2, x3, x4) = 100(x2− x21) 2 + (1 − x 1)2+ 90(x4− x23) 2+ (1 − x 3)2 +10.1((x2− 1)2+ (x4− 1)2) + 19.8(x2− 1)(x4− 1).
The search range is −10 ≤ xi ≤ 10 and the global minimum of fColville∗ = 0 at
x∗i = 1, i = 1, ..., 4.
Considering the results of systematic studies on parameters of G3 algorithm (Deb, 2004), in all above cases, a population size of N = 100, a parent size µ = 3, number of offspring λ = 2 and r = 2 (Step 2 and 3) are used. For the PCX different values of ση and σζ is implemented.
In addition, for Spherical, Ellipsoidal, Schwefel and Rosenbrock functions two cases are considered for initial population:
• Normal Case: the distribution of initial population is surrounding the optimal solution. The population is generated by uniform random numbers in the region below:
−1 < xi < 1, i = 1, ..., n
• Offset Case: the distribution of initial population is faraway from the op-timal solution.The population is generated by uniform random numbers in the region below:
−10 < xi < −5, i = 1, ..., n
3.1.2 Constrained Test Functions
Three constrained optimization problems, two of them real engineering problems, is used to evaluate the performance of the G3 algorithm and selected constraint handling method.
Test Function 1
This test problem is a two dimensional constrained optimization problem:
Minimize f (x) = (x21+ x2− 11)2+ (x1+ x22− 7) 2,
Subject to g1(x) = 4.84 − (x1− 0.05)2− (x2− 2.5)2 ≥ 0,
g2(x) = x21+ (x2− 2.5)2− 4.84 ≥ 0,
0 ≤ x1, x2 ≤ 6.
The optimal solution is x∗ = (2.246826, 2.381865) with a function value equal to f∗ = 13.59085.
Welded Beam Design
The objective is to minimize the cost of the welded beam subject to the constraints on shear stress (τ ), bending stress in the beam (σ), bucking load on the bar (Pc),
end deflection of the beam (δ), and side constraints. The problem has four design variables h(x1), l(x2), t(x3), b(x4) and five inequality constraints as follows:
Minimize f (x) = 1.1047x21x2+ 0.04811x3x4(14.0 + x2), Subject to g1(x) = τM AX − τ (x) ≥ 0, g2(x) = σM AX − σ(x) ≥ 0, g3(x) = x4− x1 ≥ 0, g4(x) = Pc(x) − P ≥ 0, g5(x) = δM AX − δ(x) ≥ 0, 0.1 ≤ x1, x4 ≤ 2, 0.1 ≤ x2, x3 ≤ 10 where τ (x) = r ((τ0(x))2+ (τ”(x))2+ x 2τ 0 (x)τ”(x))/ q 0.25 [x2 2+ (x1+ x3)2], σ(x) = 6P L x2 3x4 , δ(x) = 4P L 3 Ex2 3x4 , Pc(x) = 4.013E q x2 3x64 36 L2 " 1 − x3 2L r E 4G # , where τ0(x) = √ P 2x1x2 , τ”(x) = M R J 0 , M = PhL +x2 2 i , R = r x2 2 4 + ( x1+ x3 2 ) 2.
P = 6000 lb, L = 14 in, E = 30 × 106psi, G = 12 × 106psi, τM AX = 13600 psi, σM AX = 30000 psi, δM AX = 0.25 in.
The optimized solution reported in literature (Deb, 1991) is x = (0.2489, 6.1730, 8.1789, 0.2533) with f = 2.43 using binary GA.
Figure 3: The welded beam problem.
Minimization of the Weight of a Tension/Compression Spring
The problem consists of minimizing the weight of a tension/compression spring subject to constraints on minimum deflection, shear stress, surge frequency, limits on outside diameter and on design variables. The design variables are the mean coil diameter D(x2), the wire diameter d(x1) and the number of active coils N(x3).
The problem can be expressed as follows:
Minimize f (x) = x21x2(x3+ 2), Subject to g1(x) = (x32x3) (71785x4 1) − 1, g2(x) = 1 − 4 ∗ x2 2− x1∗ x2 12566 ∗ x3 1∗ (x2− x1) − 1 5108 ∗ x2 1 , g3(x) = 140.45 ∗ x1 xx3∗ x22 − 1, g4(x) = 1 − x1) + x2 1.5 − 1 , 0.05 ≤ x1 ≤ 2, 0.25 ≤ x2 ≤ 1.3, 2 ≤ x2 ≤ 15.
The best optimal solution obtained by using static penalty function is f∗ = 0.012729 and the lowest optimal found by (Mezura-Montes et al., 2003) is f∗ = 0.012688 using an approach based on Evolution Strategy.
Figure 4: Tension/compression string problem.
3.2
Multi objective
Similar to single objective, two sets of test functions, one for unconstrained and the other for constrained, are utilized to assess the performance of SPEA2 and the proposed constrained handling method. The Algorithm is coded with MATLAB.
3.2.1 Unconstrained Test Functions
Exercise 14
A non-convex function presented in Stromberg (2011) with one variable.
Minimize f1(x) = 1 − 2x + x2, f2(x) = √ x, 0 ≤ x ≤ 1.
The pareto-front is obtained by using two different methods; 1) Single objective G3 algorithm (the function is transformed to single-objective by weighted sum method), 2) Multi-objective SPEA2 algorithm.
Kursawe’s Test Function
Kursawe (1991) used a complicated two-objective, three variable function with a non-convex and disconnected pareto-optimal set.
KUR: Minimize F = (f1(x), f2(x)) , f1(x) = Pn−1 i=1 −10e−0.2 √ x2 i+x2i+1 , f2(x) =Pni=1 |xi|0.8+ 5 sin3i . −5 ≤ xi ≤ 5, i = 1, 2, 3.
Deb (2001) illustrates the pareto-front of KUR function in figure 201. Three distinct disconnected regions create the pareto-front of the problem. Also figure 202 of the same book shows the pareto-optimal solutions in decision space.
Zitzler-Deb-Thiele’s Test Functions
Zitzler et al. (2000) introduced six set of multi-objective problems (ZDT1 to ZDT6) which are based on a unique structure with different level of difficulties. The functions have two objective with the aim of minimization:
ZDT: Minimize F = (f1(x), f2(x)) , f1(x), f2(x) = g(x)h(f1(x), g(x)).
In the thesis five of the six test problems (ZDT1, ZDT2, ZDT3, ZDT4 and ZDT6) are implemented with SPEA2 algorithm.
ZDT1 Test Function
ZDT1 is a function with 30 variables and convex pareto-optimal set as follows:
ZDT1: Minimize F = (f1(x), f2(x)) , f1(x) = x1, f2(x) = g(x) 1 −pf1/g(x) , g(x) = 1 + 9Pn i=2 xi n − 1.
All the variables are limited between 0 and 1. Figure 213 of Deb (2001) shows the search space and pareto-front in objective space. This is the easiest among all ZDT’s and the only difficulty is the large number of variables.
ZDT2 Test Function
Another 30-variable test function with a non-convex pareto-front:
ZDT2: Minimize F = (f1(x), f2(x)) , f1(x) = x1, f2(x) = g(x) 1 − (x1/g(x))2 , g(x) = 1 + 9 n − 1 Pn i=2xi.
The range for all the variables is [0, 1]. Pareto-front and the search region in objective space is shown in figure 214 of Deb (2001). Non-convexity of pareto-optimal set is the only difficulty of this problem.
ZDT3 Test Function
ZDT3 problem with 30 variables, has a number of disconnected pareto-optimal sets: ZDT3: Minimize F = (f1(x), f2(x)) , f1(x) = x1, f2(x) = g(x) 1 −pf1/g(x) − (f1/g(x)) sin (10πf1) , g(x) = 1 + 9Pn i=2 xi n − 1.
All the variables are limited within [0, 1]. Finding all the discontinuous pareto-optimal regions with a good diversity of non-dominated solutions may be difficult for an MOEA. Deb (2001) show the search space and pareto-front in figure 215.
ZDT4 Test Function
This is a 10 variable problem with a convex pareto-front:
ZDT4: Minimize F = (f1(x), f2(x)) , f1(x) = x1, f2(x) = g(x) 1 −px1/g(x) , g(x) = 1 + 10 (n − 1) +Pn i=2(x 2 i − 10 cos (4πxi)) .
All the variables except x1, which lies in the range [0, 1], are limited within −5
and 5. Large number of multiple local pareto-fronts, shown in figure 216 of Deb (2001), will create a difficult convergence to global pareto-front for an MOEA.
ZDT6 Test Function
This is a problem with 10 variables and a non-convex pareto-optimal set:
ZDT6: Minimize F = (f1(x), f2(x)) , f1(x) = 1 − exp(−4x1) sin6(6πx1), f2(x) = g(x) 1 − (f1/g(x)) 2 , g(x) = 1 + 9 Pn i=2 xi n − 1 1/4 .
All the variables lie in the range [0, 1]. Non-convexity of pareto-front, coupled with adverse density solutions across the front, may rise some difficulty in convergence. Figure 218 in Deb (2001) shows the pareto-optimal region for this problem.
3.2.2 Constrained Test Functions
The presence of constraints may cause hurdles for an MOEA to converge to the true and global pareto-front, also maintaining diversity in the non-dominated solutions may be another problem. A number of common test problems used in literatures, are presented in this section and implemented in the SPEA2 code.
Binh and Korn Test Function
Binh and Korn (1997) introduced a problem with two-variable as follows:
BNH: Minimize f1(x) = 4x21+ 4x22, Minimize f2(x) = (x1− 5)2+ (x2− 5)2, subject to C1(x) = (x1− 5)2+ x22 ≤ 25, C2(x) = (x1− 8)2+ (x2+ 3)2 ≥ 7.7, 0 ≤ x1 ≤ 5, 0 ≤ x2 ≤ 3.
Deb (2001) illustrates the decision variable and objective space of the problem in figures 219 and 220. In BNH problem, constraints will not add any difficulty to the unconstrained problem.
Osysczka and Kundu Test Function
Osyczka and Kundu (1995) used the following six variable test function:
OSY: Minimize f1(x) = −[25(x1− 2)2+ (x2− 2)2+ (x3− 2)2+ (x4− 2)2 +(x5− 2)2], Minimize f2(x) = x21+ x22+ x23+ x24+ x25+ x26, subject to C1(x) = x1+ x2− 2 ≥ 0, C2(x) = 6 − x1− x2 ≥ 0, C3(x) = 2 + x1− x2 ≥ 0, C4(x) = 2 − x1+ 3x2 ≥ 0, C5(x) = 4 − (x3− 3)2− x4 ≥ 0, C6(x) = (x5− 3)2+ x6− 4 ≥ 0, 0 ≤ x1, x2, x6 ≤ 10, 1 ≤ x3, x5 ≤ 5, 0 ≤ x4 ≤ 6.
The pareto-front as shown in figure 221 of Deb (2001), is a line connecting some parts of five different region. Since the algorithm should maintain the solutions within intersections of constraint boundaries, this is a difficult problem to solve.
Srinivas and Deb Test Function
Srinivas and Deb (1994) suggested the following problem:
SRN: Minimize f1(x) = 2 + (x1− 2)2+ (x2− 1)2, Minimize f2(x) = 9x1− (x2− 1)2, subject to C1(x) = x21 + x22 ≤ 225, C2(x) = x1 − 3x2+ 10 ≤ 0, −20 ≤ x1, x2 ≤ 20.
Since the constraints eliminate some parts of the original pareto-front, difficulties may arise in solving the problem. Figures 222 and 223 in Deb (2001) shows the corresponding pareto-front of feasible decision variable and objective space.
Tanaka Test Function
Tanaka and Watanabe (1995) proposed a two variable test function as follows:
TNK: Minimize f1(x) = x1, Minimize f2(x) = x2, subject to C1(x) = x21+ x22− 1 − 0.1 cos 16 arctanx1 x2 ≥ 0, C2(x) = (x1− 0.5)2+ (x2− 0.5)2 ≤ 0.5, 0 ≤ x1, x2 ≤ π.
The pareto solutions lie on a surface which is non-linear. Therefore optimization algorithms may face some difficulties in finding a diverse set of feasible pareto solutions. A figure of feasible decision variable spaces for the problem can be seen in Deb (2001). (Figure 224)
3.3
Hybrid Approach
In order to test the performance of the proposed hybrid method, the archive pop-ulation generated from the hybrid technique (weighted sum single objective G3 algorithm) is compared with the random archive population for different test prob-lems. It is obvious that the test functions with non-convex pareto-fronts are not suitable for the technique, since the weighted sum method is used to transform the objectives into one objective. Therefore, the random and hybrid archive pop-ulation are plotted in objective space for ZDT1 and ZDT3 test functions. Also despite the non-convex property of ZDT6 problem, comparison of archive popu-lation has also been applied to this problem to study the performance of hybrid method on non-convex benchmarks.
4
Test Results
4.1
Single Objective
4.1.1 Unconstrained Functions
To examine the behaviour of algorithm and code, evaluation in two and multi-dimensional search space is carried out for some of the test functions as blow:
• Spherical: n = 2, 4
• Ellipsoidal: n = 2, 4
• Schwefel: n = 2, 10, 15
• Rosenbrock: n = 2, 5
Goldstein function is by default in two dimensional search space and Colville is a four variable function.
Table 1: Results of unconstrained test functions, single objective.
Function Initial Number of Number of Variance from ση,
Population Variable Evaluation Global Optimum σζ
Sphere Normal 2 98 7.32 × 10−6 0.1 Sphere Normal 2 67 5.47 × 10−6 0.4 Sphere Normal 4 279 3.54 × 10−6 0.4 Sphere Offset 2 261 3.65 × 10−6 0.1 Ellipsoidal Normal 2 99 3.75 × 10−6 0.1 Ellipsoidal Offset 2 405 3.65 × 10−6 0.1 Ellipsoidal Offset 4 553 8.10 × 10−6 0.4 Schwefel Normal 2 553 8.10 × 10−6 0.1 Schwefel Offset 10 553 8.10 × 10−6 0.3 Schwefel Offset 15 6434 7.76 × 10−6 0.3 Rosenbrock Normal 2 146 3.55 × 10−6 0.1 Rosenbrock Offset 5 2200 3.9308 × 10−6 0.5 Rosenbrock Offset 5 5095 3.55 × 10−6 0.9 Goldstein - 2 211 9.46 × 10−6 0.1 Colville - 4 no result no result 0.1 Colville - 4 1055 9.46 × 10−6 0.3 Colville - 4 1468 8.09 × 10−6 0.6
The experiment for each function runs until the best objective function of the pop-ulation reaches a minimum difference of 10−5 from the optimal solution. Number of Generation (evaluation) and best fitness are shown in table 1.
The result shows acceptable behaviour of algorithm for two dimensional search space with ση = 0.1 and σζ = 0.1 , but when the number of variable increases
or functions are more complex such as Rosenbrock, the algorithm converges in local optima or the global optima is obtained with high number of generations. Therefore by increasing the variance of zero-mean normally distributed variables in PCX operator better results are obtained as it can be seen in table of results the Schwefel’s function with 15 variable has reached the required variance from global optimal. Convergence of the best individual obtained from some of the test functions during generations can be seen in the figures 5 and 6.
Figure 5: Convergence of Variance from optimal solution for Schwefel’s function with 15 variables and optimal f∗ = 0
Figure 6: Convergence of Variance from optimal solution for Rosenbrock function with 5 variables and optimal f∗ = 0
4.1.2 Constrained Functions
The Penalty method used for constrained test functions shows a good behaviour. All three functions reached the optimal solution reported in literature, in ad-dition the optimal solution found by the algorithm in this thesis with related constraint handling method is better than some other approaches used in liter-ature. In the first test function the optimal solution of 13.590842 is found at x∗ = (2.246818, 2.381735) which is better than the optimal found at literature (Deb, 2000) with the value f∗ = 13.59085. Furthermore, figures 7, 8 and 9 illus-trate the convergence of results for the three constrained test functions.
Tables 2 and 3 compare the results of different methods for welded beam de-sign problem and the minimization of the weight of a tension/compression spring, which shows the out-performance of the approach used here to some methods.
Figure 7: Convergence of objective function for Test Function 1 with obtained optimal of f∗ = 13.590842
Figure 8: Convergence of objective function for welded beam problem with ob-tained optimal of f∗ = 1.834756
Table 2: Comparison of the results of different methods for welded beam design problem.
Method f∗(x) This Thesis 1.83475678 Coello (self-adaptive penalty approach) 1.74830941 Arora (constraint correction at constant cost) 2.43311600 He and Wang (CPSO) 1.728024 Ragsdell and Phillips (Geometric programming) 2.385937 Deb (GA) 2.433116 Coello and Montes (feasibility-based tournament selection) 1.728226 Ebehart (modified PSO) 1.72485512
Figure 9: Convergence of objective function for Tension/Compression Spring with obtained optimal of f∗ = 0.012710175
Table 3: Comparison of the results of different methods for the minimization of the weight of a tension/compression spring.
Method f∗(x) This Thesis 0.012710175 Coello (self-adaptive penalty approach) 0.01270478 Arora (constraint correction at constant cost) 0.12730274 He and Wang (CPSO) 0.0126747 Belegundu (numerical optimization technique) 0.0128334 Coello and Montes (feasibility-based tournament selection) 0.0126810 Ebehart (modified PSO) 0.01266614
4.2
Multi-Objective
SPEA2 algorithm coded with MATLAB is used to solve multi-objective test func-tions. Simulated Binary Crossover (SBX) and Polynomial Mutation (Deb, 2001) are implemented in the step 6 (Variation) of SPEA2 algorithm as recombination and mutation operators.
4.2.1 Unconstrained Functions
In exercise 14 test function, a weight vector with 20 weight factors within the range [0, 1] with step length of 0.05 is used to create 20 single objective functions and each function is optimized separately to obtain the pareto-front shown in figure 10 . In multi-objective method, the initial and archive set with population of 30 individuals after 100 generations result in the pareto-front (figure 11).
Figure 10: Pareto-front of Exercise 14 using weighted single objective algorithm
Since the function is non-convex the pareto-front obtained from single objective method does not cover the non-convex parts of pareto-optimal set. In the other hand, the pareto-front of multi-objective method clearly illustrates all parts of pareto-optimal region. Thus, an important drawback of single-objective weighted sum method for solving multi-objective optimization problems is the weakness in non-convex problems.
Table 4 shows the defined parameters of SPEA2 algorithm, such as size of initial and archive population and number of generations, for the other test functions from Kursawe to ZDT6.
Table 4: Pre-defined parameters of SPEA2 algorithm for unconstrained multi-objective test functions
Test Function Initial Population Archive Population Number of Size (N ) Size (N ) Generations
Kursawe 50 50 100 ZDT1 50 50 400 ZDT2 50 50 400 ZDT3 50 50 250 ZDT4 100 100 250 ZDT6 100 100 250
Figures 12, 13, 14, 15, 16 and 17 shows the non-dominated solutions obtained from SPEA2 algorithm for problems Kursawe, ZDT1, ZDT2, ZDT3, ZDT4 and ZDT6. Comparing the figures to the true pareto-fronts illustrated in literature (Deb, 2001), confirms the excellent performance of SPEA2 algorithm in finding the pareto-front of problems with upto 30 variables.
