Research Article  Open Access
Yalin Wang, Xiaofang Chen, Weihua Gui, Chunhua Yang, Lou Caccetta, Honglei Xu, "A Hybrid Multiobjective Differential Evolution Algorithm and Its Application to the Optimization of Grinding and Classification", Journal of Applied Mathematics, vol. 2013, Article ID 841780, 15 pages, 2013. https://doi.org/10.1155/2013/841780
A Hybrid Multiobjective Differential Evolution Algorithm and Its Application to the Optimization of Grinding and Classification
Abstract
The grindingclassification is the prerequisite process for full recovery of the nonrenewable minerals with both production quality and quantity objectives concerned. Its natural formulation is a constrained multiobjective optimization problem of complex expression since the process is composed of one grinding machine and two classification machines. In this paper, a hybrid differential evolution (DE) algorithm with multipopulation is proposed. Some infeasible solutions with better performance are allowed to be saved, and they participate randomly in the evolution. In order to exploit the meaningful infeasible solutions, a functionally partitioned multipopulation mechanism is designed to find an optimal solution from all possible directions. Meanwhile, a simplex method for local search is inserted into the evolution process to enhance the searching strategy in the optimization process. Simulation results from the test of some benchmark problems indicate that the proposed algorithm tends to converge quickly and effectively to the Pareto frontier with better distribution. Finally, the proposed algorithm is applied to solve a multiobjective optimization model of a grinding and classification process. Based on the technique for order performance by similarity to ideal solution (TOPSIS), the satisfactory solution is obtained by using a decisionmaking method for multiple attributes.
1. Introduction
Grindingclassification is an important prerequisite process for most mineral processing plants. The grinding process reduces the particle size of raw ores and is usually followed by classification to separate them into different sizes. Grindingclassification operation is required to produce pulp with suitable concentration and fineness for flotation. The pulp quality will directly influence the subsequent flotation efficiency and recovery of valuable metals from tailings. In order to improve economic efficiency and energy consumption, the process optimization objectives include product quality and output yields. Under certain mineral source conditions, the objectives are decided by a series of operation variables such as the solid flow of feed ore to ball mill, the steel ball filling rate, and the flow rates of water added to the first and the second classifier recycles. To solve the optimization model of products’ output and quality in the grindingclassification process is of great significance to improve the technical and economic specifications, and it has been a continuous endeavor of the scientists and engineers [1–3].
Grindingclassification is an energyintensive process influenced by many interacting factors with mutual restraints. The goals of grindingclassification optimization problem are decided by multiple constrained input control variables of nonlinear relationships. So, the optimization model of grindingclassification operation is a complex constrained multiobjective optimization problem (CMOP). Generally, constrained multiobjective problems are so difficult to be solved that the constraint handling techniques and multiobjective optimization methods need to be combined for optimization.
Multiobjective optimization problems (MOPs), in the case of traditional optimization methods, are often handled by aggregating multiple objectives into a single scalar objective through weighting factors. MOPs have a set of equally good (nondominating) solutions instead of a single one, called a Pareto optimum which was introduced by Edgeworth in 1881 [4] and later generalized by Pareto in 1896 [5]. The practical MOPs are often implicated in series of equations, functions, or procedures with complicated constraints. Therefore, the evolutionary algorithms are attractive approaches for low requirements on mathematical expression [6]. Since the mid 1980s, there has been a growing interest in solving MOPs using evolutionary approaches [7–10]. One of the most successful evolutionary algorithms for the optimization of continuous space functions is the differential evolution (DE) [11]. DE is simple and efficiently converges to the global optimum in most cases [12, 13]. Its efficiency has been proven [14] in many application fields such as pattern recognition [15] and mechanical engineering [16].
There have been many improvements for DE to solve MOPs. Abbass [17] firstly provided a Pareto DE (PDE) algorithm for MOPs in which DE was employed to create new solutions, and only the nondominated solutions were kept as the basis for the next generation. Madavan [18] developed a Pareto differential evolution approach (PDEA) in which new solutions were created by DE and kept in an auxiliary population. Xue et al. [19] introduced multiobjective differential evolution (MODE) and used Paretobased ranking assignment and crowding distance metric, but in a different manner from PDEA. Robic and Filipi [20], also adopting Paretobased ranking assignment and crowding distance metric, developed a DE for multiobjective optimization (DEMO) with a different population update strategy and achieved good results. Huang et al. [21] extended the selfadaptive DE (SADE) to solve MOPs by a so called multiobjective selfadaptive DE (MOSADE). They further extended MOSADE by using objectivewise learning strategies [22]. Adeyemo and Otieno [23] provided multiobjective differential evolution algorithm (MDEA). In MDEA, a new solution was generated by DE variant and compared with target solution. If it dominates the target solution, then it was added to the new population; otherwise, a target solution was added.
On the other hand, singleobjective constrained optimization problems have been studied intensively in the past years [24–28]. Different constraint handling techniques have been proposed to solve constrained optimization problems. Michalewicz and Schoenauer [29] divided constraints handling methods used in evolutionary algorithms into four categories: preserving feasibility of solutions, penalty functions, separating the feasible and infeasible solutions, and hybrid methods. The differences among these methods are how to deal with the infeasible individuals throughout the search phases. Currently, the penalty function method is most widely used, and this algorithm strongly depends on the choice of the penalty parameter.
Although the multiobjective optimization and the constraint handling problem have received lots of contribution, respectively, the CMOPs are still difficult to be solved in practice. Coello and Christiansen [30] proposed a simple approach to solve CMOPs by ignoring any solution that violates any of the assigned constraints. Deb et al. [8] proposed a constrained multiobjective algorithm based on the concept of constrained domination, which is also known as superiority of the feasible solution. Woldesenbet et al. [31] introduced a constraint handling technique based on adaptive penalty functions and distance measures by extending the corresponding version for the singleobjective constrained optimization.
In the MOP of grinding and classification process, the definitions of Pareto solutions, Pareto frontier, and Pareto dominance are in consistency with the classic definitions. Clearly, the Pareto frontier is a mapping of the Paretooptimal solutions to the objective space. In the minimization sense, general constrained MOPs can be formulated as follows where is the objective vector, is a parameter vector, is the inequality constraint, and is the equality constraint. and are, respectively, the lower and upper bounds of the decision variable .
In this paper, based on the specific industrial background of continuous bauxite grindingclassification operation, a new hybrid DE algorithm is proposed to solve complex constrained multiobjective optimization problems. Firstly, a hybrid DE algorithm for MOPs with simplex method (SMDEMO) is designed to overcome the problems of global performance degradation and being trapped in local optimum. Then, for the MOPs with complicated constraints, the proposed algorithm is formed by combining SMDEMO and functional partitioned multipopulation. In this method, the construction of penalty functions is not required, and the meaningful infeasible solutions are fully utilized.
The remainder of the paper is structured as follows. Section 2 describes the SMDEMO algorithm for unconstrained cases. The proposed algorithm of multipopulation for constrained MOPs is given in Section 3 with verification of performance by benchmark testing results. Section 4 describes the model of products’ output and quality in the grindingclassification process in detail and the application of the proposed algorithm in the optimization model. Finally, the conclusions based on the present study are drawn in Section 5.
2. SMDEMO Algorithm for Unconstrained MOPs
In order to efficiently solve multiobjective optimization problem and find the approximately complete and nearoptimal Pareto frontier, we proposed a hybrid DE algorithm for unconstrained multiobjective optimization with simplex method.
The differential evolution, with initialization, crossover, and selection as in usual genetic algorithms, uses a perturbation of two members as the mutation operator to produce a new individual. The mutation operator of the DE algorithm is described as follows.
Considering each target individual , in the th generation of size , a mutant individual is defined by where indexes , , and represent mutually different integers that are different from and that are randomly generated over , and is the scaling factor.
The simplex method, proposed by Spendley, Hext, and Himsworth and later refined by Nelder and Mead (NM) [32], is a derivativefree linesearch method that is particularly designed for traditional unconstrained minimization scenarios. Clearly, NM method can be deemed as a direct linesearch method of the steepest descent kind. The ingredients of the replacement process consist of four basic operations: reflection, expansion, contraction, and shrinkage. Through these operations, the simplex can improve itself and approximate to a local optimum point sequentially. Furthermore, the simplex can vary its shape, size, and orientation to adapt itself to the local contour of the objective function.
2.1. Main Strategy of SMDEMO
The SMDEMO algorithm is improved by the following three points compared with DE.
2.1.1. Modified Selection Operation
After traditional DE evolution, the individual may violate the boundary constraints and . is replaced by new individual being adjusted as follows:
The new population is combined with the existing parent population to form a new set of bigger size than . A nondominated ranking of is performed, and the best individuals are selected. This approach allows a global nondomination checking between both the parent and the new generation rather than only in the new generation as is done in other approaches, whereas it requires additional computational cost in sorting the combined.
2.1.2. Nondominated Ranking Based on Euclidean Distance
The solutions within each nondominated frontier that reside in the less crowded region in the frontier are assigned a higher rank, as the NSGAII algorithm [8] developed by Deb et al. indicated. The crowding distance of the th solution in its frontier (marked with solid circles) is the average side length of the cuboids (shown with a dashed box in Figure 1(a)). The crowdingdistance computation requires sorting the population according to each objective function value in ascending order of magnitude. As shown in Figure 1, and are two solutions near in the same rank, and is the crowding distance of the th solution, traditionally calculated as follows: where , are the objective vectors. For each objective function, the boundary solutions (solutions with the smallest and the largest function values) are assigned an infinite distance value.
(a)
(b)
(c)
(d)
A crowdingdistance metric is used to estimate the density of solutions surrounding a particular solution in the population and is obtained from the average distance of the two solutions on either side of the solution along each of the objectives. As shown in Figure 1, , , are the individuals of the generation on the same frontier, and we easily know that the density in Figure 1(a) is better than that in Figure 1(b). If we use (4) to calculate the crowding distance of , we only know that in Figure 1(a) it is better than in Figure 1(b); the crowding distance of in Figures 1(a) and 1(c) is equal, which is against the knowledge.
To distinguish the mentioned situations, we propose an improved crowdingdistance metric based on Euclidean distance. is the center point of the line , is the objective vector, and the crowding distance is defined as follows:
The crowdedcomparison operator guides the selection process at the various stages of the algorithm toward a uniformly distributed Paretooptimal frontier. To carry out the comparison, we assume that every individual in the population has two attributes: nondomination rank and crowding distance . Then, a partial order is defined as . If , that is, between two solutions with different nondomination ranks, we prefer the solution with the lower (better) rank, namely, . Otherwise, if both solutions belong to the same frontier, that is, , then we prefer the solution that is located in a lesser crowded region, that is, .
2.1.3. Simplex Method for Local Search
The simplex method for local search is mixed in the evolution process to enhance the searching strategy in the optimization process. The goal of integrating NM simplex method [32] and DE is to enrich the population diversity and avoid being trapped in local minimum. We apply simplex method operator to the present population when the number of iterations is greater than a particular value (like ). The individuals that achieved the single extreme value in each objective function are marked as the initial vertex points of simplex method. The present population is updated according to simplex method until the terminal conditions are satisfied.
The computation steps of the algorithm are included in Section 3.2.
2.2. Evaluation Criteria
Unlike the singleobjective optimization, it is more complicated for solution quality evaluation in the case of multiobjective optimization. Many of the suggested methods can be summarized in two types. One is to evaluate the convergence degree by computing the proximity between the solution frontier and the actual Pareto frontier. The other is to evaluate the distribution degree of the solutions in objective space by computing the distances among the individuals. Here, we choose both methods to evaluate the performance of the SMDEMO algorithm.
(1) Convergence Evaluation. Deb et al. [8] proposed this method in 2002. It is described as follows: where is the extent of convergence to a known of Paretooptimal set, is the obtained nondomination Pareto frontier, is the real nondomination Pareto frontier, is the Euclidean distance of with , and is the number of obtained solutions.
(2) Distribution Degree Evaluation. The nonuniformity in the distribution is measured by as follows: where is the Euclidean distance among consecutive solutions in the obtained nondominated set of solutions and parameter is the average distance.
2.3. Experimental Studies
Four wellknown benchmark test functions [33] are used here to compare the performance of SMDEMO with NSGAII, DEMO/Parent. These four problems are called ZDT2, ZDT3, ZDT4, and ZDT6; each has two objective functions. We describe them in Table 1.

