Abstract

Multiobjective optimization involves minimizing or maximizing multiple objective functions subject to a set of constraints. In this study, a novel constrained multiobjective biogeography optimization algorithm (CMBOA) is proposed. It is the first biogeography optimization algorithm for constrained multiobjective optimization. In CMBOA, a disturbance migration operator is designed to generate diverse feasible individuals in order to promote the diversity of individuals on Pareto front. Infeasible individuals nearby feasible region are evolved to feasibility by recombining with their nearest nondominated feasible individuals. The convergence of CMBOA is proved by using probability theory. The performance of CMBOA is evaluated on a set of 6 benchmark problems and experimental results show that the CMBOA performs better than or similar to the classical NSGA-II and IS-MOEA.

1. Introduction

In the domain of science and engineering, most of the problems are attributed to constrained multiobjective optimization problems (CMOPs), which need to optimize multiple conflicting objectives subject to various inequality and equality constraints. So the algorithms of solving CMOPs have to search the set of nondominated feasible solutions fulfilling all constraints. It is desirable that those gained solutions can approximate the true Pareto front with better diversity and even distribution. Evolutionary algorithms (EAs) are population-based search algorithms and can find multiple optimal solutions in one single run, and they are suitable to solve multiobjective problems (MOPs). But for the specific application of solving CMOPs, we find that most of the existing constrained multiobjective EAs (MOEAs) cannot effectively exploit the population because their obtained convergence and diversity are not acceptable.

Biogeography-based optimization (BBO) algorithm is a population-based search algorithm [1, 2], which had been applied to solve single objective optimization problems (SOPs) and some engineering problems [35]. In the aspect of MOPs, Ma et al. decomposed multiobjective optimization problems into several related subproblems and used parallel BBO to optimize each subproblem [6]. We successfully improved BBO for MOPs, which had proved that the migration strategy of BBO is effective for solving MOPs [7, 8]. In view of good population exploiting ability of BBO, in this study, we propose a novel constrained multiobjective biogeography optimization algorithm (CMBOA) for the first time and analyze its convergence by the probability theory.

The proposed CMBOA includes the following features. First, the individuals are classified into the feasible and infeasible ones based on their constraint violation. Second, feasible individuals are evaluated by combining their objective functions value and crowded distance. Infeasible individuals are evaluated by combining their constraint violation and Euclidean distance from the nearest nondominated feasible individual. Third, a new migration operator with additional disturbance is designed to generate diverse feasible solutions. And infeasible solutions nearby feasible regions are recombined with their nearest nondominated feasible solutions to evolve towards feasibility.

In rest of the paper, reviews of multiobjective evolutionary algorithms (MOEAs) for CMOPs are given in Section 2, and basic conception of CMOPs, the review of CMOPs, and the BBO are briefly introduced. The CMBOA is proposed in Section 3. In Section 4, compared with the classical algorithms on benchmark CMOPs, simulation results on CMBOA are analyzed and discussed. At last, conclusions are drawn in Section 5.

2.1. Problem Statement

The aim of the constrained multiobjective optimization problems (constrained MOPs) is to find multiple nondominated solutions under constraints. If these nondominated solutions are uniformly distributed and widely spread along the Pareto front, their quality is better. Without the loss of generality, we consider the minimization of CMOPs, which can be defined as follows: where is a decision vector with decision variables. is an objective vector with objects. Each dimension variable of the decision space is bounded by its upper bound and lower bound . and are inequality constraints and equality constraint, respectively. The equality constraints generally should be transformed into inequality form and combined with other inequality constraints as follows: where , , , and is a tolerance parameter for the equality constraint. In this paper, only CMOPs with inequality constraints are considered. Constraint violation function of a solution is defined as follows: where . If , then is an infeasible solution; otherwise, it is a feasible solution. By the degree of constraint violation, infeasible solutions can be compared with another one. For feasible solutions, Pareto domination is defined as follows, which is applied to evaluate their fitness.

Definition 1 (Pareto domination solution). Let , and a solution vector is said to dominate a solution and is denoted by if

