Enhancing the ϵ-constraint method through the use of objective reduction and random sequences: Application to environmental problems

Computers and Chemical Engineering 87 (2016) 36–48

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Review

Enhancing the -constraint method through the use of objective reduction and random sequences: Application to environmental problems Pedro J. Copado-Méndez a , Carlos Pozo a , Gonzalo Guillén-Gosálbez a,b,∗ , Laureano Jiménez a a b

Departament d’Enginyeria Quimica, Universitat Rovira i Virgili, Avinguda Països Catalans 26, 43007 Tarragona, Spain Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College London, SW7 2AZ, UK

a r t i c l e

i n f o

Article history: Received 15 September 2014 Received in revised form 16 December 2015 Accepted 22 December 2015 Available online 31 December 2015 Keywords: Multi-objective optimization Objective reduction Hypervolume Epsilon-constraint MOO

a b s t r a c t The -constraint method is an algorithm widely used to solve multi-objective optimization (MOO) problems. In this work, we improve this algorithm through its integration with rigorous dimensionality reduction methods and pseudo/quasi-random sequences. Numerical examples show that the enhanced algorithm outperforms the standard -constraint method in terms of quantity and quality of the Pareto points produced by the algorithm. Our approach, which is particularly suited for environmental problems that tend to contain several redundant objectives, allows dealing with complex MOO models with many objectives. © 2015 Elsevier Ltd. All rights reserved.

Contents 1. 2. 3.

4. 5.

6.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Mathematical background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Proposed approach: enhanced -constraint method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1. Dimensionality reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2. Random sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.1. Halton sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.2. Sobol sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.3. Detailed steps of the algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2.4. Remarks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .40 Comparison procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Cases of study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.1. Sustainable planning of ethanol supply chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2. Sustainable planning of hydrogen supply chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Conclusions and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

∗ Corresponding author at: Departament d’Enginyeria Quimica, Universitat Rovira i Virgili, Avinguda Països Catalans 26, 43007 Tarragona, Spain. Tel.: +34 977558618; fax: +34 977558618. E-mail addresses: [email protected] (P.J. Copado-Méndez), [email protected] (C. Pozo), [email protected] (G. Guillén-Gosálbez), [email protected] (L. Jiménez). http://dx.doi.org/10.1016/j.compchemeng.2015.12.016 0098-1354/© 2015 Elsevier Ltd. All rights reserved.

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

1. Introduction Multi-objective optimization (MOO) problems arise in all kinds of industrial application areas ranging from production to service industries, entertainment, and many others. The -constraint method is probably the most widely used approach to solve MOOs. This technique, which was first introduced by Haimes et al. (1971), relies on solving a series of single-objective problems in which one objective is kept in the objective function while the others are transferred to auxiliary constraints that bound them within some allowable levels. Each run of these single-objective problems generates a different solution that is guaranteed to be, at least, weakly Pareto efficient. One of the main advantages of this algorithm is that it can handle non-convex Pareto sets, as oppose to other approaches like the weighted-sum method (Zadeh, 1963) or goal programming (Charnes et al., 1955, 1967; Charnes and Cooper, 1961; Ijiri, 1965). Because of its advantageous characteristics, the -constraint method has found many applications in process systems engineering, system biology, and particularly in problems in which several environmental objectives must be simultaneously optimized (see Grossmann and Guillén-Gosálbez (2010)). The main drawback of the -constraint method is that it is very sensitive to the number of objectives. The standard -constraint algorithm keeps one objective as main criterion and divides the domain of the rest into equal intervals. A set of single-objective problems is then solved for the limits of those intervals, generating in each run a different Pareto point. If we choose p partitions for each objective, we will have to solve 2 + p − 1 problems for the case of two objectives, 3 + (p − 1)2 for three, and so on. In general, k + (p − 1)k−1 , being k the number of objectives being optimized, and p the number of partitions. Note that to determine the limits of the intervals, we first need to optimize each single objective separately. Hence, the complexity of this method grows exponentially in size with the number of objectives, and can get out of hand very easily for problems with several objectives. Besides, a high percentage of sub-problems lead to unfeasible or repeated Pareto solutions, thereby wasting part of the computational effort. For this reason, the overwhelming majority of works in the literature that make use of the -constraint method optimize two objective functions, and three at most, but very seldom more than three. This paper shows that the integration of objective reduction with quasi/pseudo-random sampling in the traditional -constraint method improves the performance of this method for problems with 3 or more objectives. An enhanced -constraint method is therefore introduced that integrates these two ingredients: (i) a rigorous objective reduction technique that eliminates redundant objectives from the search, thus decreasing the complexity of both the generation and analysis of the Pareto solutions; and (ii) the use of pseudo/quasi-random sequences (i.e., uniform distribution and Sobol and Halton sequences), that allow generating in a more efficient manner the values of the epsilon parameters used in the single-objective problems, thereby increasing the percentage of unique feasible solutions obtained. We illustrate the capabilities of our approach through its application to two supply chain design problems, in which we optimize the economic performance along with a set of environmental metrics quantified following LCA principles. Numerical results show that compared to the standard -constraint method the enhanced algorithm provides solutions of better quality in less computational time. Our approach can find many applications in process systems engineering, being particularly suited for environmental problems in which several environmental objectives tend to be redundant. Hence, the main novelty of this work is the development of an enhanced -constraint method that can be used to solve a wide variety of MOO problem arising in process systems engineering. Potential PSE applications could include the optimal design of

37

supply chains (Yue et al., 2014; Guillén-Gosálbez, 2008), refineries (Gebreslassie et al., 2013; Zhang et al., 2014; Santiba nez Aguilar et al., 2014; Murillo-alvarado and El-halwagi, 2013), reverse osmosis networks (Du et al., 2014), chemical plants (Azapagic et al., ˇ cek, 2013) and 1999; Buxton and Pistikopoulos, 2004; Kravanja, Cuˇ the multi-objective optimization of metabolic networks (Pozo and Guillén, 2012; Banga, 2008), among others. The article is organized as follows. In Section 2, we introduce the traditional -constraint method before describing in Section 3 an enhance version of it. Then, in Section 4, we describe how to compare both methods and evaluate their performance. The capabilities of our approach are tested next (Section 5), while in the final section (Section 6) the conclusions of the work are drawn. 2. Mathematical background Let us consider the following MOO problem: minF(x) := (f1 (x), . . ., fk (x)) x

s.t. gj (x) ≤ 0,

j = 1, 2, . . ., m,

hl (x) = 0,

l = 1, 2, . . ., e,

(1)

with k objective functions fi := X −→ R, 1 ≤ i ≤ k, where each objective function fi maps a solution x ∈ X to a value of the function vector F(x) : = (f1 (x), . . ., fk (x)), m is the number of inequality constraints, and e is the number of equality constraints. x ∈ X is a vector of design variables (also called decision variables). Definition 1. x weakly dominates y with respect to the objective set F (x  F y) if (x, y) ∈  F . Definition 2. x * ∈ X is called Pareto optimal if there is no other x ∈ X that dominates x* with respect to the set of all objectives. The -constraint method relies on solving a set of singleobjective auxiliary problems in which one objective is kept as main objective while the others are transferred to auxiliary constraints that bound them within some limits. This algorithm solves singleobjective problems of the following form: min fi (x)

