Fuzzy robust optimization for handling feed stream and model parameter uncertainties during comminution process

Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

Contents lists available at ScienceDirect

Journal of the Taiwan Institute of Chemical Engineers journal homepage: www.elsevier.com/locate/jtice

Fuzzy robust optimization for handling feed stream and model parameter uncertainties during comminution process Nagajyothi Virivinti, Kishalay Mitra∗ Department of Chemical Engineering, Indian Institute of Technology Hyderabad, Kandi, Sangareddy, Medak 502 285, Telangana, India

a r t i c l e

i n f o

Article history: Received 13 May 2016 Revised 31 August 2016 Accepted 24 October 2016 Available online 21 November 2016 Keywords: Fuzzy programming Expected value Robust programming Pareto Uncertainty Industrial grinding process

a b s t r a c t Uncertainty analysis of an industrial grinding optimization process involving various sources of uncertainties and multiple numbers of objective functions has been studied in this work. Two primary sources of uncertainties, e.g. (1) operational parameters such as feed size distributions that are subjected to uncertainty due to varied range of feed sources that the industrial grinding process handles, and (2) model parameters that are obtained by the regression analysis of experimental data and subjected to regression and experimental errors, have been considered. Moreover, obtaining statistical distributions for these uncertain parameters are practically challenging. Hence, these parameters are considered as fuzzy numbers and the embedded uncertain multi-objective optimization problem has been analyzed using credibility based fuzzy robust optimization (FRO) technique. This technique helps in converting the uncertain fuzzy optimization formulation into a deterministic equivalent form, which can be further utilized to obtain the Pareto solutions by well-developed evolutionary algorithms. Initially, PO solutions are obtained by considering the uncertainty in different sets of uncertain parameters separately. This is followed by the study of amalgamated effect of uncertain parameters on PO solutions. Along with this, different fuzzy measures e.g. credibility, possibility, necessity etc. are utilized to observe their effects on final solutions. PO solutions obtained from possibility and necessity based FRO show the optimistic and pessimistic attitudes of risk, respectively. This provides a key to a decision maker to select any point based on the existing risk appetite. As compared to the deterministic results, the robust Pareto solutions show potential improvements of 26% and 16% in throughput and mid-size fractions, respectively. This generic approach not only provides the solution for robust range of operation based on the risk appetite of the enterprise, but also helps a decision maker to decide which parameters, amidst the set of uncertain parameters, are more sensitive to the results utilizing various fuzzy measures such as credibility, possibility and necessity. Along with this, the PO solutions obtained with robust optimization are compared with the solutions obtained by Expected Value Model (EVM). © 2016 Taiwan Institute of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

1. Introduction One of the reasons behind the popularity of building a mathematical model of a chemical process is its ability to simulate various scenarios and thereby saving huge cost and time behind conducting each one of them experimentally. Mathematical models thus built often have a set of parameters embedded in them that are generally determined by tuning the model with experimental results to mimic the model behavior close to the process conditions. As these model parameters are regressed from the experimental data, they are subjected to uncertainty related to the experimental and regression errors; omission of such errors could lead to unrealistic results [1]. For example, while ∗

Corresponding author. E-mail addresses: [email protected], [email protected] (K. Mitra).

conducting optimization based studies using one of these models, the formulations are most frequently established in a deterministic fashion where it is ensured that all parameters appearing in the optimization formulation, other than the decision variables, are certainly known and their values cannot be changed during the entire course of optimization [2]. These parameters include the coefficients in the objective function as well as constraints e.g. coefficients A, b and c in an optimization formulation of type Min cx, subjected to Ax ≤b. However, if some of these parameters are exposed to real life uncertainty, the results obtained assuming them as constants may lead to infeasibility when they are changed even slightly from their existing values [3]. Conventionally, these uncertainties can be handled by overdesigning the equipment or overestimating the operational parameters and thereby absorbing the uncertainty embedded in them. Another popular approach is to replace the uncertain parameters with their nominal values and

http://dx.doi.org/10.1016/j.jtice.2016.10.046 1876-1070/© 2016 Taiwan Institute of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

412

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

solve the deterministic formulation. These ideas are not pragmatic because the former leads to a very conservative solution where full potential of the plant/equipment is never utilized, whereas in case of the latter, the benefit of the changing operating conditions can never be exploited leading to either missing opportunities or unnecessary inventory piling due to under and over production situations. The paradigm of optimization under uncertainty (OUU) has, therefore, become popular as it allows a decision maker to create reliable solutions that remain feasible in the presence of parametric uncertainty. Broadly, two types of OUU methodologies exist in the literature [4]. The first one, the probability based approaches, assumes that the uncertain parameters have to follow some statistical distribution. Stochastic programming (SP), chance constrained programming (CCP) and, expected value model (EVM) [4] belong to this category. Amidst them, the most popular is the scenario based two-stage SP (TSSP) [5] technique, where the decision variables can be deployed in two stages, one before the realization of uncertainty (“here and now”) and the other after the realization (“wait and see”). As a result, the objective function, represented as the sum of the first stage components and expectation of second stage cost components, involving both first and second stage decision variables. One of the major drawbacks of this process is that the problem size increases exponentially as the numbers of uncertain parameters and scenarios increase, thus leading to an uncontrollable situation under the finite computational resources. Apart from this, it may not be easy always to decompose the problem into two sets of decision variables. Unlike SP, CCP, first introduced by Charnes and Cooper [6], provides a competitive way to solve OUU problems in terms of problem size and execution time management. This technique introduces the probability of constraint satisfaction, which can be connected to the reliability of the solution obtained, assuming that a constraint may not be satisfied for all uncertainty realizations. To make the problem more tractable, the probabilistic formulation is converted into an equivalent deterministic optimization problem (EDOP) formulation, which can be solved by conventional optimization techniques. The size of the EDOP is generally manageable even in the presence of a large number of uncertain parameters unless the probability of constraint satisfaction is set to be very high (i.e. very close to unity). Conversely, the EVM uses expected values of the objective function and the constraints. In a nutshell, the assumption that the uncertain parameters have to follow some well-behaved, known probability distributions and that can be estimated quite accurately in the above approaches usually requires an astronomical, completely unrealistic number of observations and might not be a practical task under real life situations [3]. On the contrary, fuzzy based approaches [7] help to avoid most of the afore-mentioned lacuna. Expressing the uncertain parameters as fuzzy variables, the extent of constraint violation can be represented by defining a membership function between 0 and 1, where 0 signifies maximum extent of constraint violation and 1 signifies no violation. Intermediate violations can be interpolated using different kinds of linear or non-linear fuzzification approaches without much restrictions. Fuzzy chance constrained programming (FCCP) [8] uses possibility of constraint satisfaction instead of probability measure in CCP. Similarly, in fuzzy expected value model (FEVM), expected values are calculated using possibility measure. Along with possibility, other fuzzy measures, e.g. necessity, credibility etc., representing various attitudes of a decision maker, can also be used to analyze the uncertain scenarios. Dubois and Prade [9] defined fuzzy variable as a mapping element from the space of possibility or necessity to the space of real numbers. To define an intermediate state of decision making and overcome the drawback of not owning the property of self-duality, Liu [10] introduced the credibility measure, defined as a weighted