2.2. Reviews of CMOEA

Most of MOEAs are proposed for solving unconstrained multiobjective optimization [9]. According to different constraint handling methods adopted in MOEAs, the existing constrained multiobjective evolutionary algorithms (CMOEAs) can be categorized into five main groups.

The first group adopts the constraint handling techniques applied for single objective constraint optimization [1012]. Geng et al. proposed a constrained evolutionary multiobjective optimization with infeasible elitists and stochastic ranking selection (IS-MOEA) [10]. The algorithm conserves infeasible elitists that acts as bridges connecting disconnected feasible regions, and stochastic ranking is adopted to balance objectives and constraints. IS-MOEA especially obtains improvement on the problems with two or more disconnected feasible regions.

The second group uses the basic mechanism of MOEAs and handles constraints by optimizing them as additional objectives. Mezura-Montes and Coello put forward a naïve method to solve CMOPs by ignoring infeasible solutions [13]. The algorithm is easy to implement, but when feasible regions are small and surrounded by infeasible solutions, it is difficult to find feasible solutions.

The third group is based on ranking of priority of the feasible and infeasible solutions [1416]. Fonseca and Fleming proposed a unified approach for multiobjective optimization and multiple constraint handling [14]. Their algorithm handled constraints by assigning high priority to constraints and low priority to objective functions, when focusing on search of feasible solutions. Srinivas and Deb proposed a constrained multiobjective algorithm, in which constrained dominating relation of individuals is defined [16]. In this algorithm, all feasible solutions dominate all infeasible ones. Feasible solutions are sorted by their Pareto dominating relations and infeasible solutions are sorted based on their constraint violation. The algorithm can gain better performance but unfortunately it ignored the contribution of infeasible solutions to the Pareto front.

The forth group uses repair scheme to reproduce feasible solutions or less violated solutions from the original infeasible solutions [1719]. Jimenez et al. proposed the evolutionary algorithm of nondominated sorting with radial slots (ENORA) [17], which employs the min.-max. formulation for constraint handling. Feasible individuals evolve toward optimality, while infeasible individuals evolve toward feasibility. Harada et al. proposed Pareto descent repair (PDR) operator that searches feasible solutions out of infeasible individuals in the constraint function space [19].

The fifth group designs new mechanisms to evolve feasible solutions towards Pareto front and evolve the infeasible solutions towards feasible regions [2023]. Ray et al. suggested using three different nondominated rankings of the population [20]. The first ranking is performed by using the objective function values; the second is performed by using different constraints; and the last ranking is based on the combination of all objective functions and constraints. Depending on these rankings, the algorithm performs according to the predefined rules. Chafekar et al. proposed two novel approaches for solving constrained multiobjective optimization problems [21]. One method called objective exchange genetic algorithm of design optimization (OEGADO) runs several GAs concurrently with each GA optimizing one objective and exchanging information about its objective with others. Another called objective switching genetic algorithm for design optimization (OSGADO) runs each objective sequentially with a common population for all objectives. Deb proposed GA's population-based approach that does not require any penalty parameter. Once sufficient feasible solutions are found, a niching method (along with a controlled mutation operator) is used to maintain diversity among feasible solutions [23].

2.3. Biogeography-Based Optimization (BBO)

Biogeography is the science of the geographical distribution of biological organisms. In BBO, each problem solution is considered as a “habitat” with habitat survival index (HSI), which is similar to the fitness of EAs to evaluate an individual. High HSI habitats share their features with low HSI habitats. The process of sharing good features among solutions is denoted as migration. BBO adopts the migration strategy to share information among solutions. Good individuals’ information can be conserved during the evolutionary process to ensure the population convergence. A mutation operator is used to generate diverse solutions to promote the diversity of the population. The detailed operations are described as follows.