where

1≤i≤k

.. .

.. .

.. .

.. .

s.t. fj (x) ≤ rj

where i = / j

1≤j≤k

1≤r≤p

x ∈ X

where

LBj ≤ rj ≤ UBj

(2)

As already mentioned, this method is quite sensitive to the number of objectives. In general, the number of iterations is equal to k + (p − 1)k−1 , being p the number of partitions defined for every objective and k the number of objectives. It is therefore clear that the computational burden grows rapidly (in fact, exponentially) in size with the number of objectives. In the sections that follow we introduce some improvements that expedite the application of this algorithm. 3. Proposed approach: enhanced -constraint method In this work, the -constraint algorithm is expedited through the integration of two main ingredients: a dimensionality reduction method and the use of pseudo/quasi-random sequences. Each of these techniques is next explained in detail before providing a detailed description of the whole algorithm.

38

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

(a)

(b)

Fig. 1. Illustrative example of objective reduction, where we consider 3 Pareto optimal solutions that minimize 4 different objectives: (1, . . ., 4). As observed in (a), objectives 2 and 3 can be omitted (orange rectangles) without altering the dominance structure. This is because xf1 ,f4 y is satisfied if and only if xf1 ,f2 ,f3 ,f4 y is satisfied. On the other hand, it can be observed in (b) that further reductions are not possible without modifying the dominance structure. In this case, the delta value is 0.5 because it is the maximum difference between solutions x1 and x2 is in the objective f4 . Note that, in order to compute delta value only omitted objectives are considered. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

3.1. Dimensionality reduction Dimensionality reduction methods eliminate redundant objectives in MOO problems (see Fig. 1). Consider the general MOO model introduced before. The goal of dimensionality reduction algorithms is to find a subset F of objectives functions pertaining to the original set F with the following property: when we optimize the problem in the reduced domain of objectives F (rather than in the original domain F), we will generate (in less CPU time) a Pareto front that is very close to the “true” Pareto front of the original problem. Different dimensionality reduction strategies have been proposed in the literature. In a seminar work, Deb and Saxena (2005) introduced an algorithm for objective reduction based on principal component analysis. An alternative algorithm to reduce the number of objectives that is based on a feature selection technique was introduced by López Jaimes et al. (2008). Brockhoff and Zitzler were the first to formally state the problem of reducing the number of objectives in MOO. Their seminar work introduced the concept of error of the approximation obtained after removing objectives in an MOO problem, and formally stated the following two problems: (1) computing a minimum objective subset (MOSS) of a multi-objective problem that does not exceed a certain approximation error (denoted as the ı-MOSS problem); and (2) identifying a minimum objective subset of size k with minimum approximation error (k-MOSS problem). The authors presented also an exact and an approximation algorithm to tackle these problems that were applied to well-known test case studies, showing that substantial dimensionality reductions are possible while keeping the approximation error low (Brockhoff and Zitzler, 2006). Bearing these ideas in mind, Guillén-Gosálbez (2011) introduced an approach for dimensionality reduction based on a mixed-integer linear program (MILP) that solves both the k-MOSS and ı-MOSS problems, and takes advantage of the powerful branch-and-cut algorithms available for MILP (Guillén-Gosálbez, 2011; CopadoMéndez et al., 2014). More recently, Thoai (2011) proposed ways to reduce the number of objectives and dimension of a linear multiple criteria optimization problem using the concept of so-called representative and extreme criteria (Thoai, 2011). In this paper, we use the dimensionality reduction method introduced by Guillén-Gosálbez (2011) and later enhanced by Copado-Méndez et al. (2014), to expedite the performance of the -constraint method. We next provide a very simple example of dimensionality reduction for illustrative purposes. Further details of this method can be found in the original publication. Consider the 3 Pareto optimal solutions depicted in Fig. 1 that minimize 4 different objectives: (f1 , f2 , f3 , f4 ). This figure is a

parallel coordinates plot (Purshouse and Fleming, 2003) that depicts in the x axis the set of objectives and in the y axis the normalized value attained by each solution. Each line in the figure represents a different Pareto solution. As seen, all these lines intersect in at least one point, as no solution is dominated by any of the others. As observed, two objectives can be omitted (i.e., objective 2 and 3) without changing the dominance structure. This is because xf1 ,f4 y is satisfied if and only if x  f1,f2,f3,f4 y is satisfied. Further reductions are not possible without modifying the dominance structure. For instance, if we remove f1 and f4 , then x1 f2 ,f3 x2 is satisfied although x1  / f1 ,f2 ,f3 ,f4 x2 . As seen, solution x2 would dominate x1 in the original 4-dimensional Pareto space {f1 , f2 , f3 , f4 } if it showed the same value of f4 than x1 . The difference between the true value of f4 in x2 and that required to dominate x1 in the original space of objectives can be used as a measure to quantify the change in the dominance structure. For this example, this value, referred to as ı value by Brockhoff and Zitzler, is 0.5. The delta value quantifies the change in the dominance structure of a MOO problem that takes place after removing objectives. Our goal is to determine subsets of objectives of a MOO model that minimize the ı value. This problem was formally defined in a pioneering work by Brockhoff and Zitzler (2006a). Alternatively, we may also be interested in calculating the minimum number of objectives for a given allowable approximation error. We refer to these problems as ı-MOSS and k-MOSS problems, respectively. As it was said before, in this paper we use our algorithmic framework for dimensionality reduction developed in previous publications, which can solve efficiently the two aforementioned problems.

3.2. Random sequences Objective-reduction techniques might be unable to eliminate a significant amount of objectives in cases where few redundant criteria exist. If we partition the domain of each objective into equal intervals during the application of the -constrain method, we might still need to solve a very large number of single-objective models. Moreover, most of them might be either unfeasible or produce the same Pareto solution. To overcome this limitation to the extent possible, we propose to generate the epsilon parameters using random sequences. Particularly, in this paper we focus on the following strategies: a pseudo-random sequence (uniform distribution), and the Halton and Sobol sequences (quasi-random or low-discrepancy). In the context of our work, the aforementioned sequences are used to generate the values of the epsilon parameters in a more effective manner. The goal is to spread these values as

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

39

Fig. 2. Illustrative example of how to define discrepancy, which gives an approximation error of the area of the union of these rectangles. The area of union of the rectangles is known (0.44), yet it can approximated by means of counting how many points amongst 100 fall into the rectangles. When these points are distributed equidistantly, we 61 39 − 0.44| = 0.17, but when they are distributed according a pseudo-random sequence, discrepancy decreases to | 100 − 0.44| = 0.05, and finally obtain a discrepancy of | 100 44 −0.44| = 0. when these points are distributed according to Halton sequence, the discrepancy computed is | 100