average of the possibility and necessity measures. Applications of credibility based FCCP and FEVM are available in the literature [11– 13]. However, going by the very definition of a robust solution, which is a fixed decision variable vector that should remain feasible irrespective of the realization of the uncertainty in the parameters [3], the solutions obtained by above approaches may not be robust. So, robust optimization is an alternative way to handle OUU problems and might be extremely important under situations where the cases of constraint violation are strictly restricted. Credibility based fuzzy robust optimization is one technique which can be used to transform the uncertain optimization problem into its equivalent deterministic optimization problem. In mineral processing industries, grinding is one of the most popularly used size reduction processes. However, it is highly energy-intensive process, where any possibility for reduction of energy consumption can provide both economic edge and carbon footprint advantages. Thus, modeling and optimization of grinding have been of utmost importance for finding out the best operating conditions. Significant work have been done by several researchers [14–15] on modeling of grinding circuit. As per these works, the grinding circuit comprises operations such as rod mill, ball mill etc., as well as separation units such as hydro-cyclones, where modeling of the grinding operation is the most crucial. Bond [16] attempted to summarize the calculation methods in the grinding process using power based models. Despite being advantageous in terms of prediction accuracy as well as simplicity in approach, these models suffer from extrapolation across various dimensions of design changes. Austin et al. [17] and Morrell et al. [18] presented advanced modeling concept by introducing population balance methods encompassing selection functions, breakage functions and discharge functions. These parameters allow the designer to monitor and track particle breakage and transition within several size classes with time inside the grinding mill. While exploiting the improvements in the computational power, Chakraborti et al. [19,20] investigated the modeling of particle–particle interactions and particle flow patterns inside the grinding unit using computational fluid dynamics. Mishra [21], Mishra and Rajamani [22] and Datta and Rajamani [23] further analyzed particle–particle interaction behavior by discrete element modeling method. For the above computational techniques, the decision maker has to specify the grid points, boundary conditions and a set of starting points. Based on the starting solutions, states of the particle are calculated at all grid points honoring the boundary conditions and accurate results are obtained at the cost of computational time and rigor. Similarly, modeling of hydro-cyclones was presented by several researchers [24,25], which have been successful in predicting the split of the material coming inside the hydro-cyclone unit as a feed into the overflow and underflow streams. Many of these unit operations are considered simultaneously in an integrated grinding circuit, where a population balance approach can be used for simulating the circuit performance [26–28] and time dependent behavior of different size classes can be predicted within reasonable computational time. Wei and Craig [15] conducted a survey on the development of control philosophy of grinding where various sets of manipulated variables (mill solids flow rate, sump and mill water flow rate, sump slurry flow rates), control variables (cyclone overflow product density, feed ratio, slurry level in the sump, cyclone feed pressure, product particle size) and control objectives (throughput maximization) were surfaced. PID controllers are predominately used in the grinding mills. Though the studies on optimization are ubiquitous among various fields [29–32], most of them are deterministic in nature, where the effect of uncertain parameters on the optimization results is ignored for the sake of simplicity. In this work, credibility based fuzzy robust optimization (FRO) technique has been applied to an industrial grinding case study

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

considering uncertainties in two kinds of most generic parameter sets such as operational and model parameters. One of the uncertain operational parameters is represented by the feed size distributions, whereas uncertainties present in model parameters are considered for parameters such as grindability indices, grindability exponents and sharpness indices. These parameters are crisply known when the problem is solved in deterministic fashion. However, the model parameters are subjected to uncertainty related to experimental as well as regression errors, whereas the operational parameters are subjected to uncertainty due to varied range of feed sources that the industrial process is exposed to. At the same time, distribution information of these uncertain parameters are difficult to obtain. Several researchers [12,13,33–35] utilized FCCP and FEVM on real life applications and found that solutions of an uncertain optimization problem using these approaches address different aspects of uncertainty. However, none of these techniques may give the robust solution because they are formulated to tackle the optimization problem differently. This is primarily because some of the above approaches e.g. CCP, FCCP etc. work on the principle of varied degree of constraint violations, whereas EVM does not consider the system variance arising due to uncertainty in the parameters. Hence, robust optimization is one of the important uncertainty handling techniques used to obtain robust solution with respect to uncertainties involved in any optimization problem. Credibility based robust optimization is used here to convert the uncertain optimization problem into an EDOP [36]. The case study considered here is a multi-objective optimization problem, which has two conflicting objective functions. The EDOP contains four objective functions, which are expected values and standard deviation values of the original two objective functions. The EDOP is solved to analyze the Pareto trade-off among the objectives, where uncertain parameters influence both the constraints and the objective functions. Moreover, the uncertain parameters are related in a nonlinear fashion with one another. Along with uncertainty analysis, Pareto optimal (PO) fronts obtained by possibility and necessity measures provide a key to decision maker to select any solution based on the prevailing condition and risk appetite. Additionally, the aspect of how to conduct parametric sensitivity in a formal way using OUU approach has become an outcome of this study. Nondominated sorting genetic algorithm (NSGA II) [37] has been used in this work to solve the converted deterministic multi-objective optimization problem equivalent to the fuzzy uncertain formulation. To the best of the knowledge of the authors, fuzzy logic based evolutionary multi-objective robust optimization of grinding operation under nonlinear parametric uncertainty is very rare.

2. Formulation 2.1. Process and model description The industrial comminution process under consideration has the following four units: rod mill, ball mill, hydro-cyclones, and water sumps. The ore drawn out from the mine is crushed in a crusher and stored in a storage bin. Fresh feed from the bin is fed to the rod mill along with water. The slurry generated from the rod mill is mixed with the slurry from the ball mill in a primary sump. The primary sump outlet stream is sent to the primary cyclone. The overflow from the primary cyclone goes to the secondary sump and the underflow is taken as a feed to the ball mill. The slurry generated in the secondary sump is taken to another hydro-cyclone, called as secondary cyclone. The underflow of the secondary cyclone is recycled back to the ball mill for grinding and the final product is the secondary cyclone overflow which goes to a flotation circuit as feed. Water is added to both sumps to facil-