Suppose that the species number of each individual is , and then its immigration rate and emigration rate can be calculated as follows [1]: where is the most species number of all habitats. and represent the maximization of immigration rate and emigration rate, respectively. In migration operator, the individuals’ immigration rate and emigration rate are used to decide whether a solution should share its feature value with the other solutions. A better solution has a higher immigration rate and a lower emigration rate. By the migration, the solutions with high emigration rate tend to share their information with those with high immigration rate. Solutions with high immigration rate accept a lot of features from solutions with high emigration rate. With the aid of migration, BBO shows good exploitation ability in the search space.

Consider that species number change with species migrating; the probability that the habitat contains exactly species can be calculated using the following differential equation: Then the mutation rate is defined as [1] where is a predefined parameter, is calculated according to (6), and . The mutation operator is implemented based on . A solution with low probability is likely to mutate other solutions. Conversely, some solutions with high have very little chance to mutate. By the mutation operator, the diverse solutions are produced. The detailed operator on migration and mutation can refer to [1].

3. The Proposed Constrained Multiobjective Biogeography-Based Optimization Algorithm

3.1. CMBOA Description

In CMBOA, infeasible solutions recombine with nondominated feasible individuals and evolve towards feasibility. Firstly, the initial population is produced stochastically, and then the population is classified into the feasible and infeasible ones based on each individual’s constraint violation. Secondly, depending on whether the feasible population is empty or not, infeasible population will adopt two types of operators. If feasible population is empty, infeasible population will implement differential evolution operator until feasible individuals present; otherwise, infeasible solutions nearby feasible regions recombine with their nearest nondominated feasible solutions to obtain feasibility. Diverse nondominated feasible solutions are generated from feasible individuals by applying the novel migration operator. With the increasing of nondominated feasible solutions, update operator is used to limit their number and ensure their even distribution. Both the feasible and infeasible solutions are combined in an external archive. The proposed CMBOA is described as Algorithm 1.

Step  1: Parameter setting: population size , the size of feasible elitist archive , the size of
  infeasible elitist archive , maximum generation . Generate an initial population ,
  set the iterative generation , and archive ,
Step  2: Update of the archive
Step  2.1: Divide the population into the feasible set and the infeasible set
    based on their constraint violation.
Step  2.2: Select individuals with small fitness from by individuals’ fitness sorting to
    update feasible archive , and select individuals with small constraint violation
    from to update infeasible archive .
Step  2.3: Combine and to gain archive set , ;
    If is satisfied, output and the algorithm stops; otherwise go the next step.
Step  3: Generate the offspring population
Step  3.1: If , then , and perform the differential evolution operator on
    infeasible population to obtain . Otherwise, go to Step  3.2
Step  3.2: Implement selection operation on to gain the breeding pool , and then
    execute migration operation on to generate ;
Step  3.3: Implement the crossover and mutation on infeasible population to generate .
Step  3.4: Combine and to obtain the offspring population ,
   
Step  4: let and return Step  2.

The procedure of CMBOA is described as follows.

Step 1 (initialization). Initialize the iterative number ; the size of feasible elitist and infeasible elitist archive are and , respectively. Generate randomly the initial population with individuals; that is , the external archive .

Step 2. Update the external archive.
Step 2.1. Divide the combination populations into the feasible and infeasible ones.
Computing the constraint violation of the individuals in according to (3), we have
Depending on whether the value of is zero or not, the population is divided into the feasible subpopulation : and the infeasible subpopulation : Note that .
Step 2.2 (elitist feasible and infeasible archive). According to Definition 1, identify nondominated individuals of to form the temporary set :
If the size of is smaller than the predefined size , let . Otherwise, the crowding distance of individual , is computed as follows [24]: where denotes the kth objective function value of individual , . According to the sequencing of crowding distance, select largest crowding distance individuals from to form the elitist feasible archive :
For the infeasible population , if its size is smaller than the predefined size , then keeps invariable. Otherwise, the fitness of its individual is calculated as follows: where is the proportion of nondominated feasible individuals in current population, is constraint violation of the individual , and denotes its Euclidean distance away from the nearest nondominated feasible solution. And then the proceeding individuals with small fitness from are conserved in elitist infeasible archive .
Step 2.3 (formation of the archive). Combine elitist feasible archive and elitist infeasible archive to gain the archive :
If is satisfied, export as the output of the algorithm and the algorithm stops; otherwise, and go to Step 3.