much as possible so they cover a wider region of the search space, which will eventually lead to a better Pareto front. The pseudo-random sequence is well known, and therefore we will not get deeper into its details. The Halton and Sobol sequences (also called low-discrepancy or quasi-random sequences) are useful in numerical integration, as well as in simulation and optimization. Surveys of applications of low-discrepancy sequences can be found in Fox (1986), Bratley and Fox (1988), Bratley et al. (1992). The concept of low-discrepancy, which is introduced next, will be used in the ensuing sections for justifying the use of the sequences: The discrepancy Dpk of point set {xi }i=1,...,p ∈ [0, 1)k

Definition 3. is:

  A(E, p)

Dpk = sup  E

p

 

− (E)

therefore likely to lead to a better approximation of the Pareto set when used in the context of the epsilon constraint method (Uhlmann, 1965). We next describe each of the sequences used in our work in detail. 3.2.1. Halton sequence The Halton sequence, which was introduced by Halton (1960), can be used like a random number generator to produce points in the interval [0, 1]. The standard approach shown below is employed in our work. More details on this strategy can be found in Faure and Lemieux (2010). In one dimension, the standard Halton sequence is generated by choosing a prime number r (r ≥ 2), and decomposing an integer g in terms of the base r:

(3)

where E = [0, t1 ) × · · · × [0, tk ), 0 ≤ tj ≤ 1, j = 1, . . . k, (E) is the hypervolume of E, A(E, p) is the number of xi contained in E, p is number of points and, k is the number of dimensions. In other words, E is a set of hyperectangles, (E) is the hypervolume, A(E,p) is the percentage of points xi , that fall in E, and p discrepancy is the largest difference between A(E,p) and (E), that p is, the discrepancy measures the error in the hypervolume estimation (see Fig. 2). Let us assume that the area of rectangles in Fig. 2 represents a Pareto set and that our goal is to approximate this area counting the points contained in it. The equidistant sampling might not perform well because it produces points that are not necessary for computing the area. Similarly, the -constraint method using equidistant sampling produces repeated solutions. In contrast, the pseudo-random and low-discrepancy sequences are more efficient, as seen in Fig. 2, thereby leading to errors below that of the equidistant approach. Definition 4. A sequence of p points is considered a lowdiscrepancy sequence, if its Dk (p) is in O(log(p)k ).1 A sequence of p points of dimension k which are distributed in a set of hyperrectangles E approximates the hypervolume with a given error. This error is called discrepancy of p and if it is at most log(p)k , then the sequence p is a low-discrepancy sequence. Pseudo/quasi-random sequences are expected to be more efficient than the simple equidistant sequence, as approximating Pareto sets is analogous to the widely known numerical integration of a function (see Uhlmann (1965)). The rate of convergence of pseudo/quasi-random sequences as applied to numerical integration is better than that of equidistant sequence, and they are

1 A function f(n) is in O(g(n)) if f(n) grows at most at the same order of g(n). e.g. n and n2 are in O(n2 ).

g=

L 

bl (g)r l

(4)

l=0

where 0 ≤ bl (g) ≤ r − 1 and rL ≤ g ≤ rL+1 . L is the length of the integer g in base r. The Halton value ϕ(g) is generated by multiplying each digit in base r of g in reverse order by the power of r−l−1 : ϕ(g) =

L 

bl (g)r −l−1

(5)

l=0

where 0 . b0 (g). . .bL (g) are digits of g in base r in reverse. For instance, suppose we want to obtain the fourth element (g = 4) of the Halton sequence in base 2. The integer 4 is written in base 2 as 4 = 1 *22 + 0 *21 + 0 *20 . According to Eq. (5), this number in reverse is 0.001, and then 0 ∗ 2−1 + 0 ∗ 2−2 + 1 ∗ 2−3 = 18 . Therefore 1 forms the fourth number of the Halton sequence. 8 The Halton sequence for k dimensions is generated by the method shown above, but taking k first prime numbers as base. 3.2.2. Sobol sequence The other low-discrepancy sequence applied here is the Sobol sequence, which was introduced by Sobol (1967). We provide an informal description below, while theoretical aspects can be found in Sobol (1967) and Bratley and Fox (1988). For simplicity, we assume the case of one dimension in which we aim to generate the g Sobol value ϕ(g), with low discrepancy over the unit interval. To begin, we need to obtain the set of direction numbers v1 , v2 , . . ., where each vi is a binary fraction with this form vi = 0.vi1 vi2 . . .. Alternatively, vi = mii , where mi is an odd integer, 2

0 < mi < 2i . The direction numbers come from a primitive polynomial belonging to the finite field Z2 ,2 which presents the form

2

The field Z2 is the finite set {0,1}.

40

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

P ≡ xd + a1 ∗ x1d−1 + . . . + ad−1 ∗ x + 1 where ai ∈ {0, 1} and P is a primitive polynomial of degree d in Z2 . The polynomial coefficients are used to define the recurrence for calculating vi , thus: vi =a1 ∗ vi−1 ⊕ a2 ∗ vi−2 ⊕ . . . ⊕ ad−1 ∗ vi−d+1 ⊕ vi−d ⊕

v

i−d 2d

 , i>d (6)

where ⊕ denotes exclusive-or operator. Therefore, this recurrence can be rewritten in terms of mi mi = 2 ∗ a1 ∗ mi−1 ⊕ 22 ∗ a2 ∗ vi−2 ⊕ . . . ⊕ 2d−1 ∗ ad−1 ∗ mi−d+1 ⊕ 2d ∗vi−d ⊕ vi−d ,

i>d

(7)

heuristic approach consisting of solving a set of bi-criteria problems in each of which we trade-off one objective (when dealing with environmental problems, we will choose the economic performance), against each of the remaining criteria (the environmental indicators) separately. For each of these bi-criteria problems, the -constraint method is run for a given number of iterations, producing a set of Pareto solutions. The objective reduction method is then applied to identify and eliminate redundant objectives. A set of values of the epsilon parameters are next generated using the pseudo-random sequences. The single-objective problems are finally solved for these parameters values. Hence, the detailed steps of the enhanced -constraint algorithm are shown in Fig. 3. 3.2.4. Remarks

x3 + x + 1

• The outcome of the objective-reduction algorithm depends on the number of Pareto points used in the analysis. Experimental results generated for different problems indicate that a set of Pareto points of small/medium size suffices to identify the redundant objectives. • The general approach presented in this article can be easily extended to work with other MOO solution algorithms. • It is possible to define an outer loop in the calculations in order to repeat the execution of the dimensionality reduction algorithm with more Pareto points as iterations proceed. In practice, we found that such a loop does not improve significantly the performance of the overall method, since the algorithm for dimensionality reduction provides consistent results for a small number of Pareto points.

then

4. Comparison procedure

where ⊕ denotes exclusive-or operator. The g Sobol value ϕ(g) is obtained as follows: ϕ(g) = ϕ(g − 1) ⊕ vc