The simulation is carried out under the environment of Intel Pentium 4, CPU 3.06 GHz, 512 MB memory, Windows XP Professional, Matlab7.1. Initialization parameters are set as follows: population size , scaling factor , cross rate , maximum evolution generation , and number of SM evolution iterations .
All of the three algorithms are real coded, with equal population size and equal maximum evolution generation. Each algorithm independently runs 20 times for each test function. Because we cannot get the real Paretooptimal set, we will take 60 Paretooptimal solutions obtained by the three algorithms as a true Paretooptimal solution set.
We evaluated the algorithms based on the two performance indexes and . Table 2 shows the mean and variance of and using three algorithms: SMDEMO, NSGAII, and DEMO/Parent. We can learn from Table 2, for the ZDT2 function, that all of the three algorithms have a good performance, while the SMDEMO is slightly better than the other two algorithms. In terms of convergences, for ZDT3, ZDT4 and ZDT6, which are more complex than ZDT2, SMDEMO is significantly better than DEMO/Parent and NSGAII.

Figure 2 shows a random running of SMDEMO algorithm. It is clear that SMDEMO algorithm can produce a good approximation and a uniform distribution.
(a) ZDT2
(b) ZDT3
(c) ZDT4
(d) ZDT6
3. Proposed Hybrid Algorithm for CMOP
The space of constrained multiobjective optimization problem can be divided into the feasible solution space and the infeasible solution space, as shown in Figure 3, where is the search space, is the feasible solution space, and is the infeasible solution space. is the feasible solution, and is the infeasible solution. Assume that is the global optimal solution and is the closest one to . If the infeasible population is not excluded by the evolution algorithm, it is permitted to explore boundary regions from new directions, where the optimum is likely to be found.
3.1. General Idea of the Proposed Algorithm
Researchers have gradually realized the merit of infeasible solutions in searching for the global optimum in the feasible region. Some infeasible solutions with better performance are allowed to be saved. Farmani et al. [34] formulated a method to ensure that infeasible solutions with a slight violation become feasible in evolution. Based on the constraints processing approach of multiobjective optimization problems, the proposed hybrid DE algorithm avoids constructing penalty function and deleting meaningful infeasible solutions directly.
Here, the proposed algorithm will produce multiple groups of functional partitions, which include an evolutionary population of size , an intermediate population to save feasible individuals, an intermediate population to save infeasible individuals, a population to save the optimal feasible solution found in the search process and a population to save the optimal infeasible solution. The relationship of multipopulation is shown in Figure 4.
With the description of (1), equality constraints are always transformed into inequality constraints as , where and is a positive tolerance value. To evaluate the infeasible solution, the degree of constraint violation of individual on the constraint is calculated as follows:
The final constraint violation of each individual in the population can be obtained by calculating the mean of the normalized constraint violations.
In order to take advantage of the infeasible solutions with better performance, we proposed the following adaptive differential mutation operator to generate individual variation learning from the mutation operators in DE/randtobest/1/bin, according to rules defined by Price et al. [11]. Considering each individual vector , a mutant individual is defined by where and represent different integers and also different from , randomly generated over ; is the scaling factor; is randomly generated from , is randomly generated from ; and is the mutation factor as follows: where is the initial value of the variability factor, is a small constant, to ensure that the fractional is meaningful, and is defined as follows:
3.2. Framework of the Proposed Algorithm
The proposed algorithm is described as follows.
Step 1 (initialization). Generate the population , , and of size , , and . Set the value of (crossover probability), (the number of function evaluations), (the iterative number of evolution by simplex method), (the current generation number), and positive control parameter for scaling the difference vectors , . Randomly generate the parent population from the decision space. Set the , and , and let the intermediate populations and be empty.
Step 2 (DE reproduction). By (3) and (9) for mutation, crossover, and selection, an offspring is created. Judge the constraints of all individuals in . In accordance with (8), we first calculate constraint violation degree of all of the individuals. If , the solution is feasible and preserved to the intermediate set ; if , the solution is infeasible and preserved to the intermediate set .
Step 3 (simplex method). Apply NM simplex method operator to the present population if . Update the present population when the number of iterations exceeds maximum iterations.
Step 4 ( construction). Rank chromosomes in based on (5), and generate the elitist population (the size is ) from the ranked population .
Step 5 ( construction). Add the chromosomes in with slight constraint violations to the .
Step 6 (mixing the population). Combine with the existing parent population to form a new set . Remove the duplicate individuals in and the existing parent population.
Step 7 (evolution). Randomly choose chromosomes from , , and . Use the adaptive differential mutation and uniform discrete crossover to obtain the offspring population .
Step 8 (termination). If the stopping criterion is met, stop and output the best solution; else, go to Step 2.
3.3. Experimental Study
In this section, we choose three problems CTP, TNK, and BNH, as shown in Table 3, to test the proposed method, and compare the method with the current CNSGAII [35].