Step 3. Generate the offspring population.
Step 3.1 (operation on feasible solutions). In CMBOA, when there are no feasible solutions in the current population, that is, , we use the mutation operator of differential evolution to produce feasible individuals [25]. That is, three individuals , , and are selected randomly from and the mutation operator is performed as (16) until feasible solutions set is not empty: where is a mutation constant and is a random number in the region . Otherwise, go to the next step.
Step 3.2 (selection operation). For feasible population , in order to ensure its convergence and even distribution, we define the fitness value of each individual by combining the nondominated rank and crowed distance of each individual , as where and denote the individual ’s crowed distance and nondominated rank, respectively, and is defined in (14). By this fitness, when the number of nondominated feasible solutions is small, individuals with lower ranks have high fitness so that they have more chance to be selected. With the number of nondominated feasible solutions increasing, more individuals with large crowded distance are selected with high probability.
Perform tournament selection operator on to form the breeding pool :
Step 3.3 (migration operation). The original migration operator of BBO has good exploitation ability of the population, but it is designed for the integer encoded individuals and single optimization problem. For continuous MOPs, the migration operator cannot ensure to produce the diverse solutions. So we propose a new migration operator. During the process of species migration, an individual is often affected by the other individuals. So we introduce a disturbance term in the migration operation to promote the diversity of the population. The detail process is shown in Algorithm 2.

For to
Randomly select two individuals , from the population
Select individual based on its immigration rate
For to n
If rand < then
Else
End if
End for
End for

In Algorithm 2, the disturbance factor is defined as where is the th variable of the individual , denotes the maximum iteration number, and is the number of iteration at current generation. The amplitude of disturbance factor decreases constantly with the increasing of generation . At the beginning, large disturbance makes the population explore a wide region in decision space. Diverse solutions will be generated to promote the diversity of population because of difference of . At the end, a small disturbance is used to exploit effectively the local regions to guarantee its convergence.

The migration operator on the population is defined as Step 3.4 (crossover and mutation operation on the infeasible population). It had been noticed that infeasible solutions can contribute to the diversity of solutions on the Pareto front. When feasible solutions exist in the current population, an individual is selected randomly from and recombined with the nearest individual of . The crossover operator is described as where is a recombination parameter in the region . By the operator, infeasible individuals nearby the feasible region will approximate the feasibility.

The above crossover operation on is

Step 3.5. Combine and to obtain the offspring population ; namely, .

Step 4. If the stopping criteria is not satisfied, and return to Step 2.

3.2. Time Complexity Analysis of CMBOA

The objectives of optimization problem are , the size of population is , the size of feasible archive is , the size of infeasible archive is , and the maximum of iterative times is . Time complexity for computation of constraint violation is . For migration and mutation operators on feasible individuals, its time complexity is , while time complexity for crossover operator on infeasible individuals is ; time complexity for updating of feasible archive is , updating of infeasible archive is , and then the worst time complexity of CMBOA is

3.3. Convergence Analysis of CMBOA

According to the description of CMBOA, it can be considered as an evolution Markov chain:

Let be feasible solution space, represents a state space composed of populations whose size is not more than , and denotes the ith state in state space. denotes that the population is in the state , and means the transformation probability from to . According to the description of CMBOA, we know that the series is an inhomogeneous Markov chain [26]. By using probability theory, the convergence of CMBOA is analyzed as follows.

Lemma 2. There exists , s.t. .

Proof. where , , and .

Lemma 3. There exists , s.t. .

Proof. where   and denotes the nondominated feasible solutions of the population .

Lemma 4. Let be the nondominated feasible solutions set; if and , then there exists , s.t.   .

Proof. Using K-C equation, we can obtain
Note that ; based on Lemmas 2 and 3, we derive