(8)

where c is right-most zero bit of Gray code representation of g. The Gray code is a binary numeral system where each number differs from the previous one in only one bit, and it can be computed as follows: g = . . .g3 g2 g1 := . . .b3 b2 b1 ⊕ . . .b4 b3 b2

(9)

To start the recurrence, we take ϕ0 = 0. For example, we take as primitive polynomial

mi = 2 ∗ mi−2 ⊕ 8 ∗ mi−3 ⊕ mi−3

In this section we describe the different approaches tested and how we assessed them. Particularly, the following approaches are tested in the case studies (see Fig. 4):

and take m1 = 1, m2 = 3 and m3 = 7. Therefore m4 = 12 ⊕ 8 ⊕ 1 = 1100 ⊕ 1000 ⊕ 0001 = 0101 = 5, v1 =

1 =binary 0.1 21

v2 =

3 =binary 0.11 22

v3 =

7 =binary 0.111 23

v4 =

5 =binary 0.0101 24

ϕ(g(1)) = ϕ(g(0)) ⊕ v1 = 0 ⊕ 0.1 = 0.1=decimal

1 2

• • • • •

Standard -constraint (equidistant epsilon values, see Fig. 5). Standard -constraint with pseudo-random epsilon values. Standard -constraint with Halton epsilon values. Standard -constraint with Sobol epsilon values. -constraint integrated with objective reduction and using equidistant epsilon values (see Fig. 5). • Our -constraint integrated with objective reduction and using pseudo-random epsilon values. • Our -constraint integrated with objective reduction and using Halton epsilon values. • Our -constraint integrated with objective reduction and using Sobol epsilon values.