413

itate the flow of the slurry smoothly within the circuit. Complete circuit configuration can be found in Fig. 1. Modeling of individual unit operations of the grinding circuit is performed separately using an amalgamated approach of population balance [17,18] and empirical correlations. A simulation of an entire circuit is done by using a connectivity matrix which connects all the unit operations in terms of binary numbers. Here 0 denotes no connection and 1 denotes existence of a connection. In general material balance encompassing all unit operations, the population balance approach, is considered where several important size classes across the entire circuit are considered. However, for calculating breakage functions in the rod and the ball mill and sharpness indices in multiple numbers of parallel hydrocyclones, empirical correlations are used. Selection functions used in a number of different unit operations are supposed to follow nth order kinetics where n is to be determined from the industrial data. Overall construction of the dynamic model of the grinding circuit goes like this. Dynamic equations of individual size classes (five in number) along with the overall mass balance equation define the six differential equations for the rod mill. Similarly, there are six such equations each for ball mill, primary sump and secondary sump altogether leading to 24 differential equations. There are five corrected efficiency equations for calculating efficiency factors for each of the five size classes in each of the primary and the secondary hydro-cyclones separately giving rise 10 more algebraic equations. In addition to this, there are two more algebraic equations (pump factor states for the primary sump and the secondary sump), which makes it a set of total 36 differential algebraic equations (DAE) describing the behavior of the grinding circuit. This entire set of DAEs is solved using well tested public domain software, called DASSL [38]. Details on these model equations can be found elsewhere [8] and not attached here for the sake of brevity.

2.2. Deterministic multi-objective optimization problem Grinding unit feeds the final product slurry to the following flotation circuit where size plays a very important role for selective flotation of intended metals. To maintain the correct functioning of the flotation circuit, supplying slurry material within the desired range of size (mid-size, MS, middle size class among the five size classes considered), representing quality, consistently from the grinding circuit is extremely important. Productivity, represented by the circuit throughput (TP), is another important issue that has to be considered at the same time. It has been already established [28] that there lies a trade-off between these two entities (MS and TP), when the aim is to maximize both of them. Simultaneous maximization of the above mentioned two objectives builds the basis for deterministic multi-objective optimization. Upper bound constraints are selected on the control variables such as percent solids (PS), size distribution (percent passing of coarse (CS) and fine (FS) size classes i.e. the largest and the smallest among the five size classes considered), and recirculation load (RCL) of the final output stream. Raw ore feed flowrate to rod mill (S), and flow rate of water to primary sump (W1 ) and secondary sump (W2 ), the manipulated variables, can be changed within their respective upper and lower bounds to achieve the objectives. Unnecessary circulation of slurry material within the circuit can be kept in check with the help of the upper bound constraint for recirculation load, which also means energy saving and reducing the carbon count. On the other hand, the upper bound constraints for percent solids and different size classes ensure the desired mineral liberation in the following flotation operation. With this background, the deterministic multi-objective optimization formulation can be represented as follows:

414

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

Fig. 1. An industrial grinding circuit block diagram.

Objective functions

Max

T P (,  )

Max

MS (,  )

S,W1 ,W2

S,W1 ,W2

These uncertain parameters are, therefore, assumed as fuzzy numbers and the OUU problem has been converted into EDOP using fuzzy based robust optimization. So, the EDOP with four objectives (maximization of the expected values and minimization of variance of the objectives) can be solved by using well developed evolutionary optimization techniques. It is worth mentioning here that uncertain parameters in this case are related nonlinearly with one another. Applications on handling uncertainty with non-linearly related uncertain parameters are relatively scarce in the literature. Approaches to handle such situations are completely different from the case where the uncertain parameters are related in linear fashion.

Subject to

C S (,  ) ≤ C SU F S (,  ) ≤ F SU P S (,  ) ≤ P SU

2.3. Robust optimization with fuzzy uncertain parameters

RC L (,  ) ≤ RCLU

Let us consider a general multi-objective optimization problem, which has m objectives and n constraints with some uncertain parameters in both objective functions and the constraints

Decision variables bounds

Maximize f1 (X, ξ ), f2 (X, ξ )... fm (X, ξ )

L

S ≤ S ≤ SU

Subjected to g j (X, ξ ) ≤ 0,

WL1 ≤ W1 ≤ WU 1 WL2 ≤ W2 ≤ WU 2

(1)

where  is a vector comprising of six uncertain parameters (φ 1 , φ 2 , φ 3 , φ 4 , φ 5 , φ 6 ) related to model uncertainty. They are grindability indices for the rod mill (RMGI) and the ball mill (BMGI), the grindability exponents for the rod mill (RMExpo) and the ball mill (BMExpo), the sharpness indices for the primary (PCSI) and secondary cyclones (SCSI), respectively.  is a vector which contains four uncertain parameters (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) related to operational uncertainty represented by the four size fractions fed to the rod mill. Out of the five size fractions considered, the fifth one can be calculated by subtracting the sum of these four size classes from unity. The deterministic optimization formulation can be solved assuming these parameters are crisply known; however, they are considered as uncertain parameters for an OUU formulation. Parameters represented by vector  are obtained from the regression of experimental data and thus are exposed to uncertainty due to experimental and regression errors, whereas parameters in vector  are related to inherent uncertainty of the operation. As the parameters are uncertain, there is certain amount of variance involved in them. However, the collection of data for quantifying the uncertainty involved in these parameters is difficult to achieve.

XLi ≤ Xi ≤ XUi ,

j = 1, ..., n

(2)

i = 1, ..., k

where ξ represents a vector of fuzzy parameters involved in the constraints and objective functions and X is the set of decision variables. The above optimization problem contains uncertain parameters in objective function and constraints. Hence, it is necessary to convert the OUU problem into its deterministic equivalent optimization problem. In this work, the OUU problem is converted to EDOP by using expected and standard deviation values of objective functions and constraints which is as shown below.

Maximize E { f1 (X, ξ )}, E { f2 (X, ξ )}...E { fm (X, ξ )} Minimize SD{ f1 (X, ξ )}, SD{ f2 (X, ξ )}...SD{ fm (X, ξ )} Subjected to [E {g j (X, ξ )} + SD{g j (X, ξ )}] ≤ 0, XLi ≤ Xi ≤ XUi ,

j = 1, ..., n

(3)

i = 1, ..., k