Theorem 5. CMBOA is weakly convergent for any initial population distribution; that is,

Proof. If , then If , the archive is applied to conserve the elitist individuals; then By (27) and (28), we can obtain Hence, CMBOA is weakly convergent for any initial population distribution.

4. Simulation Results

4.1. Experimental Setup

In order to test the validity of the proposed CMBOA, several benchmark functions with multiple features are selected including OSY [27], TNK [26], CONSTR [27], CTP1-CTP5 [28], CF1, CF2, CF4, and CF6 [29]. For OSY, its Pareto front is a concatenation of five regions and every region lies on the intersection of certain constraints; for TNK, its Pareto optimal solutions lie on a non-linear constraint surface; for CONSTR, its Pareto optimal set is concatenation of the constraint boundary and some parts of unconstrained Pareto optimal; for CTP serious functions, their Pareto optimal set is a collection of a number of discrete regions and most of solutions lie on non-linear constraint boundary. OSY, CONSTR, CTP1, CF4 and CF6 have continuous Pareto fronts, while the remaining ones have disjoint Pareto fronts (TNK, CTP2-CTP5, CF1, and CF2).

4.2. Performance Metrics

In this experiment, two performance metrics are selected to do quantitatively comparison.

Cover Metric C [30]. Suppose that U and V are two approximate Pareto optimal sets obtained by Algorithms 1 and 2: where denotes the dominated or equal relation. The value represents that all individuals in are weakly dominated by individuals in . denotes no individuals in which is weakly dominated by . Note that ; hence two directions must be considered simultaneously.

Hyper Volume HV [30]. The indicator calculates the volume covered by all nondominated solutions in the objective space. For each solution , a hypercube is constructed with a predefined reference point and the solution as the diagonal corners of the hypercube . All hypercubes are found and HV is calculated as follows:

The indicator is related to the approximation and diversity of the nondominated solution set. The higher the value of HV is, the better the diversity and approximation of solution set obtained is.

4.3. Performance of CMBOA on Benchmark Functions
4.3.1. Test on Benchmark Functions

To demonstrate the effectiveness of the proposed CMBOA for CMOPs, 12 benchmark functions are chosen to show its validity. In the experiment, each individual is described as a real vector. The parameters of CMBOA are set as follows: population size = 100, feasible elitist maximum size , infeasible elitist maximum size , the maximum immigration rate and migration rate , the termination generation = 100, is a random in the interval , and .

For all benchmark function, the Pareto fronts obtained by CMBOA are shown in Figure 1. From this figure, it can be seen that the Pareto optimal solutions obtained by CMBOA are very close to the true Pareto front for all benchmark functions. For most of benchmark functions, the solutions generated by the proposed algorithm can be distributed evenly on the true Pareto front except CF2, CF4, and CF6 because they have variable linkages.

4.3.2. Comparison with Original Migration Operator

In order to demonstrate the effectiveness of the novel migration operation, OSY and CTP4 are selected. The Pareto fronts gained by the algorithms with original and novel migration operator are shown in Figure 2, where “ ” denotes the Pareto front gained by CMBOA with the novel migration and “o” is the ones gained by the algorithm with original migration operator. From Figure 2, it can be seen that the algorithm with original migration cannot converge to the true Pareto front for OSY and CTP4, and only few solutions are produced for OSY. However, CMBOA with the novel migration operator obtains good convergence and diversity for OSY and CTP4, which demonstrates that the novel migration operator is superior to the original migration operator for OSY and CTP4.

4.3.3. Parameter Sensitivity Analysis

The disturbance parameter is not tuned elaborately but is set as (20). In this section, to investigate the robustness of the disturbance parameter , simulations with different settings are performed. Benchmark functions OSY and CTP4 are selected to test the sensitivity of . The Pareto fronts gained under different settings are shown in Figure 3. From the results, it is observed that the algorithms with different settings can all converge to the true Pareto front for OSY and CTP4, which illustrates that the disturbance parameter is capable to perform consistently and effectively for OSY and CTP4. So the disturbance parameter is reliable and robust to produce better solutions.