(c = 1 and g(1) = 0001

ϕ(g(2)) = ϕ(g(1)) ⊕ v2 = 0.1 ⊕ 0.11 = 0.01=decimal

1 4

(c = 2 and g(2) = 0011)

ϕ(g(3)) = ϕ(g(2)) ⊕ v1 = 0.01 ⊕ 0.1 = 0.11=decimal

3 4

(c = 1 and g(3) = 0010)

ϕ(g(4)) = ϕ(g(3)) ⊕ v3 = 0.01 ⊕ 0.111 = 0.101=decimal

5 (c = 3 and g(4) = 0110) 8

And so on. This method can be extended to k dimensions by choosing k different primitive polynomials, computing their direction numbers and generating each component ϕ(g, k) separately. 3.2.3. Detailed steps of the algorithm Having presented the ingredients of our algorithm, we next describe in detail how it works. The first step is to generate an initial set of Pareto points on the basis of which we will perform the objective-reduction analysis. To this end, we apply a

In order to assess the performance of each algorithm, we used two different indicators that quantify the quality of the Pareto front produced by them. The first is the number of unique Pareto solutions generated in a given time. The second is the quality of the Pareto front, which is quantified using the hypervolume indicator. The hypervolume indicator (or S-metric, Lebesgue measure), which was introduced by Zitzler and Thiele (1998), is regarded as a fair measure of the quality of a Pareto front, given its convenient mathematical properties (Zitzler et al., 2003). Formally, the hypervolume indicator is defined as follows: Definition 5. Given a set P of Pareto points , the hypervolume indicator is the size of the hyperrectangle d where d = {x ∈ Rd+ :   x for some  ∈ P} This hyperrectangle corresponds to the space which is dominated by at least one point in the set P. The dominated hypervolume

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

41

Fig. 3. Algorithm for the proposed methodology. Note that steps 1b and 2b make use of the pseudo/quasi-random sequences to generate the required epsilon parameters in an efficient, whereas step 1c resorts to an objective reduction technique to eliminate redundant objectives that will only waste computational time and complicate the interpretation of the results if retained. These are the two main novelties of the enhanced -constraint method introduced in this manuscript.

is calculated with respect to a reference point which is chosen to coincide with the nadir point. In summary, the hypervolume indicator measures the area dominated by the Pareto optimal solutions and therefore reflects the dominance of sets. Besides, the hypervolume also promotes diverse sets. Hence, this indicator unifies the two criteria: dominance and diversity (see Fig. 6). The calculation of the hypervolume indicator is a particular case of the Klee’s measure problem, which was formulated by Klee (1977). The best algorithm currently known for a d-dimensional d

space is O(N log(N + N 2 )), and was obtained by Beume (2006). The complexity of computing the hypervolume grows rapidly in size with the number of objectives and Pareto points. For our experiments, the approximation algorithm developed by Everson et al. (2002) and implemented in Matlab was used. This approach uses Monte Carlo sampling to approximate the hypervolume, which avoids expensive calculations. The algorithm employed to measure the hypervolume indicator defines a hyperrectangle for each Pareto point, where each hyperrectangle has as lower coordinate the Pareto point and as upper coordinate the reference point.

The reference point is the same for all the hyperectangles (see Fig. 6). In essence, the algorithm distributes ˛ random points in the space and counts how many points are contained within the hyperrectagles, which corresponds to the number of dominated points. Finally, the hypervolume indicator is calculated as the quotient between the number of dominated points and ˛. Note that the accuracy of this method depends on the number of points distributed ˛, and for this reason ˛ is an input parameter of the algorithm. Obviously, more points lead to more accuracy, but also to larger CPU times. A brief outline of this algorithm is given in Fig. 6 5. Cases of study We test the capabilities of our method through its application to two supply chain design problems, where we optimize the economic performance against a set of life cycle assessment metrics. These problems are formulated as multi-objective MILPs. Details on them can be found in Kostin et al. (2012) and Sabio et al. (2010, 2012). The case studies that we deal with are based on a superstructure of alternatives that accounts for a set of available technologies to produce, store and deliver ethanol/hydrogen (see Fig. 7).

42

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

Fig. 4. This is the procedure followed to analyze the performance between standard -constraint and our approach. In step 1, our approach was run for each sampling techniques (steps 1a–1c) and making use of the objective redaction technique. In step 2, -constraint method was run for each sampling techniques but this time the objective reduction algorithm was not applied. In steps 3 and 4, we carried out the feasibility analysis and the computation of hypervolume indicator.

Fig. 5. The pseudo-code for the equidistant sampling method is shown in (a). Roughly speaking, the algorithm requires itmax as input parameter which specifies the maximum number of iterations. In line 2 the algorithm chooses a word as the first one. In the main loop while the end condition is not satisfied, the algorithm solves the related constraint problem (line 3a) and computes the successor word (line 3b) as long as the end condition is not satisfied. An example of words generation in lexicographic order is shown in (b). Given an alphabet , each letter is assigned to an interior point of the objective intervals. In the example of figure (b), each objective is divided into 4 partitions giving rise to 3 interior points. Hence, 9 words (||k−1 = 33 = 9) are generated from aaa to ccc. The word bbc corresponds to the second interior points of objectives 2 and 3 (2b and b3 ) and the third interior point of the objective 4 (4c ).

Two multi-objective MILPs were derived, one for each example, and implemented in GAMS, where the objective reduction algorithm was also coded. The MILPs were solved by CPLEX 12.5.1.0 on an Ubuntu 13.10 with an Intel i7-4770, 3.4 GHz 8-cores processor, and 16 GB of RAM. The Halton and Sobol sequences were implemented in GAMS as an extrinsic function. The objective reduction was carried out using the algorithm presented in a previous publication (see Section 3.1). To run this algorithm, we used a Pareto set generated by solving a set of bi-objective problems in each of which we applied pseudo-random sequences to generate the values of the epsilon parameters as explained in Section 4.

5.1. Sustainable planning of ethanol supply chain The first example determines the configuration of a threeechelon bioethanol network and associated planning decisions that maximize the net present value and minimize the environmental impact. Decisions to be made include the number, location, and capacity of the production plants, and warehouses to be set up in each region, their capacity expansion policy for a given forecast of prices and demand over the planning horizon, the transportation links and transportation modes that need to be established in the network, and the production rates and flows of feedstocks, wastes, and final products (Kostin et al., 2012).

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

43

Fig. 6. The pseudo-code for computing the hypervolume indicator is outlined in (a): in short, the algorithm distributes ˛ random points in the space (line 1) and counts how many points are contained inside of the rectangles (dominated points, line 2). Finally, the hypervolume indicator is calculated by dividing the number of dominated points by ˛ (line 3). An illustrative example of the hypervolume computation is shown in (b): given 3 Pareto points (diamond points), the algorithm considers the rectangles associated with each Pareto point (blue, red and orange rectangles), where each rectangle has as lower coordinate the Pareto point and as upper coordinate the reference point. The reference point is the same for all the rectangles. The hypervolume indicator in this example is 0.42, namely 42 dominated points divided by 100. Note that the accuracy of this computation depends on ˛ random points: the greater ˛ is, the more accuracy this method achieves. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 7. Superstructure for the supply chain design problem, which is formulated as an MILP. We derive an MILP for the optimal design of the network that optimizes simultaneously the net present value/total cost and environmental performance. The MILP provides the optimal location of the production and storage facilities, the technologies to be implemented and the transportation links between the SC entities.

The multi-objective optimization problem contains 6 objectives: the net present value (economic objective and abbreviated as NPV), and the individual categories considered in the Eco-indicator 99, namely, damage to human health (DHH), damage to eco-system quality (DTE), and damage to resources (DTR), along with the global warming potential (GWP), and the Eco-indicator 99 itself (EI99).

The MILP contains 36,629 continuous variables, 10,962 discrete variables, and 48,588 equations. Following our approach (see Figs. 3–6), we first generated 100 Pareto solutions solving bi-criteria problems in the original search space. More precisely, 20 Pareto solutions were generated with the bi-objective model NPV-DTE, other 20 with the bi-objective

44

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

Fig. 8. Feasibility analysis and hypervolume for the case of sustainable planning of ethanol supply chains in each sampling techniques for 6 objectives and 4 objectives: (a) number of unique feasible solutions vs number of iterations. (b) Number of unique feasible solutions vs time. (c) Hypervolume vs CPU time.

model NPV-DTR, and so on. We then applied the objective reduction (ı-MOSS) algorithm (Copado-Méndez et al., 2014), which took 2 s to identify 2 redundant objectives (DHH, EI99) considering an approximation error ı equal to 0. We next ran 2000 iterations of the enhanced -constraint method as well as the standard -constraint method using each sampling technique (equidistant, random, Halton and Sobol). The equidistant sampling was carried out by optimizing the NPV against the different environmental objectives considering 5 partitions for each objective. The 4 interior partition points of each objective j were labeled with alphabet letters aj . . . dj . Each subproblem was then defined as a permutation with repeats of the letters a . . . d in lexicographic order3 (see Fig. 5). As observed in Fig. 5, firstly, we choose a random permutation as seed, and then, we solve the next 2000 sub-problems considering the environmental objectives in the order discussed above. The methods were applied to both, the original full space problem and the reduced model without redundant objectives. The results are shown in Fig. 8 and Tables 1 and 2. Particularly, in Table 1 we show the number of unique feasible solutions, the total time required to run 2000 iterations and the hypervolume obtained for each approach, whereas in Table 2 some intermediate results (i.e., at 2000 s) are shown. Note that the results in both Tables (1 and 2) have been sorted in first place according to the hypervolume indicator and in second place by the number of unique feasible solutions (both in descending order). As seen in Tables 1 and 2 and Fig. 8, the best approach according to the hypervolume obtained after the 2000 iterations is the constraint method with Halton sampling run in the reduced space (RS), followed by the Sobol and random methods executed in the

3

Lexicographic order is the same as dictionary order.

Table 1 Final results obtained for each approach in ethanol case. 2000 iterations

# solutions

Total time (s)

Hypervolume

Halton RS Sobol RS Random RS Sobol FS Halton FS Equidistant RS Random FS Equidistant FS

936 943 925 557 531 127 537 69

5.45 × 103 5.43 × 103 5.40 × 103 4.93 × 103 4.58 × 103 9.74 × 103 4.61 × 103 12.43 × 103

0.16 0.15 0.15 0.05 0.05 0.05 0.04 0.03

Table 2 Intermediate results computed up to 2000 s in ethanol case. 2000 s

# solutions

Hypervolume

Halton RS Sobol RS Random RS Sobol FS Halton FS Equidistant RS Random FS Equidistant FS

340 347 341 233 222 48 232 38

0.15 0.14 0.14 0.05 0.05 0.05 0.04 0.03

RS, the Sobol and Halton run in the full space (FS), the equidistant sampling in the RS, the random sampling in the FS and the equidistant sampling in the FS. We next proceed to analyze in detail these results: • Full space (FS) against reduced space (RS): the -constraint method run in the RS solves single-objectives problems with less auxiliary constraints than the FS approach. This leads to significant CPU time reductions for each iteration, as seen in Fig. 8 and in the column “Total time” of Table 1. Having fewer constrained

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

objectives reduces the chances of having one constraint acting as a bottleneck, which prevents the algorithm from calculating repeatedly the same solution over the 2000 iterations. This increases significantly the amount of new unique feasible solutions, having in turn a positive effect on the hypervolume indicator. Note that the hypervolume value gained by reducing objectives is higher than that produced by using only pseudo/quasi-sampling sequences (see Tables 1 and 2). • Equidistant sampling against pseudo/quasi-random sampling: as we observe in Tables 1 and 2 and Fig. 8, the equidistant sampling is the slowest even in the RS. A possible reason why this might happen is that solvers firstly try to identify solutions involving active constraints. Defining equidistant values for the epsilon parameters makes it more difficult to find feasible integer solutions satisfying the constraints as equalities (Ehrgott and Ryan, 2003). In practice, the branch and bound method needs to perform a large number of iterations until it can guarantee that the problem is unfeasible. On the other hand, when the epsilon parameters values are chosen by pseudo/quasi-random sequences, it is more likely that there will be an integer solution satisfying the auxiliary epsilon constraint as a strict equality (Ehrgott and Ryan, 2003). Note that, even in the FS, the number of unique feasible solutions obtained with pseudo/quasi-random sequences is greater than that obtained with the traditional (i.e., equidistant) -constraint in the RS. • Pseudo-random against quasi-random sampling: Halton and Sobol sequences (see Section 3.2) choose epsilon parameters values with random properties, but more uniformly distributed than pseudo-random sequences. This is because the quasi-random sequences choose the next parameter value taking into account the previous one, and because of this not all the values have the same probability of occurrence. On the other hand, the epsilon parameters chosen by pseudo-random sequences are all equiprobable. If we compare the percentage of epsilon parameters values giving rise to Pareto frontier points, we will see that it is higher for quasi-random sampling than for pseudo-random sampling, as we can observe in Tables 1 and 2, where quasirandom sampling identify more feasible solutions. • Number of feasible solutions against hypervolume indicator: the enhanced -constraint method obtained more unique feasible solutions, yielding in turn a Pareto front with a better hypervolume value (see Tables 1 and 2 and Fig. 8). Note however, that the hypervolume value does not depend exclusively on the number of Pareto solutions, but rather on their quality (i.e., the fraction of the space that they cover). Therefore, a few well situated solutions may involve a higher hypervolume than lots of them located in the wrong place. This is reflected in Tables 1 and 2, where the Sobol RS method shows more unique solutions than the Halton RS, yet the Halton RS yields a better hypervolume (Tables 1 and 2). In summary, we find that in this case study, both strategies (i.e., reducing objectives and using pseudo-random sampling) improve the performance of the algorithm from the viewpoints of number of feasible solutions generated and quality of the Pareto front. Note that, as already mentioned, in this example there are 2 redundant objectives that can be eliminated. This is because some environmental impacts tend to be equivalent (when one is minimized the other one is reduced and vice-versa). A detailed discussion of this topic can be found in Kostin et al. (2012). 5.2. Sustainable planning of hydrogen supply chains This example deals with the optimal design of a hydrogen supply chain for vehicle use in Spain taking into account economic and environmental concerns. The problem, which was first proposed by Almansoori and Shah (2009), considers different technologies for

45

Table 3 Final results obtained for each approach in hydrogen case. 300 iterations

# solutions

Total time (s)

Hypervolume

Sobol RS Sobol FS Random RS Halton RS Random FS Equidistant RS Halton FS Equidistant FS

48 203 59 53 204 24 300 0

10.60 × 102 533.00 × 102 16.80 × 102 12.30 × 102 529.00 × 102 5.62 × 102 538.00 × 102 6.10 × 102

0.955 0.948 0.946 0.946 0.944 0.939 0.936 0

production, storage and transportation of hydrogen to be established in a set of geographical regions distributed all over the country. The goal is to determine the optimal network configuration in terms of its economic and environmental performance. The problem can be formulated as a multi-objective MILP that seeks to minimize the total cost of the network and its environmental impact. In this formulation, integer variables indicate the number of plants and storage facilities to be opened in a specific region (i.e., grid), whereas binary variables are employed to denote the existence of transportation links connecting the SC entities. The environmental impact was quantified via 15 life cycle assessment indicators based on the Eco-indicator 99 methodology (Hischier et al., 2010). The objectives considered are the following: total cost, acidification potential, climate change, eutrophication potential, freshwater aquatic ecotoxicity, freshwater sediment ecotoxicity, human toxicity, ionizing radiation, land use, malodors air emissions, marine aquatic ecotoxicity, marine sediment ecotoxicity, photochemical oxidation, resources, stratospheric ozone depletion and terrestrial ecotoxicity (Hischier et al., 2010). Further details on this case study can be found in a previous publication (Sabio et al., 2010, 2012). The MILP contains 11,804 continuous variables, 7304 discrete variables and 29,495 equations. We proceeded in the same manner as before (see Figs. 3–6). That is, we first generated 100 Pareto solutions using a bi-criteria algorithm in the original search space. A total of 15 bi-objective models were constructed by combining the cost and each environmental objective. From these solutions we selected 100 randomly. The objective reduction method (Copado-Méndez et al., 2014) was run considering 100 Pareto points. It identified 13 redundant objectives (keeping finally 3 objectives: the cost, the photochemical oxidation and terrestrial ecotoxicity) and requiring 4 s for this task. The full-space equidistant sampling was carried out by optimizing the cost against the environmental objectives considering 3 partitions for each objective. The 2 interior partition points of each objective j are labeled as point aj and point bj . According to Fig. 5, an initial sequence is chosen as seed, then from this seed (i.e. abb . . . ab), we solve the 300 sub-problems in lexicographic order considering the environmental objectives in the order presented above. In this case study, we only ran 300 iterations in contrast to the 2000 iterations we executed in the previous case, since in this case each iteration requires higher CPU times. Hence, we ran 300 iterations of the enhanced and standard -constraint methods using each sampling technique (pseudo-random, Halton and Sobol). Fig. 9, which is equivalent to Fig. 8, summarizes the results obtained. Note that Table 3 is equivalent to Table 1 and Table 4 is equivalent to Table 2. It can be seen in Table 3 and Fig. 9 that the best approach, in terms of hypervolume indicator, after 300 iterations is the -constraint method with Sobol sampling run in the reduced space (RS), followed by the Sobol in the full space (FS), random in the RS, Halton in the RS, random in the FS, equidistant in the RS, Halton in the FS and finally equidistant in the FS. We next proceed to analyze the results in further detail.

46

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

Fig. 9. Feasibility analysis and hypervolume for the case of optimal design of a hydrogen supply chains in each sampling techniques for 16 objectives and 3 objectives. (a) Number of unique feasible solutions vs number of iterations. (b) Number of unique feasible solutions vs time. (c) Hypervolume vs CPU time. -Constraint was not able to obtain feasible solutions with equidistant technique in full space.

Table 4 Intermediate results computed up to 500 s in hydrogen case. 500 s

# solutions

Hypervolume

Sobol RS Random RS Halton RS Equidistant RS Sobol FS Halton FS Random FS Equidistant FS

18 22 15 24 3 2 2 0

0.949 0.942 0.941 0.939 0.656 0.557 0.396 0

• Full space (FS) against reduced space (RS): in this case, during the first 500 s, all the RS approaches outperform the FS methods in terms of both feasible solutions and hypervolume, as seen in Table 4. Besides, the RS approaches are able to perform the 300 iterations in a fraction of the CPU time required by the FS methods, as shown in Tables 3 and 4. This might be due to the fact that the RS models show fewer auxiliary constraints. However, consdering the 300 iterations, the RS approaches identify fewer unique feasible solutions than the corresponding FS counterparts as shown in Table 3. This happens because the objective reduction in this case study is very strong, going from 16 to only 3 criteria, which during the optimization gives rise to many repeated solutions and fewer unique feasible solutions. The only exception for this tendency is the -constraint with equidistant sampling, which in the FS was unable to obtain a single feasible solution. Note that running the FS equidistant method considering 3 partitions for every objective would lead to 215 iterations (out of which we only performed 300). Because of this, the probability to obtain a feasible solution in these 300 iterations is extremely low. • Equidistant sampling against pseudo/quasi-random sampling: the -constraint with pseudo/quasi-random sampling

outperforms the equidistant sampling in terms of hypervolume indicator and number of feasible solutions, as seen in Table 4. The reason why this happens is the same as in the previous case. That is, solvers firstly try to identify solutions involving active constraints. Hence, defining the epsilon parameters by pseudo/quasi-random sequences increases the chances to find an integer solution satisfying the epsilon constraints as strict equalities, thereby expediting the solver execution (Ehrgott and Ryan, 2003). On the other hand, when epsilon parameters are chosen by the equidistant method, it is more difficult to find integer solutions for which the epsilon constraint is active. Remarkably, in Table 3, the Halton sequence in the FS is the worst in terms of hypervolume, yet it identifies the largest number of feasible solutions. • Pseudo-random against quasi-random sampling: as seen in Table 3, the Sobol sequence outperforms the pseudo-random method in both, number of unique solutions and hypervolume indicator. We already explained in the previous case study the advantages of quasi-random sequences over the pseudo-random sequence. Note, however, that the Halton sequence did not work so well in this case (i.e., it obtained the lowest hypervolume), while in the FS it identified the highest number of feasible solutions. Recall that obtaining many solutions does not necessarily guarantee a better hypervolume value, as shown in Tables 3 and 4, where it can be seen that the Halton FS identified more solutions than the other methods, despite leading to the second lowest hypervolume value. • Number of feasible solutions against hypervolume indicator: in the first 500 s, the enhanced -constraint with pseudo/quasirandom sampling in the RS produced more unique feasible solutions, leading to higher hypervolume values than the FS methods (see Table 4). On the contrary, if more time is allowed so that the 300 iterations can be completed, then

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

some pseudo/quasi-random approaches in the FS identify more unique solutions than the corresponding RS approaches (yet their hypervolume is still lower). Hence, there is no evidence that a higher number of solutions necessarily lead to a higher hypervolume. For instance, the Halton sequence in the FS identifies the highest number of feasible solutions but leads to the second lowest hypervolume (see Table 3). Finally, we find that the best approach in terms of the hypervolume indicator obtained after 300 iterations is the Sobol in the RS, followed by the Sobol in the FS. In terms of CPU time, we would choose the pseudo-random or Halton methods in the RS instead of the Sobol in the FS. In this second case study, the performance patterns observed in the first case are not so well defined, probably due to numerical issues, yet in general the same tendencies are still observed. In both cases, the use of objective reduction or/and pseudo/quasi-random sampling leads to significant improvements when compared to the standard -constraint method, mainly in terms of CPU time and quality of the Pareto set. Note that a detailed discussion on the existence of redundant environmental indicators in this problem can be found in Sabio et al. (2012). 6. Conclusions and future work This work introduced an enhanced -constraint method that incorporates two main ingredients: an objective reduction technique, and the use of pseudo/quasi-random sampling methods for generating the values of the epsilon parameters. To test the capabilities of our approach, we applied it to the multi-objective optimization of supply chains considering both, economic and environmental concerns. For this comparison, we used several sampling techniques (pseudo-random, Halton and Sobol, and the standard equidistant) combined with an objective reduction algorithm. The performance of each method was assessed in terms of number of unique feasible solutions and hypervolume indicator values. The numerical experiments reveal that the best sampling techniques are the quasi-random sequences (Halton or Sobol), since they improve the performance from the viewpoints of number of feasible solutions and hypervolume value attained in a given time. This is due to the property of low-discrepancy, which distributes the epsilon parameters values in a more efficient manner compared to the pseudo-random and equidistant techniques. Concerning the objective reduction approach, the results show that it expedites the -constraint method by eliminating auxiliary epsilon constraints. Furthermore, this strategy provides valuable insight into the physical problem by identifying redundant objectives and facilitates in turn the post optimal analysis of the solutions by concentrating only on a reduced set of non-redundant objectives (Pozo and Guillén, 2012). Finally, we conclude that combining quasi-random sampling methods with objective reduction improves in general the overall performance of the standalone -constraint method in terms of number of feasible solutions and quality (i.e. hypervolume value) of the Pareto front obtained in a given time frame. Our enhanced method can find many applications in a wide variety of multi-objective problems arising in many domains of science and engineering. Acknowledgments Pedro J. Copado wishes to acknowledge support of this research work from the Spanish Ministry of Education and Science (DPI200804099/DPI). The authors also wish to acknowledge support from the Spanish Ministry of Education and Science (Projects CTQ201237039-C02, DPI2012-37154-C02-02 and ENE2011-28269-C03-03).

47

References Almansoori A, Shah N. Design and operation of a future hydrogen supply chain: multi-period model. Int J Hydrogen Energy 2009];34(19):7883–97. Azapagic A, Clift R. The application of life cycle assessment to process optimisation. Comput Chem Eng 1999];23(December (10)):1509–26 http://linkinghub. elsevier.com/retrieve/pii/S0098135499003087. Banga J. Optimization in computational systems biology. BMC Syst Biol 2008];2(1):47. Beume N. Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem; 2006]. Bratley P, Fox BL. ALGORITHM 659: implementing Sobol’s quasirandom sequence generator. ACM Trans Math Softw 1988];14(March (1)):88–100 http://portal. acm.org/citation.cfm?doid=42288.214372. Bratley P, Fox BL, Niederreiter H. Implementation and tests of low-discrepancy sequences. ACM Trans Model Comput Simul 1992];2(July (3)):195–213 http:// portal.acm.org/citation.cfm?doid=146382.146385. Brockhoff D, Zitzler E. Are all objectives necessary? On dimensionality reduction in evolutionary multiobjective optimization. Parallel problem solving from nature – PPSN IX; 2006a]. p. 533–42. Brockhoff D, Zitzler E. Dimensionality reduction in multiobjective optimization with (partial) dominance structure preservation: generalized minimum objective subset problems. TIK Report 247; 2006]. Buxton A, Pistikopoulos EN. Environmental impact minimization through material substitution: a multi-objective optimization approach; 2004]. p. 407–17. Charnes A, Clower RW, Kortanek KO. Effective control through coherent decentralization with preemptive goals. Econometrica 1967];35(2):294–320 http:// www.jstor.org/stable/1909114. Charnes A, Cooper WW. Charnes Q, Cooper WW, editors. Management models and industrial applications of linear programming, vol. 1. John Wiley & Sons; 1961]. Charnes A, Cooper WW, Ferguson RO. Optimal estimation of executive compensation by linear programming. Manage Sci 1955];1(January (2)):138–51, http://pubsonline.informs.org/doi/abs/10.1287/mnsc.1.2.138. Copado-Méndez PJ, Guillén-Gosálbez G, Jiménez L. MILP-based decomposition algorithm for dimensionality reduction in multi-objective optimization: application to environmental and systems biology problems. Comput Chem Eng 2014, April. http://linkinghub.elsevier.com/retrieve/pii/S0098135414001112. Deb K, Saxena DK. On finding Pareto-optimal solutions through dimensionality reduction for certain large-dimensional multi-objective optimization problems. KanGal Report Number 2005011; 2005]. Du Y, Xie L, Liu J, Wang Y, Xu Y, Wang S. Multi-objective optimization of reverse osmosis networks by lexicographic optimization and augmented epsilon constraint method. Desalination 2014];333(January (1)):66–81 http://linkinghub. elsevier.com/retrieve/pii/S0011916413005109. Ehrgott M, Ryan DM. Constructing robust crew schedules with bicriteria optimization. J Multi-Criteria Decis Anal 2003];150(2002):139–50. Everson RM, Fieldsend JE, Singh S. Full elite sets for multi-objective optimisation (Acdm); 2002]. p. 1–12. Faure H, Lemieux C. Improved Halton sequences and discrepancy bounds. Monte Carlo Methods Appl 2010];16(January (3–4)):1–18 http://www.degruyter.com/ view/j/mcma.2010.16.issue-3-4/mcma.2010.008/mcma.2010.008.xml. Fox BL. ALGORITHM 647: implementation and relative efficiency of quasirandom sequence generators. ACM Trans Math Soft 1986];12(4):362–76. Gebreslassie BH, Waymire R, You F. Sustainable design and synthesis of algaebased biorefinery for simultaneous hydrocarbon biofuel production and carbon sequestration. AIChE J 2013];59(5):1599–621. Grossmann IE, Guillén-Gosálbez G. Scope for the application of mathematical programming techniques in the synthesis and planning of sustainable processes. Comput Chem Eng 2010];34(September (9)):1365–76 http://linkinghub. elsevier.com/retrieve/pii/S0098135409002841. Guillén-Gosálbez G. Optimal design and planning of sustainable chemical supply chains under uncertainty; 2008]. Guillén-Gosálbez G. A novel MILP-based objective reduction method for multiobjective optimization: application to environmental problems. Comput Chem Eng 2011]. Haimes YY, Lasdon LS, Wismer D. Integrated Optimization. IEE Trans Syst Man Cybern 1971];47(July):296–7. Halton JH. On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals. Numer Math 1960];2:84–90. Hischier R, Bauer C, Doka G, Dones R, Frischknecht R, Hellweg S, Jungbluth N, Loerincik Y, Margni M, Nemecek T. Implementation of life cycle impact assessment methods (3); 2010]. Ijiri Y. Management goals and accounting for control. Studies in mathematical and managerial economics. Amsterdam: North Holland; 1965]. Klee V. Can the measure of n 1 [ai, bi] be computed in less than O(n log n) steps? Am Math Mon 1977];84(4):284–5 http://www.jstor.org/stable/2318871. Kostin A, Guille G, Mele FD, Jime L. Identifying key life cycle assessment metrics in the multiobjective design of bioethanol supply chains using a rigorous mixedinteger linear programming approach; 2012]. ˇ cek L. Multi-objective optimisation for generating sustainKravanja Z, Cuˇ able solutions considering total effects on the environment. Appl Energy http://linkinghub.elsevier.com/retrieve/pii/ 2013];101(January):67–80 S030626191200311X. López Jaimes A, Coello Coello CA, Chakraborty D. Objective reduction using a feature selection technique. In: Proceedings of the 10th annual conference on Genetic and evolutionary computation. ACM; 2008]. p. 673–80.