For CTP problem, there are the six parameters , , , , , and that must be chosen in a way so that a portion of the unconstrained Paretooptimal region is infeasible. Each constraint is an implicit nonlinear function of decision variables. Thus, it may be difficult to find a number of solutions on a nonlinear constraint boundary. We take two sets of values for six parameters in CTP problem, which are determined as CTP1: , , , , and ; CTP2: , , , , , and . The Pareto frontiers, the feasible solution spaces, and the infeasible solution spaces are shown in Figure 5.
(a) CTP1
(b) CTP2
The parameters are initialized as follows. The size of population is , size of is , size of is , scaling factors and are randomly generated within , cross rate is , maximum evolution generation is , and number of SM evolution iterations is . All of the proposed algorithms and CNSGAII are real coded with equal population size and equal maximum evolution generation. Each algorithm runs 20 times independently for each test function.
Figure 6 shows the result of a random running of the proposed algorithm and NSGAII, the smooth curve “—” represents the Pareto frontier, and “◆” stands for the solution achieved by the proposed algorithm or NSGAII.
(a) CTP1 calculated by the proposed algorithm
(b) CTP1 calculated by CNSGAII
(c) CTP2 calculated by the proposed algorithm
(d) CTP2 calculated by CNSGAII
(e) BNH calculated by the proposed algorithm
(f) BNH calculated by CNSGAII
(g) TNK calculated by the proposed algorithm
(h) TNK calculated by CNSGAII
It is obvious that the proposed algorithm returns a better approximation to the true Paretooptimal frontie and a distribution of higher uniformity. We also evaluated the algorithms based on the two criterions and , as shown in Table 4. It can be observed from the data in Table 4 that the proposed algorithm performs significantly better than the classical CNSGAII algorithm in convergence and distribution uniformity. The simulation results show that this algorithm can accurately converge to global Pareto solutions and can maintain diversity of population.