The expected value operator is used to calculate expected values and standard deviation values of the objective functions and constraints by using membership function values. Let X be a fuzzy number belonging to fuzzy set A and it has some value of membership function μA (X)→[0,1], which represents the degree of belongingness of a particular element into the fuzzy set A. Under this

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

condition, the fuzzy set A can be represented as [39]

A = {( X |μA (X ) ) : X ∈ R : μA (X ) ∈ [0, 1]}

(4)

Possibility measure shows the measure of the best case of the fuzzy event, which means the decision maker is optimistic and pos(A) gives the maximal chance that A holds. Similarly, necessity measure shows the measure of the worst case of the fuzzy event, which means the decision maker is pessimistic and nec(A) gives the minimal chance that A holds. Similarly, let f(X, ξ L ) be a function which contains a fuzzy variable ξ and r be a real number; possibility(f(X, ξ L ) ≥ r) is defined as degrees to what extent f(X, ξ L ) is not smaller than r and is calculated by

pos{ f (X,ξ ) ≥ r} = max

1≤L≤N

{ μL |{ f (X, ξL ) ≥ r}}

techniques are used to calculate the expected values of uncertain functions (objective functions and constraints) [40]. By using expected values of objective functions and constraints, standard deviations of objective functions and constraints can be calculated by using Eq. (11) shown above. 2.3.1. Fuzzy simulation algorithm for expected value model based on credibility

E = E { f (X, ξ )} or E {g(X, ξ )} Step 1: Assume e = 0; Step 2: Randomly generate N points of μL (L = 1, 2, …., N) from ε-set of ξ , where ε is a sufficiently small number and μ is the membership function; Step 3: Set the two numbers a = min (f(X, ξ 1 ), f(X, ξ 2 ), ......., f(X, ξ N )) and b = max (f(X, ξ 1 ), f(X, ξ 2 ), ......., f(X, ξ N )); Step 4: Randomly generate r between [a, b]; Step 5: If r ≥ 0, then e = e + cr(f(X,ξ ) ≥ r)); If r < 0, then e = e − cr(f(X,ξ ) ≤ r)); Step 6: Repeat the above two steps for N  times  Step 7: E = max(a, 0 ) + min(b, 0 ) + e.(b − a )/N

(5)

Possibility measure can be calculated by the maximum membership function value among all the constraint satisfaction cases inside the uncertain space. Necessity measure is defined as the impossibility of the complement set Ac which is the dual of possibility measure. With this background, necessity(f(X, ξ L ) ≥ r) is defined as degrees to what extent f(X, ξ L ) is not smaller than r and is calculated by

nec{ f (X,ξ ) ≥ r} = min { 1 − μL |{ f (X, ξL ) < r}} 1≤L≤N

(6)

Necessity measure can be defined as the minimum among complement of membership function values whenever constraint is not satisfied However, possibility and necessity fail to satisfy the property of self-duality. To fulfill this, credibility measure is used to represent the uncertainty in a fuzzy function. Credibility measure is defined as weighted average of possibility and necessity measures

cr { f (X,ξ ) ≥ r} = w1 × pos{ f (X,ξ ) ≥ r} + w2 × nec{ f (X,ξ ) ≥ r} (7) w1 and w2 are the weights of possibility and necessity measures which are in between 0 and 1 and the sum of w1 and w2 should be equal to one. These three different measures are used to calculate expected values and standard deviation values of objective functions and the constraints. Expected value can be calculated by using expected value operator which is defined here based on possibility measure

E Pos (ξ ) =



+∞ 0

pos{ξ ≥ r} −



0 −∞

pos{ξ ≤ r}

(8)

Similarly, necessity measure is used to calculate expected value of uncertain objective function that can be shown as below

E Nec (ξ ) =



+∞ 0

nec{ξ ≥ r} −



0 −∞

nec{ξ ≤ r}

ECr (ξ ) =

+∞ 0

cr {ξ ≥ r} −



0 −∞

cr {ξ ≤ r}



E[( f (X, ξL ) − e ) ] 2

As mentioned earlier, we assume six model parameters and four feed size fractions as uncertain in this formulation. Based on the description of robust optimization in Section 2.3, the OUU with fuzzy parameters can be converted [41] into equivalent deterministic non-fuzzy multi-objective optimization problem, which has four objective functions and four constraints. The four objectives are maximization of expected values of objective functions and minimization of standard deviation of objective functions. Hence, the converted EDOP can be as written as:

Max E (T P (,  ) ) Min SD(T P (,  ) ) Max E (MS (,  ) ) Min SD(MS (,  ) ) subject to

E (P S(,  )) + κ SD(P S(,  )) ≤ P SU E (CS(,  )) + κ SD(CS(,  )) ≤ C SU E (F S(,  )) + κ SD(F S(,  )) ≤ F SU E (RCL(,  )) + κ SD(RCL(,  )) ≤ RC LU

(10)

Standard deviation of uncertain functions can be calculated by using the expected values of uncertain functions. Let e be the expected value of f(X,ξ L ), then SD of f(X,ξ L ) can be calculated by using e value calculated earlier as

SD{ f (X, ξL )} =

2.4. Grinding optimization with fuzzy parameters

(9)

The expected value of uncertain function can be calculated by using credibility measure as given by



415

(11)

where e = E{f(X,ξ L )} In this work, the uncertain parameters are related to one another in non-linear manner and hence fuzzy simulation based

Decision variable bounds

S ≤ S ≤ SU L

WL1 ≤ W1 ≤ WU 1 WL2

≤ W2 ≤

(12)

WU 2

where the value of κ is assumed to be unity in this work. The value of κ can be decided based on the intended relative weightage on expected value and variance components in the constraints. The above deterministic equivalent can be solved with the help of a multi-objective optimization algorithm as well as the fuzzy simulation algorithm mentioned in Section 2.3.

416

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