48

P.J. Copado-Méndez et al. / Computers and Chemical Engineering 87 (2016) 36–48

Murillo-alvarado PE, El-halwagi MM. Optimization of pathways for biorefineries involving the selection of feedstocks, products, and processing steps; 2013]. Pozo C, Guillén G. Identifying the preferred subset of enzymatic profiles in nonlinear kinetic metabolic models via multiobjective global optimization and Pareto filters. PLOS ONE 2012];7(9):1–11. Purshouse R, Fleming P. Conflict, harmony, and independence: relationships in evolutionary multi-criterion optimisation. Evol Multi-Criterion Optim 2003]:67. Sabio N, Gadalla M, Guillén-Gosálbez G, Jiménez L. Strategic planning with risk control of hydrogen supply chains for vehicle use under uncertainty in operating costs: a case study of Spain. Int J Hydrogen Energy 2010];35(13):6836–52. Sabio N, Kostin A, Guillén-Gosálbez G, Jiménez L. Holistic minimization of the life cycle environmental impact of hydrogen infrastructures using multiobjective optimization and principal component analysis. Int J Hydrogen Energy 2012];37(March (6)):5385–405. Santiba nez Aguilar JE, González-Campos JB, Ponce-Ortega JM, Serna-González M, ElHalwagi MM. Optimal planning and site selection for distributed multiproduct biorefineries involving economic, environmental and social objectives. J Clean Prod 2014];65(February):270–94 http://linkinghub.elsevier.com/retrieve/pii/ S0959652613005234. Sobol I. On the distribution of points in a cube and the approximate evaluation of integrals. USSR Comput Math Math Phys 1967];7(January (4)):86–112 http:// www.sciencedirect.com/science/article/pii/0041555367901449.