4. Optimization of Grinding and Classification Process
4.1. Bauxite Grinding and Classification Process
The grinding and classification process is the key preparation for the bauxite mineral processing. Here, we consider a bauxite grinding process in a certain mineral company with single grinding and twostage classification, as shown in Figure 7.
The process consists in a grinding ball mill and two spiral classifiers. First classifier recycle will be put back to the ball mill for regrinding, and the firststage overflow will be put into second spiral classifier after being mixed with water; the second classifier recycle will be prepared for Bayer production as the rough concentrate, and the secondstage overflow will be sent to the next flotation process. The production objectives are composed of the production yields, technically represented by the solid flow of feed ore since the process is nonstorable, and the mineral processing quality, represented by percentage of the smallsize fractions of mineral particles in the secondstage overflows.
4.2. Predictive Model of the Grinding and Classification Process
Here, we establish the mathematical predictive model of each unit process in the bauxite grinding and classification process. The notations of the indexes, decision variables, and parameters are listed in Table 5. These notations will be used for the model of the grinding and classification process.

4.2.1. Ball Mill Circuit Model
Here, is the particle percentage of th size fraction in the ball mill overflow, is the particle percentage of th size fraction in feed ore, rate of the first classifier recycle is known, and is the efficiency of the first spiral classifier. According to a technical report of field investigation and study, we have that where is the mean residence time of minerals, is the internal concentration in ball mill, and where () is the solid flow of feed ore, is the water addition of the first classifier recycle, and is the classifier water addition. is the breakage distribution function; is the breakage rate function, and it satisfied the following equation: where is the particle with the th size, it is a unit, when per millimeter is a unit, , , and , , , and are four key parameters to control the breakage rate function.
In a concrete grinding and classification process, the ball mill size is fixed, and the speed of ball mill is constant. Through data acquisition and testing of grinding and classification steadystate loop, the regression model between , , , and condition variables, size fraction distribution can be established. The input variables are ball filling rate , solid flow of feed ore , water addition of the first classifier recycle , and parameters of feed ore size fraction distribution , . The regression model is
The value of can be obtained by the experimental data regression, , can be obtained from feed ore size fraction distribution, and is the cumulative particle percentage less than the th size fraction in feed ore, and it is represented as follows:
4.2.2. Spiral Classifier Model
is the particle percentage of the th size fraction in the first classifier overflow, and is the particle percentage of the th size fraction in ball mill overflow. The spiral classifier model is as follows: where is the efficiency of the first spiral classifier and the mechanism formula of is shown as follows: where is the particle with the th size, and represent maximum and minimum particle sizes, is the particle size fraction after correction separation, is separation accuracy, and is intermix index.
Through data acquisition and testing of grinding and classification steadystate loop, the regression model between classification parameters and condition variables, size fraction distribution can be established. The input variables include the solid flow of feed ore , the classifier water addition , and the parameters of ball mill overflow size fraction distribution , . The regression model is shown as follows: where the value of can be obtained by data regression.
The firststage overflow calculation formula is as follows:
Similarly, we can get the second spiral classifier model as follows: where is the particle percentage of the th size fraction in the second classifier overflow, is the particle percentage of the th size fraction in the first classifier overflow, and is the efficiency of the second spiral classifier. The spiral classifier model is as follows: where, , , , and are key parameters to the efficiency of the second spiral classifier. Through data acquisition and testing of grinding and classification steadystate loop, the regression model between classification parameters and condition variables, size fraction distribution can be established. The input variables include solid flow of feed ore and parameters of the firststage classifier overflow size fraction distribution , , which are solved by similar equation to (20). The regression model is shown as follows: where the value of can be obtained by experimental data regression.
4.3. Optimization Model of Grinding and Classification Process
Two objective functions in the process are identified: one is to maximize output , and the other is to maximize the smallsize fractions (less than 0.075 mm fractions) in the secondstage overflow . It is also necessary to ensure that the grinding product meets all of other technical requirements and the least disturbance in the following flotation circuit. As the constraints, the feed load of the grinding circuit , the steel ball filling rate , the first and the second overflows and , and the particle percentage of fine size fraction in the first and the second classifier overflows and should be within the user specified bounds.
The operation variables are the solid flow of feed ore , water addition of the first classifier recycle , ball filling rate , and water addition of the second classifier . Based on all of the above, grinding and classification process multiobjective optimization model is as follows:
In (24), implicates the model of grinding and classification represented by (12)–(23). The proposed algorithm is applied to solve the problem, and the optimization results are shown in Table 6.

With the practical process data from a grinding circuit of a mineral plant, the simulation of this hybrid intelligent method adopted the same parameters on the variation in fresh slurry feed velocity, density, particle size distribution, and cyclone feed operating configurations.
The comparison of production data and optimization results in Table 6 is shown in Figure 8, where “◆” represents the proposed algorithm optimization results and “○” represents the original data collected from the field without optimization of raw data. According to the objectives, the data point closer to the upper right edge is more beneficial. Obviously, the proposed optimization result is far better than the original data, indicating the effectiveness of the optimization approach.
4.4. TOPSIS Method for Solution Selection
The resolution of a multiobjective optimization problem does not end when the Paretooptimal set is found. In practical operational problems, a single solution must be selected. TOPSIS [36] is a useful technique in dealing with multiattribute or multicriteria decisionmaking (MADM/MCDM) problems in the real world. The standard TOPSIS method attempts to choose alternatives that simultaneously have the shortest distance from the positiveideal solution and the farthest distance from the negativeideal solution. According to the TOPSIS method, the relative closeness coefficient is calculated, and the best solution in Table 6 is the solution number 10 as and