4.4. Comparison with Other Algorithms

To show the superior performance of the proposed algorithm, it is compared with the most popular multiobjective algorithms including NSGA-II [24] and IS-MOEA [1]. For NSGAII, the parameters are set as population size = 100, crossover probability = 0.9, mutation probability = , SBX crossover parameter = 20, polynomial mutation parameter = 20, and the termination generation = 100. For IS-MOEA, the parameters are set as population size = 100, crossover probability = 0.9, mutation probability = , SBX crossover parameter = 20, polynomial mutation parameter = 20 comparison probability = 0.45, penalty parameters , , and the termination generation = 100. For the proposed CMBOA, the parameters are set the same as the previous section. For all algorithms, 30 independent runs are conducted on each of the benchmark functions to get the statistical results in cover metric and hypervolume HV. Their distribution of simulation results is shown in Tables 16.

In Table 1, if the value of C (CMBOA, NSGA-II) is larger than that of C (NSGA-II,CMBOA), it indicates that the proposed CMBOA has better convergence than NSGA-II; otherwise, it indicates that the proposed CMBOA is inferior to NSGA-II in term of convergence. From Table 1, it can be seen that for CONSTR, CF2, CF4, and CF6, and NSGA-II is better than CMBOA in term of convergence. However, for the other eight test functions, CMBOA has better convergence than NSGA-II. Wilcoxon rank-sum test is used to examine their difference [31] and the results are shown in Table 2. The alternative hypothesis is and . If is met, the algorithms have significant difference; otherwise, they have no difference. From Table 2, we can see that, compared with NSGA-II, CMBOA is significantly superior on OSY, CTP1, CTP3, CTP4, CTP5, CF1, and CF6 in converging close to Pareto front. In order to analyze their difference in convergence, the distribution of their cover metric values is studied by Wilcoxon rank-sum test which is summarized in Table 3. In Table 3, the CMBOA is significantly superior to IS-MOEA in most benchmark functions except TNK, CTP1, and CTP2 on convergence. In Table 4, we can see that, except TNK, CMBOA is superior to IS-MOEA in convergence, while IS-MOEA is better than CMBOA in TNK.

In order to evaluate the convergence and the diversity of solutions obtained by the proposed CMBOA, statistical results of hypervolume metric are summarized in Table 5. In this table, higher hypervolume value indicates that the algorithm has better diversity. From Table 5, it can be seen that CMBOA has better diversity than the other two algorithms for almost all test functions except CONSTR and CF4. From the variance of metric HV, we can see that CMBOA has the smallest variance which indicates that it is more reliable and robust than NSGA-II and IS-MOEA in producing better solutions. In order to analyze the distribution of HV value in further, its Wilcoxon rank-sum test value is summarized in Table 6. From Table 6, we can conclude that CMBOA is superior to NSGA-II and IS-MOEA in terms of the distribution and diversity of solutions except CONSTR, CTP2, and CF4. Experiment results above show that the CMBOA has competitive performance with NSGA-II and IS-MOEA in the terms of convergence and diversity.

5. Conclusions

In this paper, we propose a new constrained multiobjective biogeography-based optimization algorithm. A new disturbance migration operator is designed to generate diverse feasible solutions. Infeasible solutions nearby the feasible region are recombined with their nearest feasible ones to change the feasibility. Theoretically, the weak convergence of CMBOA is proved by the probability theory and its time complexity is analyzed. Experimentally, CMBOA is tested on several typical benchmark functions and compared with classical NSGA-II and IS-MOEA. The statistical results show that the proposed CMBOA is highly competitive in terms of convergence and diversity. In future work, we may improve CMBOA to obtain better performance on variable linkage problems.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The work is supported by National Nature Science Foundation of China (no. 61075113), Excellent Youth Foundation of Heilongjiang Province (no. JC12012), the Fundamental Research Funds for the Central Universities (no. HEUCFZ1209), and Harbin Excellent Discipline Leader (no. 2012RFXXG073).