3. Results and discussions The case study described here is motivated by an industrial Pb–Zn ore beneficiation operation. The model is validated with the steady state data collected from industry under various operating conditions and the results can be seen from the work of Mitra and Gopinath [28]. The deterministic optimization problem has two conflicting objectives and four constraints on the control variables as given in Section 2.3. The objective functions and constraints are function of uncertain parameters. Often, the statistical distribution information about these uncertain parameters is not available. Hence, the techniques based on statistical distributions cannot be applied to analyze the uncertainty involved in the parameters. One way out is to solve the uncertain formulation using fuzzy approach, where uncertain parameters can be represented by the fuzzy numbers. Optimization under uncertainty problem can be converted into a non-fuzzy equivalent deterministic formulation and the new deterministic problem can be solved using welldeveloped multi-objective optimization techniques. To observe the impact on the PO front, uncertainty has been considered 25% on the both sides, upper or lower, of the deterministic value of the uncertain parameters. The value of the membership function at the deterministic value is assumed to be 1, which decreases to zero at maximum allowable deviation values (e.g. 25% in this case) on either side. Fuzzy membership functions represented in this way can help a decision maker to understand whether any values of the uncertain parameters other than their deterministic values can improve or deteriorate the final PO solutions or not. The nature of uncertainty could have been represented by other fuzzy membership functions, e.g. trapezoidal, bell shaped etc. and the choice for linear membership function has been taken up as an example only. Moreover, the span of uncertainty level is kept same 25% on both sides for ease in analysis and can be changed as per decision maker’s choice. The deterministic equivalent of the FRO optimization problem has four objectives (mean and standard deviation of two objectives) and four constraints with three manipulated variables. The importance of considering the standard deviation component of the objectives and the constraints is that the decision maker aims to improve the performance of the operation with the minimum variation in the parameters. There lies a trade-off between the mean and the standard deviation values of the objective functions. 3.1. Effect of model related uncertain parameters First, the uncertainties in all six (φ 1 , φ 2 , φ 3 , φ 4 , φ 5 , φ 6 ) model parameters have been considered. PO fronts for lower and upper side uncertainties, called as “Type a” and “Type b”, respectively, can be observed in Fig. 2. Concentrating on the mean values of the objective functions, upper side uncertainty (“Type b”) is observed to provide with better PO front. The four-objective plot, where the fourth objective is presented by means of the shades in the color, has been presented by three subplots to observe the effect of uncertainty on the standard deviation. For upper side uncertainty (“Type b”), the favorable right and upward shift in the PO front has been observed (see Fig. 2b). The direction is favorable because the aim is to maximize both the objectives. For a fixed value of throughput, there is an increase in mid-size fraction passing values, i.e. a case of finer grinding has been recommended, compared to that of the deterministic case. Improvement in the mid-size fraction to the extent of ∼16.35% for a fixed throughput value can be observed while making a transition from the deterministic case to the case of upper side uncertainty (“Type b”). Similarly, for a constant mid-size fraction value, there is an increase in throughput when compared with the deterministic formulation i.e. to the extent of ∼20.43% improvement in throughput can be observed for a fixed mid-

size fraction value while making a transition from the deterministic case to upper side uncertainty (“Type b”). However, for the lower side uncertainty (“Type a”), it can be observed that the movement of the PO front is in the unfavorable left and downward direction. It indicates that the PO front obtained with deterministic values is better than that obtained with lower side of uncertain values. With this result, a decrease in both the objective function values has been observed, e.g. ∼26.9% decrease in the value of mid-size fraction for a fixed throughput value can be seen while making a transition from the deterministic case. Similarly, to an extent of ∼30.5% decrease in throughput value can be observed for a fixed mid-size fraction value while making a transition from the deterministic case. Lower values of standard deviation can be observed for lower side uncertain values from throughput point of view (Fig. 2c). However, expected values are low for the case of lower side uncertain values compared to expected values obtained with upper side uncertain values. Hence, there lies a trade-off between the objectives. Assuming the values of standard deviation are not significant, the effect of standard deviation can be hopefully neglected for throughput. On the contrary, lower standard deviation and higher expected values are observed for the case of upper side uncertain values from mid-size point of view (Fig. 2d). Hence, the values obtained with upper side uncertain values give better PO solutions compared to lower side uncertain (“Type a”) values from quality point of view (mid-size fraction). From throughput point of view, if we are allowed to ignore significantly small value of standard deviations, upper side uncertain values (‘Type b”) gives better PO solutions. A decision maker can select any point depending upon the selection of minimum risk with maximum performance. This means some values of uncertain parameters other than their deterministic values can also lead to better PO front—hence their uncertainty analysis is extremely important. This also shows if uncertainty is present in parameters, there is no merit in solving a deterministic optimization problem with the deterministic values, which is a common industry practice, avoiding a systematic optimization under uncertainty analysis like this. 3.2. Effect of operational related uncertain parameters Next, uncertainty in operational parameters ( ) has been considered to carry out the analysis with lower and upper side uncertainties (“Type a” and “Type b”) with 25% deviation from the deterministic values. As opposed to the previous case, “Type a” uncertainty has been observed to provide better PO solutions compared to “Type b” uncertainty (Fig. 3a and b). At the same time, standard deviation of throughput values are low for “Type a” uncertainty (Fig. 3c), though the corresponding magnitudes are insignificant, whereas “Type b” uncertainty has high standard deviation of mid-size fraction values (Fig. 3d). When the PO fronts between the expected values for two types of uncertainties are compared with that obtained with deterministic values, it has been observed that the expected values obtained with lower side uncertainty values, show the favorable right and upward movement in the PO front (Fig. 3b). For fixed value of throughput, an increase in value in the mid-size fraction passing has been observed and vice versa for “Type a” uncertainty and respective % improvement values can be observed in Table 1. However, the length of the PO front for “Type a” uncertain parameters has become short because of the constraint violation in the model (Percent solids and coarse size feed fraction constraints). Similarly, for “Type b” uncertain values, the unfavorable and the opposite trend emerges as there is deterioration in the Pareto front movement as given in Table 1. At the same time, the “Type a” uncertainty is found to be favorable from standard deviation point of view for both throughput and mid-size fraction (Fig. 3c and d). Hence in this case, it can be concluded that “Type a” uncertainty provides better PO solutions compared

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

417

Fig. 2. Credibility based four-objective Pareto optimal solutions (uncertainty in model parameters).

to “Type b” uncertain values considering both mean and standard deviation aspects. This basically means as compared to the current operation, if further finer feed can be provided as a feed to the plant, the performance can not only be improved further, but also be made more robust.

3.3. Amalgamated effect of model and operational related uncertain parameters Next, the combined effect of uncertain parameter sets (,  ), has been considered. Two types of uncertainty cases have been

418

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

Fig. 3. Credibility based four-objective Pareto optimal solutions (uncertainty in operational parameters).

considered again in the similar fashion. From Fig. 4a, it is very difficult to declare which type of uncertainty gives rise to better PO solutions due to the trade-off involved among them. However, from Fig. 4b, right upward movement in the Pareto front can be observed for lower side (“Type a”) uncertain parameters with respect to the deterministic case. As reported in Table 1, by considering the uncertainty in both the sets of parameters, for a fixed