Thoai NV. Criteria and dimension reduction of linear multiple criteria optimization problems. J Glob Optim 2011]:1–10. Uhlmann W. Introduction in the Monte Carlo methods. Biom Z 1965];7(2):116–9. Yue D, Slivinsky M, Sumpter J, You F. Sustainable design and operation of cellulosic bioelectricity supply chain networks with life cycle economic, environmental, and social optimization. Ind Eng Chem Res 2014];53(March (10)):4008–29, http://pubs.acs.org/doi/abs/10.1021/ie403882v. Zadeh L. Optimality and non-scalar-valued performance criteria. IEEE Trans Autom Control 1963];8(January (1)):59–60 http://ieeexplore.ieee.org/lpdocs/epic03/ wrapper.htm?arnumber=1105511. Zhang Q, Gong J, Skwarczek M, Yue D, You F. Sustainable process design and synthesis of hydrocarbon biorefinery through fast pyrolysis and hydroprocessing. AIChE J 2014];60(3). Zitzler E, Thiele L. Multiobjective optimization using evolutionary algorithms – a comparative case study. Lecture notes in computer science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics). Vol. 1498 LNCS; 1998]. p. 292–301 http://www.scopus.com/inward/record. url?eid=2-s2.0-84878545782&partnerID=tZOtx3y1. Zitzler E, Thiele L, Laumanns M, Fonseca C, da Fonseca V. Performance assessment of multiobjective optimizers: an analysis and review. IEEE Trans Evol Comput 2003];7(April (2)):117–32 http://www.scopus.com/inward/record.url?eid=2s2.0-0037936618&partnerID=tZOtx3y1.