value of throughput, there is increase in the value of mid-size fraction (15.9%). Similarly, for a fixed value of mid-size fraction, there is increase in throughput value (26.45%). However, it is favorable from the maximization point of view of both the objectives. Upper side (“Type b”) uncertain parameters show left and downward movement in the PO front. As reported earlier, the length of the PO front with “Type a” uncertain parameters has become shorter than

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

419

Table 1 Movement of Pareto front for different types of uncertain parameters. Type a Uncertainty in model related uncertain parameters

Uncertainty in operational related uncertain parameters

Amalgamated effect of uncertain parameters

Type b

Deterioration Mid-size fraction 26.9% Improvement Mid-size fraction 33.4% Improvement Mid-size fraction 15.9%

the previous case because of the constraint violation in the model (percent solids and coarse size feed fraction constraints). Further shortening of the PO front can be attributed to the opposing effect of “Type a” uncertainty for parameter sets  and  as reported before. From Fig. 4b and c, it can be observed that the optimal solutions obtained with “Type a” uncertain values gives better PO solutions because of the low standard deviation and high mean values. No significant difference could be observed between the standard deviation values for mid-size % passing for “Type a” and “Type b” uncertainty (Fig. 4d). This is closer to what is observed above while considering uncertainty only in the operational parameters. This clearly shows that the effect of uncertainty in operational parameters is more impactful in this case than that of the model parameters as the effect of the model parameters is suppressed when the uncertainty in  and  are considered together. On the basis of overall discussion, one may conclude that the solutions presented in Fig. 4 are more robust than that of Fig. 3 because the results presented in Fig. 4 account for the parameter uncertainty in the model, which itself is used to generate the results presented in Fig. 3. This shows that the robust optimization results predicted with a model that itself has uncertainty in its prediction should be analyzed properly and corrected against the uncertainty present in the parameters to generate more reliable results (clearly visible from the numbers presented in Table 1). This is precisely the reason for which some of the solutions taken from the “Type a” PO front of Fig. 3b might turn out to be infeasible if compared with the similar results presented in the Fig. 4b. 3.4. Sensitivity analysis However, it might be very difficult to identify a particular operational/model parameter that improves/deteriorates the circuit performance because the effects of individual parameters are all muddled together. Hence, it is necessary to consider the effect of individual parameter on the overall optimization problem. For that purpose, the uncertainty in parameters is considered separately. For example, the uncertainty only in the rod mill grindability index is considered first. All other nine uncertain parameters are assumed to be fixed at their deterministic values. Similarly, other model and operational parameters are considered one at a time. The effect of individual uncertain parameters on the PO front with “Type a” uncertain values is shown in Fig. 5. The improvement in the PO front is high in case of coarse size feed uncertainty, whereas the deterioration in the PO front is the maximum in case of grindability index of ball mill. Concentrating on the movement of the PO front from its deterministic position, it can be concluded that the coarse size feed input and ball mill grindability index are the most sensitive parameters among the overall set of ten uncertain parameters. From Fig. 5c, standard deviation values of mid-size fraction are high in case of uncertainty in grindability index of ball mill and coarse size feed input. Here also the same two parameters are surfaced as having higher standard deviation values, thereby showing the evidence of their being more sensitive on the overall

Throughput 30.5% Throughput 42.6% Throughput 26.45%

Improvement Mid-size fraction 16.35% Deterioration Mid-size fraction 38.5% Deterioration Mid-size fraction 22.88%

Throughput 20.43% Throughput 40.5% Throughput 24.89%

optimization results. The sensitivity of PO solutions with respect to other eight uncertain parameters is almost of similar and negligible order. This is another way of conducting sensitivity analysis for parameters in any model, which is definitely more systematic as compared to the general ad hoc approach of conducting parametric sensitivity by perturbing parameters randomly and observe their effects on the optimization result. 3.5. Effect of different measures (possibility, necessity and credibility measures) on the Pareto optimal solutions The uncertainty and sensitivity analysis carried out till now are based on the credibility measure, which owns the property of selfduality. Along with credibility measure, possibility and necessity measures are available in the literature. To observe the effect of possibility and necessity measures, credibility measure has been replaced by them. The PO fronts obtained with all these three measures with uncertainty in model parameters have been shown in Fig. 6. The PO front obtained with necessity measure gives better PO solutions compared to possibility and credibility measures based on expected values. The PO front obtained with necessity measure gives the aggressive nature of the Pareto front whereas the PO front obtained with possibility measure gives the conservative nature of the PO front (Fig. 6a). The credibility based PO front is found to lie between the possibility and the necessity PO fronts (Fig. 6a). However, the standard deviation values are found to be higher for the PO solutions obtained with necessity measure (Fig 6b and c). Hence, it can be concluded that the optimal solutions obtained with necessity measure are aggressive based on expected values whereas the optimal solutions obtained with possibility measure are better based on standard deviation values. Next, uncertainty in operational parameters has been considered to find out possibility, necessity and credibility based PO fronts. Unlike the previous case, the effect can’t be clearly seen, which might be due to the non-linearity present in the model (Fig. 7). 3.6. Comparison of PO solutions obtained with EVM and RO While converting the uncertain optimization problem into an equivalent deterministic optimization problem, FRO approach uses both the expected and standard deviation values of the constraints. In the above analysis, the relative importance of the expected and standard deviation values are kept same and for that the к value is assumed to be unity. To observe the effect of к on the optimization problem, optimal solutions have been found out for three different values of к (0, 0.5 and 1) and the results are shown in Fig. 8. With increase in к value, constraint violation (percent solids and coarse size feed fraction constraints) increases in the uncertain parameter space leading to shorter in length PO front, in terms of the length of it, emerge for high values of к. With increase in к value, constraints become tighter leading to shorter in length PO front. To carry out a comparative study of robust optimization and EVM, the PO solutions obtained by robust optimization are compared with

420

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

Fig. 4. Credibility based four-objective Pareto optimal solutions (uncertainty in model and operational parameters together).

the solutions obtained by EVM. The amalgamated effect of all ten parameters has been considered. EVM and RO techniques lead to better PO solutions compared to the deterministic PO solutions, represented by ‘certain’ in Fig. 9. However, the PO front obtained by RO got restricted because of constraint violation involved in the problem. Robust optimization technique is conservative compared to the EVM (neglects standard deviation involved in the problem).

However, robust optimization method is more advantageous as it considers both mean and standard deviation of objective functions and constraints. Based on the results presented so far, it is observed that multiobjective optimization problems have a set of Pareto solutions and this entire set might change its location in presence of the uncertainty in the problem. Though a single solution is to be selected

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

421

Fig. 5. Sensitivity analysis of uncertain parameters on the overall optimization result.

for implementation from this PO set of solutions, in presence of several objectives having trade-offs, the PO set has to be identified first. Since the location of the PO set itself is not fixed under uncertainty, it should be identified first in robust fashion. Through this study, it is shown that identification of the robust location of the PO set can be achieved through FRO based uncertainty handling technique. Once the robust location of PO set is identified, a single solution can be selected based on the selective preference of the objectives by a decision maker. This study could show that it has a potential of 26% and 16% improvement in throughput and mid-size fractions, respectively, which means the plant was being under-utilized without considering uncertainty and therefore was behaving conservatively. 3.7. Trends of manipulated variables Pareto solutions are arranged in the increasing order of throughput values and the trends of the manipulated variables are plotted with respect to the chromosome numbers. Considering uncertainty in model and operational parameters separately, these trends are plotted in Figs. 10 and 11, respectively. For the model parameter uncertainty, water flow rates to the primary sump are found to be confined within the middle to upper section of the entire range given (0.45–0.99 in the normalized scale), whereas the water flow rates to the secondary sump is more dispersed in the entire range given (0.27–0.93 in the normalized scale). Similarly, for the operational parameter uncertainty, water flow rates to the

primary sump are found to vary within the middle to upper section of the entire range given (0.58–0.98 in the normalized scale), whereas the water flow rates to the secondary sump are spread over the entire range given (0.07–0.97 in the normalized scale). The solid flowrates to the rod mill for operational parameters corresponding to the robust Pareto solutions are observed to vary in a smaller band (0.6–0.88 in the normalized scale) than that of the model parameters (0.33–0.88 in the normalized scale). So, based on the solid flowrates to the rod mill and the water flow rates to the primary sump, the case of operational parameter is relatively easier to control. The diversity in the nature of model parameters and their nonlinear relationship lead to more variation in the system, whereas the operational parameters are related to only size. Generally, these kinds of trends are obtained by any plant operator after running the plant for so many years. However, these trends have been found here as a result of the optimization exercise which can be used by operators for running the plant efficiently. 3.8. Size distribution of different steams The Pareto fronts for credibility based robust optimization approach for model and operational parameters are divided into two regions: lower throughput region and higher throughput region. Among different streams in the grinding circuit, some important streams such as ball mill discharge (BM Dsch), rod mill discharge (RM Dsch), primary cyclone overflow (PC OFlo), primary cyclone

422

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

Fig. 6. Pareto optimal solutions based on Credibility, Possibility and Necessity measures (uncertainty in model parameters).

Fig. 7. Pareto optimal solutions based on Credibility, Possibility and Necessity measures (uncertainty in operational parameters).

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

423

Fig. 8. Pareto optimal solutions showing relative importance of mean and standard deviation on the constraints under RO formulation.

Fig. 9. Comparative study of Robust Optimization and Expected Value Model. Fig. 10. Trends of manipulated variables for uncertainty in model parameters alone (Upper side uncertain values).

underflow (PC UFlo), secondary cyclone overflow (SC OFlo) and secondary cyclone underflow (SC UFlo) are chosen and the stream properties are determined primarily in terms of particle size distribution. For the terminal points of the Pareto fronts obtained with credibility measure for both the cases (uncertainty in operational parameters and uncertainty in model parameters), size distributions are shown in Figs. 12 and 13 for the above mentioned important streams. The lower throughput Pareto points (1 and 3 corresponding to Figs. 6 and 7) result in more finer grinding performance compared to the higher throughput Pareto points (2 and 4 corresponding to Figs. 6 and 7). As we move from higher to lower

throughput regions, more % passing values for different size classes for all different streams are reported. 4. Conclusions In this work, the multi-objective optimization of an industrial grinding problem has been analyzed in an uncertainty framework using fuzzy robust optimization approach, resembling one of the practical scenarios when the distribution information of the uncertain parameters cannot be made available. The uncertainties

424

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

Fig. 11. Trends of manipulated variables for uncertainty in operational parameters alone (Lower side uncertain values).

Fig. 13. Normalized cumulative size distributions in different streams of the grinding circuit corresponding to (a) minimum and (b) maximum throughput PO points (3 and 4) in Fig. 7 corresponding to the credibility measure Pareto curve (Uncertainty in operational parameters alone).

Fig. 12. Normalized cumulative size distributions in different streams of the grinding circuit corresponding to (a) minimum and (b) maximum throughput PO points (1 and 2) in Fig. 6 corresponding to the credibility measure Pareto curve (Uncertainty in model parameters alone).

in model predictions and operation have been assumed to appear through different parameters such as grindability indices of ball mill and rod mill, sharpness indices of primary cyclone and secondary cyclone, grindability exponents of ball mill and rod mill and four feed side fractions, which are considered as constants during deterministic optimization. The deterministic model is adopted from the work of Mitra and Gopinath (2004) and required modifications have been carried out on it to accommodate uncertainty in the model parameters. The following points can be concluded from the study:

1. It has been shown that deterministic multi-objective optimization problems have a set of Pareto solutions and this entire set might change its location in presence of the uncertainty in the problem. Through this study, it has been shown that identification of the robust location of the PO set can be achieved through FRO based uncertainty handling technique that considers the average change and variability in the uncertain parameters. 2. The robust location of PO set has shown scope of improvements of 26% and 16% in throughput and mid-size fractions, respectively as compared to the current deterministic PO solutions. This shows a tremendous scope of improvement of the robust plant operation under operational uncertainty. 3. The above results of robust optimization have been thoroughly analyzed by considering the uncertainty in model and operational parameters separately as well as together. This enabled attainment of the robust Pareto optimal solutions separating out the effect of uncertainty present in the parameters of the model itself. 4. To identify the most sensitive parameter for the Pareto movement, sensitivity analysis has been carried out in a novel fashion and found that the coarse size feed input is the most sensitive parameter. Moreover, the operational set of parameters is found to have more sensitivity than the model parameters. 5. Possibility and necessity measure based PO fronts provide a key to decision maker to operate the plant with the maximum per-

N. Virivinti, K. Mitra / Journal of the Taiwan Institute of Chemical Engineers 70 (2017) 411–425

formance or minimum risk as they represent two extremes of the uncertain space involved. 6. Along with this, the comparative study of robust optimization with expected value model shows that robust optimization results are more conservative compared to that of the expected value model, which neglects the component of variance in the uncertain parameters. 7. Trends among the decision variables have been established based on robust PO solutions which can help an operator to run the plant in optimal fashion. Generally, these trends are extremely useful for the operators, but hard to find without running the plant for a very long duration. References [1] Xu W, Diwekar UM. Improved genetic algorithms for deterministic optimization and optimization under uncertainty. Part II. Solvent selection under uncertainty. Ind Eng Chem Res 2005;44:7138–46. [2] Valipour M. Optimization of neural networks for precipitation analysis in a humid region to detect drought and wet year alarms. Meteorol Appl 2015;23:91–100. [3] Ben-Tal A, El Ghaoui L, Nemirovski A. Robust optimization. 41 William Street, Princeton, New Jersey: Princeton University Press; 2009. [4] Sahinidis NV. Optimization under uncertainty: state-of-the-art and opportunities. Comput Chem Eng 2004;28:971–83. [5] Liu ML, Sahinidis NV. Optimization in process planning under uncertainty. Ind Eng Chem Res 1996;35:4154–65. [6] Charnes A, Cooper WW. Chance-constrained programming. Manage Sci 1959;6:73–9. [7] Zimmermann HJ. Fuzzy set theory and its application. Boston: : Kluwer Academic Publishers; 1991. [8] Mitra K. Multi-objective optimization of an industrial grinding operation under uncertainty. Chem Eng Sci 2009;64:5043–56. [9] Dubois D, Prade H. Possibility theory: an approach to computerized processing of uncertainty. New York: Plenum Press; 1988. [10] Liu B. A survey of credibility theory. Fuzzy Optim Decis Making 2006;5:387–408. [11] Tian G, Chu J, Liu Y, Ke H, Zhao X, Xu G. Expected energy analysis for industrial process planning problem with fuzzy time parameters. Comput Chem Eng 2011;35:2905–12. [12] Virivinti N, Mitra K. Fuzzy expected value analysis for an industrial grinding process. Powder Technol 2014;268:9–18. [13] Virivinti N, Mitra K. A comparative study of fuzzy techniques to handle uncertainty: an industrial grinding process. Chem Eng Technol 2016;39:1031–9. [14] Powell MS, Morrison RD. The future of comminution modeling. Int J Miner Process 2007;84:228–39. [15] Wei D, Craig IK. Grinding mill circuits—a survey of control and economic concerns. Int J Miner Process 2009;90:56–66. [16] Bond FC. The third theory of comminution. Trans SME/AIME 1952;193:484–94. [17] Austin LG, Klimpel RR, Luckie PT. Process Engineering of size reduction: ball milling. New York: SME of AIME; 1984. [18] Morrell S, Sterns UJ, Weller KR. The application of population balance models to very fine grinding in tower mills. In: XVIII Int. Min. Proc. Cong. Sydney, 1; 1993. p. 61–6.

425

[19] Chakraborti N, Shekhar A, Singhal A, Chakraborty S, Chowdhury S, Sripriya R. Fluid flow in hydrocyclones optimized through multiobjective genetic algorithms. Inverse Probl Sci Eng 2008a;16:1023–46. [20] Chakraborti N, Siva Kumar B, Satish Babu V, Moitra S, Mukhopadhyay A. A new multi-objective genetic algorithm applied to hot-rolling process. Appl Math Model 2008b;32:1781–9. [21] Mishra BK. Study of media mechanics in tumbling mills by the discrete element method. Ph.D. Thesis, Department of Metallurgical Engineering, University of Utah; 1991. [22] Mishra BK, Rajamani RK. Simulation of charge motion in ball mills: numerical simulation. Int J Miner Process 1994;40:171–97. [23] Datta A, Rajamani RK. A direct approach of modelling batch grinding in ball mills using population balance principles and impact energy distribution. Int J Miner Process 2002;64:181–200. [24] Lynch AJ, Rao TC. Modeling and scale up of hydrocyclone classifiers. In: Proceedings XI Int Min Proc Congress, Cagliari; 1975. p. 245–69. [25] Plitt LR. A mathematical model of the hydrocyclone classifier. CIM Bull 1976;69:114–23. [26] Herbst JA, Siddique M, Rajamani K, Sanchez E. Population based approach to ball mill scale-up: bench and pilot scale investigations. Trans AIME 1983;272:1945–54. [27] Kinneberg DJ, Herbst JA. A comparison of linear and nonlinear models for open-circuit ball mill grinding. Int J Miner Process 1984;13:143–65. [28] Mitra K, Gopinath R. Multiobjective optimization of an industrial grinding operation using elitist nondominated sorting genetic algorithm. Chem Eng Sci 2004;59:385–96. [29] Valipour M. Sprinkle and trickle irrigation system design using tapered pipes for pressure loss adjusting. J Agrc Sci 2012;4:125–33. [30] Valipour M, Ali Gholami Sefidkouhi M, Eslamian S. Surface irrigation simulation models: a review. Int J Hydrol Sci Technol 2015;5:51–70. [31] Yannopoulos SI, Lyberatos G, Theodossiou N, Li W, Valipour M, Tamburrino A, et al. Evolution of water lifting devices (pumps) over the centuries worldwide. Water 2015;7:5031–60. [32] Khasraghi MM, Ali Gholami Sefidkouhi M, Valipour M. Simulation of openand closed-end border irrigation systems using SIRMOD. Arch Acker Pfl Boden 2015;61:929–41. [33] Huang X. Chance-constrained programming models for capital budgeting with NPV as fuzzy parameters. J Comput Appl Math 2007;198:149–59. [34] Mitra K, Gudi RD, Patwardhan SC, Sardar G. Midterm supply chain planning under uncertainty: a multi-objective chance constrained programming framework. Ind Eng Chem Res 2008;47:5501–11. [35] Xu J, Zhou X. Approximation based fuzzy multi-objective models with expected objectives and chance constraints: Application to earth-rock work allocation. Inf Sci 2013;238:75–95. [36] Marano GC, Quaranta G. Fuzzy-based robust structural optimization. Int J Solids Struct 2008;45:3544–57. [37] Deb K., Multi-objective optimization using evolutionary algorithms. Chichester, UK: Wiley; 2001. [38] Petzold LR, et al. A description of DASSL: a differential/algebraic system solver. In: Stepleman RS, et al., editors. Scientific Computing. Amsterdam: North-Holland; 1983. p. 65–8. [39] Zadeh LA. Fuzzy sets, Inf. Control 1965;8:338–53. [40] Liu B, Liu YK. Expected value of fuzzy variable and fuzzy expected value models. IEEE Trans Fuzzy Syst 2002;10:445–50. [41] Marano G C, Quaranta G. Robust optimum criteria for tuned mass dampers in fuzzy environments. Appl Soft Comput 2009;9:1232–43.