Bi-objective integrating sustainable order allocation and sustainable supply chain network strategic design with stochastic demand using a novel robust hybrid multi-objective metaheuristic

Author's Accepted Manuscript

Bi-objective integrating sustainable order allocation and sustainable supply chain network strategic design with stochastic demand using a novel robust hybrid multi-objective metaheuristic Kannan Govindan, Ahmad Jafarian, Vahid Nourbakhsh www.elsevier.com/locate/caor

PII: DOI: Reference:

S0305-0548(15)00009-X http://dx.doi.org/10.1016/j.cor.2014.12.014 CAOR3713

To appear in:

Computers & Operations Research

Cite this article as: Kannan Govindan, Ahmad Jafarian, Vahid Nourbakhsh, Biobjective integrating sustainable order allocation and sustainable supply chain network strategic design with stochastic demand using a novel robust hybrid multi-objective metaheuristic, Computers & Operations Research, http://dx. doi.org/10.1016/j.cor.2014.12.014 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


                                    Kannan Govindana1, Ahmad Jafarianb, Vahid NourbakhshC a

b

Department of Business and Economics, University of Southern Denmark, Odense, Denmark. Faculty of Management and Accounting, Allameh Tabataba’i University Business School, Tehran, Iran. c The Paul Merage School of Business, Irvine, CA 92617, United States Corresponding author, Email address: [email protected]



 Sustainability has been considered as a growing concern in supply chain network design (SCND) and in the order allocation problem (OAP). Accordingly, there still exists a gap in the quantitative modeling of sustainable SCND that consists of OAP. In this article, we cover this gap through simultaneously considering the sustainable OAP in the sustainable SCND as a strategic decision. The proposed supply chain network is composed of five echelons including suppliers classified in different classes, plants, distribution centers that dispatch products via two different ways, direct shipment, and cross-docks, to satisfy stochastic demand received from a set of retailers. The problem has been mathematically formulated as a multi-objective optimization model that aims at minimizing the total costs and environmental effect of integrating SCND and SMP, simultaneously. To tackle the addressed problem, a novel multi-objective hybrid approach called MOHEV with two strategies for its best particle selection procedure (BPSP), minimum distance, and crowding distance is proposed. MOHEV is constructed through hybridization of two multi-objective algorithms, namely the adapted multi-objective electromagnetism mechanism algorithm (AMOEMA) and adapted multi-objective variable neighborhood search (AMOVNS). According to achieved results, MOHEV achieves better solutions compared with the others, and also crowding distance method for BPSP outperforms minimum distance. Finally, a case study for an automobile industry is used to demonstrate the applicability of the approach. Key words: Supply chain network design, Multi-objective optimization, Sustainability, Stochastic demand, Robust metaheuristic, Greenhouse gas emissions.





























































1

Corresponding author, Email address: [email protected]

1. Introduction Nowadays, global concerns about sustainability have increased substantially. Activists, media, consumers, and non-governmental organizations (NGOs) are stimulating the corporations to extend their responsibility and consider sustainability in their decisions and operations. According to Brundtland (1987), sustainable development is to meet the needs of the present without compromising the ability of future generations to meet their own needs. In turn, to operationalize sustainability, Elkington (1998) defined the three pillars of sustainability as economy, environment, and society. In this way, business sustainability is defined as the ability to do business with the aim of economic, environmental, and social well-being (Hassini et al., 2012). Regarding the environmental aspect of sustainability, governmental legislations and policies — for instance, in Europe, Japan, and North America — are obliging the corporations to minimize their environmental impacts (Fleischmann et al., 2002; Robeson et al., 1992). As our focus in this paper is on supply chain management (SCM), we use the following definition of sustainable supply chain management (SSCM) taken from Seuring (2013): “SSCM is the management of material, information and capital flows as well as cooperation among companies along the supply chain while taking goals from all three dimensions of sustainable development, i.e., economic, environmental and social, into account which are derived from customer and stakeholder requirements.” Supply chain network design (SCND) and supplier management are considered as the most crucially important strategic decisionsin SCM that play a role in overall sustainable performance of the supply chain (Pishvaee and Razmi, 2012; Seuring and Müller, 2008; Shen et al., 2012). SCND contains the determination of locations, numbers, and capacities of network facilities and the material flow between them (Pishvaee and Razmi, 2012). Sustainable SCND (SSCND) tries to define the best supply chain configuration that enables an organization to maximize its longterm economic profitability as well as environmental and social performance (Chaabane et al., 2012). Also, supplier management is a source of increasing concern in connention with training managers in sustainable management, legislation, expanding of organizations responsibility, and paving the way to support sustainable managerial decision making (Dao et al., 2011). In an SSCM perspective, the supplier management decision is one of the critical issues that help organizations to maintain a strategically competitive position through operations and purchasing management. In the supplier management process, the potential suppliers are reviewed, evaluated, and selected to become part of the company’s supply chain, and subsequently the optimum order quantity of each supplier is determined (Weber et al., 1991). Sustainable supplier management is clearly a critical activity in purchasing management due to the environmental, social sustainability, and also ecological performance that can be demonstrated by suppliers (Godfrey, 1998). In other words, companies must cooperate with their suppliers about socially and environmentally friendly practices for purchasing and materials management(Hsu and Hu, 2009). Accordingly, supply chains are shifting toward SSCM with various motivations, such as gaining public image (Fombrun, 2005), satisfying activists’ requirements (Spar and La Mure,

2003), and maintaining customers over the long term (Bhattacharya and Sen, 2004). In this regard, there is a gap in the previously published papers in the area of sustainable network design (Chaabane et al., 2012; Elhedhli and Merrick, 2012; Hu and Li, 2009; Pishvaee and Razmi, 2012; Shen and Daskin, 2005; Srivastava, 2007; Vachon and Klassen, 2008; Wang et al., 2011). These papers lack the simultaneous consideration of network design and the order allocation problem (OAP), which can be attributed to the complexity of modeling economical and, to some extent, sustainable aspects. In other words, many studies discuss SSCN and sustainable supplier selection but rarely develop supplier management methods to solve the integration of SSCND problems and multi-sourcing sustainable order allocation problems in SSCM. Therefore, the main contribution of this study includes developing a SSCND model integrated with sustainable order allocation by considering both traditional aspects, i.e. economical, and environmental aspects. Based on the aforementioned considerations, this paper addresses the issue of multiobjective, multi-echelon supply network design including suppliers, plants, distribution centers, and retailers. The major innovations and features that distinguish this study from the abovementioned ones are: (i) the economical and environmental impacts of suppliers and opening facilities are considered; (ii) the decision maker is allowed to embed his/her priority against different suppliers in the model in terms of a grading scheme; (iii) the inherent uncertainty of customers’ demands is considered; (iv) two types of shipment are allowed in the model: direct shipment from distribution centers to retailers and indirect shipment from distribution centers to cross-docks to retailers; and (v) it should be noted that logistics network design is considered as an NP-hard problem (Jo et al., 2007; Syarif et al., 2002). Thus, we have developed and compared three novel metaheuristics to find Pareto-optimal solutions. The rest of the paper has the following structure. Section 2 briefly delineates the literature review. In Section 3, the problem is described in detail. The mathematical formulation is explained in Section 4. Section 5 describes the definitions of Pareto optimality and the structure of the algorithms. Numerical problems developed to study the performance of the proposed algorithms, parameter settings, and comparisons among algorithms are presented in Section 6. The proposed model is implemented for the case study in Section 7, and finally Section 8 is devoted to conclusions and future directions for research. 2. Literature review In recent years, sustainability has fascinated practitioners and scholars based on simultaneous pressure from various stakeholders including consumers, managements, nongovernmental organizations (NGOs), governmental legislation, community activists, and global competition (Govindan et al., 2013; Hassini et al., 2012). Along the way, Tang and Zhou (2012) indicated that at the mercy of recent considerations to sustainability, the operations research and management science development lies within the domain of environmentally and socially sustainable operations. Also, Ageron et al. (2011) have pointed out that the implementation of sustainability in supply networks is considered as a critical factor for the success of the whole

supply chain. They indicated that sustainable responsibility can lead to companies’ development and competitive advantage, which is not possible without incorporating SSCM practices. SSCM is defined as the management of information and material flows and cooperation among organizations along the supply chain with triple bottom line (TBL) selection criteria that contain all three dimensions of sustainable development, i.e. economic, environmental, and social dimensions. SSCM aims to minimize environmental concerns and improve social impacts while, at the same time, improving long-term economic performance(Erol et al., 2011; Seuring and Müller, 2008). TBL suggests that in addition to economic performance, organizations need to engage in activities that positively affect the environment and society (Elkington, 1997). Among different publications about SSCM, our paper is about sustainable SCND which aims at considering the above-mentioned three pillars of sustainability in the design phase (Frota Neto et al., 2008). SCND is a new, emerging approach that seeks to embed economic, environmental, and social decisions into supply chains during the design process (Chaabane et al., 2012). Despite the fact that the field of SSCM is considered as being quite new, the interest in SSCM has increased rapidly over recent years (Ageron et al., 2011). Seuring and Müller (2008) comprehensively reviewed a total of 191 papers about SSCM published from 1994 to 2007 and outlined major streams of research in this field. Successively, Seuring (2013) has reviewed in detail 36 papers which have utilized quantitative models for SSCM among more than 500 papers published until 2010. A large number of papers on sustainability and SSCM have recently been published in different journals (e.g. Bai and Sarkis, 2010; Basurko and Mesbahi, 2012; Baumgartner, 2011; Bergenwall et al., 2012; Gimenez et al., 2012; Gunasekaran and Spalanzani, 2011; Linton et al., 2007; Pagell et al., 2008; Seuring and Gold, 2013; Smith and Ball, 2012; Tseng et al., 2013; 



). For instance, Carter and Rogers (2008) concentrated their research on preparing a literature review with a conceptual framework for SSCM through reviewing journal papers and related publications from previous years. Teuteberg and Wittstruck (2010) published a systematic review of SSCM. The published papers about SSCM mainly focus on greenness and cost-effectiveness of networks. Only four out of 36 surveyed papers (Cruz and Matsypura, 2009; Cruz, 2008; Cruz, 2009; Hsueh and Chang, 2008) have stated that social aspects as well, as economic and environmental aspects are considered. However, these papers have misused the term corporate social responsibility (CSR), while their models have only encompassed economic and environmental issues. Altogether, social impacts are inherently qualitative, which causes the rareness of quantitative models in this dimension of SSCM. Consequently, we are obliged to narrow down our literature review to green supply chain management (GSCM). For a more detailed and complete review of SSCM, readers are referred to Srivastava (2007), Carter and Rogers (2008), Seuring and Müller (2008), Ageron et al. (2011) and Brandenburg et al., 2014.

According to reviews conducted by Srivastava (2007) and Zhu and Sarkis (2004), GSCM includes green selection of raw material suppliers, product shipment through the chain from suppliers to plants to distributors to customers, and even reverse logistics. The neighboring field of our paper is the reverse SCND problem, which deals with decisions about the number,location, and capacity of collection, recovery, recycling, and disposal centers along with decisions about material flows between them (Fleischmann et al., 2004). Fleischmann et al. (1997) provided the first review on quantitative models for reverse logistics. For a more detailed and complete review of reverse and closed loop supply chain management, readers are referred to Govindan et al., 2014 and Sasikumar and Kannan (2008a-b, 2009). Reverse SCND has attracted a lot of researchers so far (Du and Evans, 2008; Ko and Evans, 2007; Lieckens and Vandaele, 2007; Pishvaee et al., 2010a; Salema et al., 2007; Srivastava, 2008; Soleimani and Govindan, 2014; Govindan and Popiuc, 2014; Kannan et al., 2012). Closed-loop SCND is another topic that is closely related to our study (Chaabane et al., 2012; Pishvaee et al., 2010b; Qiang et al., 2013; Wang and Hsu, 2010; Kannan et al., 2009). Among various proposed quantitative models of closed-loop SCND, the most relevant ones are those proposed by Pishvaee and Razmi (2012) and Pinto-Varela et al. (2011), as they both consider a trade-off between environmental impacts and economic performance of network configurations.As mentioned above, there is a large body of literature about reverse and closed-loop SCND; however, according to Chaabane et al. (2012), collecting used products and re-processing might not only increase operating costs but also contribute to an increase in greenhouse gas (GHG) emissions, which defeats long-term sustainability. Accordingly, our proposed model is focused on designing a forward GSCND. Among the plethora of studies of forward GSCND, in terms of supply chain configuration, the most relevant one is a SCND proposed by Bachlaus et al. (2008), as it consideres five echelons and utilizes the cross-dock concept as an intermediate level between distribution centers and end customers. In terms of modeling approach, objectives and the trade-off between objectives, the most relevant work is a mixed integer linear programming model proposed by Wang et al. (2011) which has two objectives: (i) minimization of total costs (i.e., fixed setup costs, environment protection investment, transportation cost, and handling cost) and (ii) minimization of total environmental impacts(i.e., CO2 emission produced by product shipment and manufacturing).In this paper, different levels of investments at plants with the aim of environmental protection are considered.


3. Problem description In this study, a mixed integer programming model is proposed for a single product forward supply chain with five echelons: suppliers, plants, distribution centers, cross-docks, and retailers (Fig.1). [Please insert Figure about here]

One of the major contributions in/features of the proposed model is that suppliers have to be selected from a pool of graded potential suppliers. In other words, it is assumed that suppliers are classified into some limited categories called grades based on their different attributes (e.g., quality, delivery, environmental impacts, and performance history). Different methods can be used for ranking suppliers, such as analytic hierarchy process (AHP) firstly proposed by Nydick and Hill (1992)for vendor selection, analytical network process (ANP) proposed by Sarkis and Talluri (2002) for supplier selection, and data envelopment analysis (DEA) firstly proposed byWeber (1996) for vendor selection. In this study, we assume that the grade of each supplier is given. In this paper, the environmental impact of all the members of a supply chain, including the most harmful GHG emissions like CO2, CFC, and NOx, is taken into account. In addition, environmental impacts related to open/established plants, distribution centers, and cross-docks are taken into account. In the proposed model, products are shipped from suppliers to plants to distribution centers to cross-dock to retailers. However, products can be shipped either directly from distribution centers to retailers or via cross-docks, and the model tries to find the amount to be shipped directly and indirectly. Last but not least, the inherent uncertainty of retailers’ demands is considered by introducing a scenario-based stochastic model. The proposed SCND model is bi-objective: minimization of both costs and environmental impacts. In the cost objective,seven types of costs are taken into account: (i) developing long-term partnerships with the supplier(s), (ii) opening the facilities (i.e., plants, distribution centers, and cross-docks), (iii) technology establishment at plants. In supply chain strategic design problems, different production technologies/production lines are available to be established at plants as a strategic decision which affects both costs and environmental impacts (Pishvaee et al., 2012; Shen, 2007), (iv) purchasing raw materials from suppliers, (v) manufacturing products at plants, (vi) handling products at cross-docks and distribution centers, and (vii) shipment of products through the network. It will be discussed later in detail that there may be a trade-off between the above mentioned costs and the following environmental impacts in the environment objective: (i) environmental impacts of shipping products through the network, (ii) environmental impacts caused by opening facilities (i.e., plants, distribution centers, and cross-docks), and (iii) manufacturing environmental impacts at plants related to technology. In the proposed SCND model, the goal is to find the supply chain configuration that minimizes both costs and environmental impacts. The decisions to be made are as follows: (i) choosing suppliers from a set of potential suppliers that are categorized into some grades by the policy makers according to their attributes,

(ii) locating facilities (i.e., plants, distribution centers, and cross-docks), (iii) choosing the production technology of plants, (iv) determining the amount of orders from chosen suppliers, and (v) determining the flow of products through the network from suppliers to retailers. The assumption, notations, and proposed model are presented in Section 4. 4. Mathematical formulation The proposed model is based on the following common assumptions in the literature (Syarif et al., 2002; Yao and Hsu, 2009): • The demand of customers must be satisfied • A single type of product is considered • Facilities cannot order more than received orders (storage is not allowed) • Back order is not permitted in facilities • The number of facilities in each echelon as well as their potential sites is restrained by predefined values • There is no flow between the facilities of the same echelon • Suppliers, plants, distribution centers, and cross-docks have limited capacity • Distribution centers and cross-docks must satisfy retailers’ demand • Manufacturers cannot send production directly to retailers or cross-docks. The notations used in the model are presented in Tables 1 to 3. [Please insert Table 1, 2 and 3about here]


In this model, the uncertainty of retailers’ demands is taken into account by formulating a scenario-based two-stage stochastic programming (TSP) model. In the first stage, decisions (i)(iii) about choosing the supply base, location of facilities, and technology of plants are made prior to knowing the demands. In the second stage, decisions (iv)-(v) about the order quantities from suppliers and the flow of products are made after the demands are known. The uncertainty of demands is captured by postulating a number of scenarios, {1, 2, ..., |ξ |} , each with a possibility, P(ξ ) , that probabilities satisfy the condition of

¦ξ

∈Ξ

P (ξ ) = 1 . We begin the model

by introducing the cost objective: Min OB1 = ¦ ¦ FC s g α s g + ¦ ¦ (OC p + TEC pt ) β pt + ¦ OC d γ d + ¦ OC c λc g ∈G s g ∈S g

p ∈P t ∈T

+ ¦ P (ξ ){ ¦ ξ ∈Ξ

¦ ¦ (V C

d ∈d sg

+ TC s g p ) x s g p (ξ ) + ¦ ¦ MC pt h pt (ξ )

g ∈G s g ∈S g p ∈P

p ∈P t ∈T

+ ¦ ¦ (V C d + TC pd ) x pd (ξ ) + p ∈P d ∈D

+ ¦ ¦TC cr x cr (ξ ) + c ∈C r ∈R

c ∈C

¦ ¦ (V C

c

+ TC dc ) x dc (ξ )

(1)

d ∈D c ∈C

¦ ¦TC

dr

x dr (ξ )}

d ∈D r ∈R

The cost objective minimizes total variable and fixed costs of the supply chain network. The first summation is the fixed cost of establishing a business with suppliers. The second

summation is the fixed cost associated with operating/opening plants and establishing/acquiring technologies. The subsequent two terms are the fixed opening costs of distribution centers and cross-docks, respectively. While the above-mentioned terms are associated with the first stage of the TSP model decided before knowing the demands, the remaining terms are associated with the second stage determined after the scenarios of demands have become known. The fifth summation of the objective function in this stage is the variable purchasing cost of raw materials and the variance transportation cost from suppliers of plants. The next two summations represent the variable manufacturing cost at plants, and the variable handling cost at distribution centers and the variable transportation costs from plants to distribution centers, respectively. The eighth term is the variable handling cost at cross-docks and the variable transportation cost from distribution centers to cross-docks. The subsequent summation represents the transportation cost from crossdocks to retailers. In this model, direct shipments of products from distribution centers to retailers without using cross-docks are allowed. To explain, products could be shipped from distribution centers to retailers either through cross-docks to decrease the costs of less than a truck load or directly from distribution centers to retailers when a large amount of products is supposed to be shipped from a specific distribution center to a specific retailer. This cost of direct shipment from distribution centers to retailers is represented in the final term. The environment objective is as follows: Min OB2 =

¦ ¦ ES

sg

α s + ¦ ¦ EO pt β pt + ¦ EOd γ d + ¦ EOc λc g

g ∈G s g ∈S g

p ∈P t ∈T

d ∈D

+ ¦ P (ξ ){ ¦ ¦ EM pt h pt (ξ ) + ¦ ξ ∈Ξ

p ∈P t ∈T

c ∈C

¦ ¦ ET

sg p

g ∈G s g ∈S g p ∈P

x s g p (ξ ) + ¦ ¦ ET pd x pd (ξ ) p ∈P d ∈D

(2)

+ ¦ ¦ ET dc x dc (ξ )+ ¦ ¦ ET cr x cr (ξ ) + ¦ ¦ ET dr x dr (ξ )} d ∈D c ∈C

c ∈C r ∈R

d ∈D r ∈R

The first term of this objective represents the environmental impact related to suppliers which are selected as the first echelon/level of supply chain. The second to fourth terms are the environmental impacts associated with opening facilities/plants, distribution centers, and crossdocks, respectively. The fifth summation in the second stage is the environmental impacts produced by manufacturing at plants. Finally, the last five summations are, respectively, environmental impacts related to shipping raw materials from suppliers to plants and products from plants to distribution centers, cross-docks and finally retailers or directly from distribution centers to retailers. The constraints of the model are as follows:

¦¦x

sg p

(ξ ) = ¦ x pd (ξ ) ∀p ∈ P , ∀ξ ∈Ξ

g∈G sg ∈S g

¦x

pd

p∈P

(ξ ) = ¦ xdc (ξ ) + ¦ xdr (ξ ) ∀d ∈ D , ∀ξ ∈Ξ c∈C

¦x

dc

d ∈D

d ∈D

r∈R

(ξ ) = ¦ xcr (ξ ) ∀c ∈ C , ∀ξ ∈Ξ r∈R

(3) (4) (5)

¦h

pt

(ξ ) = ¦

t∈T

¦x

sg p

(ξ ) ∀p ∈ P , ∀ξ ∈Ξ

(6)

g∈G sg ∈S g

¦x

(ξ ) + ¦ xcr (ξ ) ≥ Demr (ξ ) ∀r ∈ R , ∀ξ ∈Ξ

dr

d ∈D

(7)

c∈C

¦y

cr

(ξ ) ≤ 1 ∀r ∈ R , ∀ξ ∈Ξ

(8)

c∈C

xcr (ξ ) ≤ Demr (ξ ) ycr (ξ ) ∀c ∈ C , ∀r ∈ R , ∀ξ ∈Ξ

¦x

sg p

(9)

(ξ ) ≤ Capsg α sg ∀g ∈G , ∀sg ∈ Sg , ∀ξ ∈Ξ

(10)

p∈P

¦¦x

sg p

(ξ ) ≤ ¦ Cap pt β pt ∀p ∈ P , ∀ξ ∈Ξ

g∈G sg ∈S g

¦x

(11)

t∈T

(ξ ) ≤ Capd γ d ∀d ∈ D , ∀ξ ∈Ξ

(12)

(ξ ) ≤ Capc λc ∀c ∈ C , ∀ξ ∈Ξ

(13)

pd

p∈P

¦x

cr

r∈R

¦β

pt

≤ 1 ∀p ∈ P

(14)

t∈T

¦ ¦α

sg

≤ S Max

(15)

g∈G sg ∈S g

¦¦ β

pt

≤ PMax

(16)

p∈P t∈T

¦γ

d

≤ DMax

(17)

c

≤ CMax

(18)

d∈D

¦λ c∈C

¦

αs ≥ g

¦ ¦x

g ′∈{ g -1,g -2...,1} s g ∈S g p ∈P

sg p

(ξ ) −

¦

¦ Cap

g ′∈{ g -1,g -2...,1} s g ∈S g

sg

+ 1 ∀g ∈{2,..., G } , (19)

∀sg ∈ Sg

αs ≥ 0 ∀g ∈{2,..., G } , ∀sg ∈ Sg

(20)

g

αs ≥ g

¦

¦ ¦x

g ′∈{ g -1,g -2...,1} s g ∈S g p ∈P

sg p

(ξ ) −

¦

¦ Cap

g ′∈{ g -1,g -2...,1} s g ∈S g

sg

+ 1 − M δsg (21)

∀g ∈{2,..., G } , ∀sg ∈ Sg αs ≤ g

¦

¦ ¦x

g ′∈{ g -1,g -2...,1} s g ∈S g p ∈P

sg p

(ξ ) −

¦

¦ Cap

g ′∈{ g -1,g -2...,1} s g ∈S g

sg

+ 1 + M δs g (22)

∀g ∈{2,..., G } , ∀sg ∈ Sg

α s ≥ −M (1 − δ s ) ∀g ∈{2,..., G } , ∀sg ∈ Sg

(23)

α s ≤ M (1 − δ s ) ∀g ∈{2,..., G } , ∀sg ∈ Sg g

(24)

α s ∈{0,1} ∀g ∈G , ∀sg ∈ Sg

(25)

g

g

g

g

β pt ∈{0,1} ∀p ∈ P , ∀t ∈ T

(26)

γ d ∈{0,1} ∀d ∈ D

(27)

λc ∈{0,1} ∀c ∈ C

(28)

ycr (ξ ) ∈{0,1} ∀c ∈ C , ∀r ∈ R , ∀ξ ∈Ξ

(29)

xsg p (ξ ) ≥ 0 ∀g ∈G , ∀sg ∈ Sg , ∀p ∈ P , ∀ξ ∈Ξ

(30)

x pd (ξ ) ≥ 0 ∀p ∈ P , ∀d ∈ D , ∀ξ ∈Ξ

(31)

xdc (ξ ) ≥ 0 ∀d ∈ D , ∀c ∈ C , ∀ξ ∈Ξ

(32)

xcr (ξ ) ≥ 0 ∀c ∈ C , ∀r ∈ R , ∀ξ ∈Ξ

(33)

xdr (ξ ) ≥ 0 ∀d ∈ D , ∀r ∈ R , ∀ξ ∈Ξ

(34)

hpt (ξ ) ≥ 0 ∀p ∈ P , ∀t ∈ T , ∀ξ ∈Ξ

(35)

Constraints (3) to (5) guarantee the flow conservation. Constraints (6) should be satisfied to compute the amount of products which have to be handled/manufactured in plant p under technology t and scenario ξ , hpt (ξ ) . Constraints (7) ensure that the demand of each retailer is satisfied by a cross-dock and/or distribution center(s). Constraints (8) state that each retailer should be assigned at most to one cross-dock in order to minimize the cost of less than truck load.Constraints (9) ensure that products are allowed to be transported from cross-dock c to retailer r under scenario ξ if and only if the retailer r is assigned to cross-dock c. Constraints (10) imply that products could be shipped only from selected suppliers. In the same way, constraints (11) to (13) ensure that products could be shipped only from opened/operateing plants, distribution centers, and cross-docks, respectively. Constraints (14) ensure that in each potential node for plants at most one plant with only one technology is established. Constraints (15) to (18) limit the maximum number of selected suppliers and opened facilities (i.e., plants, distribution centers, and cross-docks, respectively) according to the decision maker’s attitude. Constraints (19) to (24) are the linear equivalents of constraints (36). We introduced additional binary decision variables, δ sg , and a suitably large number, M, to transform constraints (34) to the linear equivalents.

°­

α s ≤ max ® g

¦

¦ ¦x

°¯ g ′∈{ g -1,g -2...,1} s g ∈S g

sg p

(ξ ) −

p ∈P

¦

¦ Cap

g ′∈{ g -1,g -2...,1} s g ∈S g

sg

°½ + 1, 0 ¾ °¿

(36)

∀g ∈{1, 2,..., G −1} , ∀sg ∈ Sg , ∀ξ ∈Ξ Constraints (36) ensure that a supplier of grade g (i.e. a supplier which belongs to Sg ) can be selected if and only if the whole capacity of the suppliers of higher grades (i.e., suppliers belonging to {S g −1
S g −2...
S1} ) are consumed. Finally, constraints (25) to (35) define the variables' types.




5. Multi-objective optimization problem 5.1. Mathematical formulation A multi-objective optimization problem with k objective functions can be expressed mathematically as:

     Minimize/maximize f (x ) = [f 1 (x ), f 2 (x ), ..., f k (x )]T Subject to:

 g j (x ) ≤ 0

 g k (x ) ≤ 0  x ∈R n where

(37)

j = 1, 2,…, n

k = 1,2,…, k

 x = [x 1 , x 2 ,..., x N ]T

is

a

vector

consisting

of

N decision

variables;

each

x n ( n = 1,2, . . . , N ) can be real-valued, integer-valued or Boolean-valued. Objective functions,

f i : R n → R (i = 1, 2,. . ., k ) , and constraints, g j ( j = 1,2, . . . , n ) and g k ( k = 1,2, . . . , m ) , can be linear or nonlinear arbitrary functions. The multi-objective optimization should solve the    problem and find x ∗ where x ∗ = [x 1∗ , x 2∗ ,..., x N ∗ ]T . x ∗ must satisfy all the constraints and  minimize/maximize f (Tsou and Kao, 2008). In this regard, a variety of approaches can be used to solve this problem. One popular approach is to combine those objectives into a single composite objective so that traditional mathematical programming methods can be applied. To this end, a value or utility function needs to be identified according to the preference of decisionmaker(s). The simplest method is to assume independent preferences among those objectives and apply an additive utility function. On the other hand, instead of transforming the original problem into a single-objective one, the Pareto optimal concept based on non-dominance can be utilized. Pareto dominance and non-dominance can be determined through multiple pair wise



vector comparisons (Moslemi and Zandieh, 2011) . More specifically, let y = ( y 1, y 2 ,…, y N

) be

another vector containing N decision variables. In a maximization problem, it is said that solution x dominates solution y if and only if:

    f k ( x ) ≥ f k ( y ) ∀k , and f k (x ) > f k ( y ) for at least one k ∈{1, 2,…, K } .
(38) Let A be a set containing all thenon-dominated solutions in X. Then, the set A is called

the Pareto optimal set or the Pareto frontier of the multi-objective optimization problem. By introducing the Pareto optimality concept, more choices may be provided to decision-makers with different perspectives. In many cases, however, the number of solutions in the Pareto optimal solution set increases as the number of conflicting objectives increases, which may not be desirable to a decision-maker. The novel method proposed in this paper is expected to fill the

gap between the single solution and the Pareto optimal set by providing decision-makers with a medium-sized set of solutions (several representative solutions) from a holistic view (Li et al., 2009).

5.2. Multi-objective algorithms Regarding solution methods, we have proposed a hybrid multi-objective metaheuristic algorithm called MOHEV, which is a combination of the adapted multi-objective electro magnetism mechanism algorithm (AMOEMA) and the adapted multi-objective variable neighborhood search (AMOVNS). Furthermore, to examine the proposed method, some experiments are designed and conducted, in which the results of the hybrid algorithm are compared to the Nondominated Sorting Genetic Algorithm (NSGAII) as one of the best multi-objective evolutionary genetic algorithms based on methods (Deb et al., 2002) and multi-objective particle swarm optimization algorithms (MOPSO) as one popular swarm intelligent algorithm (Coello Coello and Lechuga, 2002). These metaheuristic algorithms are utilized to obtain the best design of SSCND that minimizes both costs and environmental effects. Below, the MOHEV method is explained in detail. 5.2.1. MOHEV The construction of hybrid metaheuristics is motivated by the need to achieve a good trade-off between the global exploration and the local exploitation during the search (Behnamian et al., 2009). Furthermore, the combination of certain aspects of different metaheuristics as a new, growing field of research becomes more and more important and feasible due to increasing computational power. Accordingly, the combination of AMOVNS as an algorithm that gives the search process a time-varying and ultimate zero probability of jumping (Zhang et al., 2007) and AMOEMA which utilizes the electromagnetism theory of physics by considering each particle to be an electrical charge (Tsou, 2009) paves the ground to avoid falling into the local minimum and eventually find the global optimum. In this regard, the MOHEV method integrates drastic features from AMOEMA and AMOVNS with the aim of benefiting from both algorithms’ strengths while covering their weak points. We expect that hybridization of such two different algorithms enriches the search process and enhances the overall or local search capability and efficiency (Behnamian and Ghomi, 2010). I) Main loop (AMOEMA) Electromagnetism mechanism algorithm (EMA) is a metaheuristic method introduced by Birbil and Fang (2003) for global optimization. It is a population-based metaheuristic structure to search for the optimal solution in continuous optimization problems (Tsou and Kao, 2006) and has rarely been used for discrete optimization problems (Vahdani and Zandieh, 2010). EMA utilizes the attraction–repulsion mechanism of electromagnetism theory, which is based on Coulomb’s law to determine the optimal solution (Su and Lin, 2011). EMA and evolutionary algorithms are similar to each other in terms of structural aspects; mainly, they are population-

based approaches with information sharing among members. It seems a natural progression to extend EMA to multi-objective optimization problems (MOOP) (Tsou et al., 2008). Tsou and Kao (2006) restructured EMA to solve a MOOP and called it MOEMA. There are four phases in MOEMA: “Initialization”, “calculation of resultant force”, “movement along the direction of the force” and “local search”. The pseudo-code of the MOHEV main loop is presented in Fig. 2.


[Please insert Figure 2 about here] II ) Encoding scheme and initialization As MOEMA is considered as a continuous optimization, the random-key (RK) technique as an encoding scheme is utilized with the aim of enabling the approach to solve discrete problems(Chang et al., 2009; Naderi et al., 2010; Tavakkoli-Moghaddam et al., 2009). RK is easily applied to parse the primary solution. The encoding of a sub-solution in the plants echelon, the most important and complicated echelon, is schematically illustrated in Fig. 3. As shown in Fig. 3, firstly, a matrix with P elements, each with uniform distribution U(0,1) is generated where RK demonstrates the desired number of plants restricted with PMax . For example, the encoded solution {0.86, 0.08, 0.23, 0.65, 0.52} represents the parsed solution {1,0,0,1,1} with Pmax = 3 . Then, according to the solution of the first step, a P2 × D

while

P2 ≤ PMax random matrix is generated as the next part of the solution. Accordingly,the flow of products from selected plants to the distribution centers is determined. In this step, we can intelligently reduce the dimensions of the second matrix by selecting the opened plants and opened distribution centers. Hence, we can save plenty of time. Subsequently, the columns of the second matrix are normalized to determinehow distribution centers assign their demandsto selected plants. For instance, the encoded solution for the first distribution center {0.49, 0.95, 0.34} represents the parsed solution {0.27, 0.54, 0.19} . In the next step, to specify a technology level for each plant, we produce a uniform distributed random number, for each plant. Then, RK shows the chosen technology level by rounding these numbersfor instance, the encoded solution {3.26, 2.71,1.16} represents the parsed solution {3,3,1} which indicates that the third technology should be installed in the first and fourth plants and the first technology in the fifth plant. [Please insert Figure 3 about here]


The initialization procedure is used to sample PopSize particles which constitute the population. These particles are randomly chosen from the feasible region, which is an N dimensional hyperspace. The initial value is assumed to be uniformly distributed between the corresponding upper and lower bounds (Su and Lin, 2011). Upon the particle generation, a conversion method is used to parse the primary solutions with the corresponding values of each particle in the discrete area as illustrated in Fig. 4.

[Please insert Figure 4 about here] III) Calculating the resultant force In MOEMA, each solution is considered as a charged particle. Each particle’s charge is calculated based on goodness of the corresponding solution that is similar to the fitness function concept in genetic algorithm (Khalili and Tavakkoli-Moghaddam, 2012). The electrostatic force between the charges of two particles is proportionally calculated based on magnitudes of the particles’ charge, and it is inversely proportional to the distance between the charges (Vahdani and Zandieh, 2010). For a minimization problem, the charge of each particle is estimated as follows (Tsou and Kao, 2006):   § · prox( x i ) − prox( x best ) ¨ ¸ , i = 1, 2,..., PopSize q = exp − N (39) j  best ¸ PopSize ¨ (prox( x ) − prox( x )) ¦ j =1 © ¹ where A is the current non-dominated set obtained from the best design of each iteration, N is the i

 dimension of the solution space, and x best is the particle from non-dominated front A selected by two strategies, nearest/minimum distance (MD) and S1,and crowding distance (CD), S2. The first strategy is as follows:  i   best  prox( x i ) = min ( f ( x ) − f (x ) ) best  x

(40)

∈A

and  i x best = arg min  i {prox( x )} ∀i

(41)

x

 In this strategy, the charge of each point x j is calculated by its minimum distance to A

 

 

calculated by its objective vector f (x j ) and non-dominated solution f (x best ) , for each   x best ∈ A . In the second strategy, x best is selected based on CD. The CD procedure is illustrated in Fig. 5. The procedures of best particle selection by strategies S1 and S2 are demonstrated in Fig. 6.
[Please insert Figure 5 about here] [Please insert Figure 6 about here] The objective value of the present best point and the charge of each point would be changed at each iteration.Obviously, particles with lower objective value have higher charge (Su and Lin, 2011). In addition, agood solution encourages the otherparticles to converge in an attractive point. In contrast, a bad solution obliges the other particles to move toward this region (Vahdani and Zandieh, 2010). By comparing the objective-function values of two particles, the  direction of a particular force between them is concluded upon. The resultant force vector, F j ,  which exerts force on the corresponding particle, x j , is determined by the following equation:

­  j i qiq j ( x x ) − °   x i −x j ° ° m i  qiq j °  F = ¦ ®( x i − x j )   j =1, j ≠ i ° x i −x j ° °do nothing, ° ¯

2

,

2

,

  ½ if prox(x j ) < prox( x i ) ° ° ° j i ° if prox( x ) > prox(x ) ¾, i = 1, 2,..., PopSize ° ° otherwise ° ° ¿

(42)

  If particle i dominates particle j (x i  x j ) , particle jis attracted by particle i , and vice

versa. However, if two particles are not dominating each other, they stay neutral in the population.Technically speaking, MOEMA does not seek to search for optimal solutions but tries to find the efficient solutions that can be expressed in terms of non-dominated points in the objective space.  The forces exerted on x i by each of the other points are calculated by vector summation   as illustrated in Fig. 7 (Tsou and Kao, 2006). F 13 is aforce exerted by particle x 1 on particle       x 3 (the repulsion phenomenon indicates that f (x 1 ) is worse than f (x 3 ) ) and F 23 is the force     exerted byparticle x 2 on particle x 3 (the attraction phenomenonimplies that f (x 2 ) is better than       f (x 3 ) ). The resultant force exerted on x 3 equals F 3 = F 13 + F 23 . The resultant force calculation procedure is shown in Fig. 8.


[Please insert Figure 7 about here] [Please insert Figure 8 about here] IV) Movement along the direction of the force



After evaluation of the total force Fi , the Update Position procedure is used to move the particle  x i in the direction of the force by a random step length according to Eq. (43). ­ i Fni x + λ (u n − x ni ), if Fni > 0 ° n i F ° x ni = ® , n = 1, 2,..., N i F i i i ° x + λ n ( x − l ), if F ≤ 0 n n n ° n Fi ¯

(43)

The random step length λ is generated with uniform distribution between (0, 1) (Tsou, 2009). In this equation, as the imposed force is normalized and only the direction of the movement is identified, it can be guaranteed that candidate particles do not move to the unvisited particles along this direction with random step length (Jamili et al., 2011). Also, l n and u n are the lower and upper bound for each dimension of particle i, respectively, that denote the allowed feasible movement for the corresponding dimension. The Update Position procedure is shown in Fig. 9.

[Please insert Figure 9 about here]


V) Updating the Pareto archive set After the movement process and local search process, the non-dominated solutions are gathered in the Pareto archive A . Therefore, the new solutions that are not dominated by any of the archive’s members are added to Pareto archive A and any members of A that are dominated  by x i are eliminated from archive A . Accordingly, as the archive set A usually has a limited number of solutions if the number of non-dominated solutions exceeds its limitation, these solutions should be sorted based on their CDs, and those with lower CDs should be eliminated. The procedure of Updating Pareto Archive Set is illustrated in Fig. 10.


[Please insert Figure 10 about here] VI) Fitness function As this model aims to strategically design SCN with two objective functions, total costs and environmental effects, we should consider some affabilities in solving the model. In this regard, to achieve an efficient search through the infeasible region and to assure that the final best solution is feasible, the following penalty terms are appended to the objective functions to consider the capacity limitation of suppliers, plants, distribution centers, and cross-docks, respectively:   f ( x i ) + ( ¦ x s g p (ξ ) − Cap s g ) × 2 ×V C s g , ∀g ∈ G , s g ∈ S g , ξ ∈ Ξ | ¦ x s g p (ξ ) > Cap s g (44) p ∈P

p ∈P

  f ( x i ) + ( ¦ ¦ h pt (ξ ) − Cap pt ) × 2 × MC pt , ∀p ∈ P , ∀t ∈T | ¦ h pt (ξ ) > Cap pt p ∈P t ∈T

d ∈D

  f ( x i ) + ( ¦ x pd (ξ ) − Cap d ) × 2 ×V C d , ∀d ∈ D , ∀ξ ∈ Ξ | ¦ x pd (ξ ) > Capd p ∈P

p ∈P

  f (x i ) + ( ¦ x cr (ξ ) − Capc ) × 2 ×V C c , ∀c ∈C , ∀ξ ∈ Ξ | ¦ x cr (ξ ) > Capc r ∈R

r ∈R

(45) (46) (47)

VII) Local search The proposed AMOEMA is hybridized with a local search to increase the performance of the algorithm. VNS is one of the most recent metaheuristics developed to solve the problems in an easier way (Hansen and Mladenoviü, 2001).The first multi-objective VNS algorithm(MOVNS) was proposed by Geiger (2008). During the VNS process, the local optimum in each neighborhood is discovered iteratively and paves the ground for reaching the global optimum at the end (Vahdani and Zandieh, 2010). MOVNS-related wealthy literature can be found in Geiger (2008) and Arroyo et al. (2011). This metaheuristic has been known and widely applied as an attractive algorithm for many optimization problems, such as transportation problems (Polacek et al., 2004), scheduling (Arroyo et al., 2011; Liang et al., 2009), etc. In this paper, we propose AMOVNS with both intensification and diversification as a reinforced local search (Van

Hentenryck and Vergados, 2007).AMOVNS profits from the ε -domination concept (Mostaghim  and Teich, 2003). In a minimization problem, it is described that a particle x i ε -dominates a    particle x j for some ε >0 ( x i ε x i ) (Moslemi and Zandieh, 2011), if:   f k ( x i ) / (1 + ε ) ≤ f k ( x j ) , ∀ k ∈ {1, 2,..., K }, ∀ i , j ∈ {1, 2,..., m } i≠j (48) i j f k ( x ) / (1 + ε ) < f k ( x ) for at least one k ∈ {1, 2,..., K }, and i , j ∈ {1, 2,..., m } i≠j

i where f k (x ) is the value of the

i k th objective for particle x . AMOVNS allows the search

algorithm to escape from the local minima. Therefore, in spite of the fact that AMOVNS is a local search algorithm, it avoids getting trapped in a local optimum by accepting the worse solutions as well. This approach firstly receives the Pareto front A and population P in each iteration from AMOEMA as the initialization population and chooses one of them by Roulette Wheel Selection Method, with more probability for selecting A . AMOVNS utilizes a set of neighborhood structures (NSs) which are symbolized by N ns (ns = 1, 2,..., ns max ) to find the near-optimal solution, where ns max indicates the maximum number of neighborhood types that their implementation sequence has determined previously. The AMOVNS structure is illustrated in Fig.11. [Please insert Figure 11 about here]


VIII) Neighborhood structures The quality of responses obtained using the local search metaheuristic approach significantly depends on the efficiency of the developed NSs utilized in the neighborhood search process. In this paper, various types of neighborhood search structures in different parts of the solution separately are used separately to generate a new solution. To explain, the structures of applied NSs in the plant echelon as the most complicated sub-solution are described as follows: • Neighborhood structures in the selection process: a) Swap: In this NS, two units of a solution are selected randomly. Then their positions are exchanged. b)Intelligent Swap: Selects the units with maximum and minimum values and then exchanges their positions. c) Reversion: Through this policy, after conducting the swap, the units located between the swapped units are reversed, too. d)Insertion: In this case, the unit in the second position is located immediately after the unit in the first location, and the other units are shifted towards the right accordingly (Rabiee et al., 2012). e)Slide: In this NS, two units of a solution are selected randomly. Then, the first unit is eliminated and others, from its neighbor to the second number, are moved toward the first position. Then the first number is located on the position of the second number.




f) Exchange two neighborhoods (ETN): Selects two neighboring units of a solution randomly, and then swaps their positions. g) Shifting a part of the solution (SPS): In this NS, a subsequence of the solution is selected and then shifted to a new position. h) Shifting a reversed part of the solution (SRPS): This NS is similar to SPS, but inverts the subsequence of the selected part while moving the process. i) Random shuffling (RS): Through this policy a part of a solution with predetermined length is selected, and then the values of its positions are exchanged with others randomly(Vahdani and Zandieh, 2010). The NSs in the sub-solution selection process are illustrated in Fig. 12.
[Please insert Figure 12 about here]




•

Neighborhood structures in the transship rate determination process: a) Exchanging the composition of demands required by lower hand customers from their upper echelon: In this NS, first two columns are selected randomly, and then their positions are exchanged with each other. b) Exchanging the composition of transships supplied by upper hand members for their lower echelon: Through this NS, first two rows are chosen randomly, and then they are exchanged with each other. c) Exchanging the amounts that a lower hand member requires from two upper hand members: In this NS, a random column is selected and two random positions are swapped. d) Exchanging the amounts that two lower hand members require from a specific upper hand member: In this NS, a random row is selected and two random positions in this row are swapped with each other. e) Two positions of the transship matrix are selected randomly and these two positions are swapped with each other. f) Elimination of the semi-idle supplier for a lower hand member: In this NS, first a column is selected randomly, and then the lowest value that this member requires from whole of its upper hand members is changed to zero. The NSs in the sub-solution transship rate determination process are illustrated in Fig. 13.




[Please insert Figure 13 about here]


• Neighborhood structures in the tech-level selection process: The NSs in this sub-solution are considered as the NSs represented in the selection process.


6. Computational results In this section, we are going to examine the proposed multi-objective algorithms, i.e., MOHEV with two strategies (S1: Minimum distance (MD), and S2: Crowding distance (CD)) and compare

it with other existing algorithms, i.e., NSGAII and MOPSO. It should be also mentioned that due to adoptions and innovations stated in the construction of MOEMA that pave the way to achieve more efficient outputs than MOEMA, and due to primary examinations done on MOHEV and MOEMA, it is concluded that MOHEV outperforms MOEMA, so the proposed algorithm in two strategies, i.e. MOHEV(md) and MOHEV(cd), is only compared with NSGAII and MOPSO. The comparisons are performed on the basis of the sets of (non-dominated) solutions obtained by each algorithm. We implemented the algorithms in MATLAB 2009b language and ran on three parallel PCs with 2.4 GHz and 2GB RAM memory.


6.1. Data generation To compare the proposed hybrid algorithm with other benchmark algorithms, an experiment was conducted to deliberate on the performance of different multi-objective metaheuristics. In this regard, a numerical experiment of the problem studied here consists of retailers (R) with their demand in different scenarios, cross-docks (C), distribution centers (D), plants (P) with their different potential technologies (T), and suppliers (S) in different grades as a strategic issue in SCND, is constructed. The test data which are generated for the problem are shown inTable 4.
[Please insert Table 4 about here] Considering the numbers of potential units in each echelon, the problems are classified into three scales. In each scale, four random problems are produced (i.e., totally 12 test problems). Four scenarios are defined for each problem, and accordingly the demands are generated in four distances (−2σ, −σ ) , (−σ,0) , (0,σ ) , and (σ ,2σ ) .The maximum desired number of units in each echelon (i.e., SMax , PMax , DMax and CMax ) are assumed to be half of potential units in each echelon. The following section explains how the opening and fixed costs of each unit ( OCc , OCd , OC p , and FCsg ) are estimated (for example in the cross-docks echelon): OC

mean c

=

α × C Max C

¦ ¦ × C

R

c =1

r =1

TC cr

C ×R

(49)

OC cmin = ª¬(1-δ ) ×OC cmean º¼ and OC cMax = ª¬(1+δ ) ×OC cmean º¼

(50)

OCc =U (OCcmin ,OCcMax )

(51)

Where α is a coefficient with U(40,50) as the opening value of a center with respect to its transportation cost, and δ is a coefficient with U(0,0.1). Moreover, the opening emission of each unit ( EOsg , EOpt , EOd and EOc ) is calculated similar to its opening cost; while transportation emission is used instead of transportation cost. In addition, the capacity of each unit ( Capc ,

Capd , Cappt and Capsg ) is estimated as follows (for example in the cross-docks echelon):

Capcmean =

ζ × (1 + β )

(52)

C Max

Capcmin = ª¬(1-δ ) ×Capcmean º¼ and CapcMax = ª¬(1+δ ) ×Capcmean º¼

(53)

OCc =U (Capcmin ,CapcMax )

(54)

where ζ is the maximum summation of demands required by the lower hand echelon in all scenarios. β is a coefficient with U(0,0.1). 6.2. Evaluation metric To compare the quality of different Pareto-optimal fronts produced by multi-objective optimization methods, some performance assessment metrics are used to make quantitative comparisons between multi-objective metaheuristic algorithms (Behnamian and Fatemi Ghomi, 2011). As these metrics assess the quality of approximated Pareto sets from different perspectives and due to their conflicting and incommensurable nature, various metrics should be used simultaneously (Behnamian et al., 2009). Accordingly, in this paper the following four known performance assessment metrics are used: • Diversification metric (DM): Is used to measure the spread of the non-dominated solution set (Hyun et al., 1998). In this metric, the algorithm with a higher value has a better capability. It can be expressed by the following equation: n

D=

¦ max ( x ′ − y ′ )
i

i

(55)

i =1

where n is the number of Pareto front members and x i′ − y i′ is the Euclidean distance between the non-dominated solutions, xi and y i . •

Spread of non-dominant solution (SNS): Is considered as a diversity measure and evaluates the standard deviation of distance of an ideal point from a/an non-dominated set (Behnamian and Fatemi Ghomi, 2011). The higher value of SNS brings the better solution quality. SNS can be calculated as:

¦ i =1 (c − c i ) n

SNS =

2

(56)

n −1    c  
c i = f i − f Ideal c = i n
f Ideal = {min(f 1 ), min(f 2 ),..., min(f k )}
•

Percent of domination (POD) was proposed on the basis of the coverage metric (CM) introduced by (Zitzler and Thiele, 1998). CM measures the extent to which one Pareto set X ′′ is covered by another solution set X ′ by comparing the number of solutions in X ′′ covered by solutions in X ′ . In this regard, if the number of algorithms is more than two, all pair combinations should be analyzed. To calculate the POD, all of the non-dominated solutions acquired by algorithms are combined in one Pareto set, and then the percentage

of solutions belonging to each algorithm is calculated. The algorithm with higher POD has better performance. The equation of POD is represented as follows:


POD (X ′, X ′′,..., X n ) = •

{x ki ∈ X i } ∃ x kj ∈ X j , i ∈ j : x ki X

≤ x kj

× 100

(57)

j

Data envelopment analysis (DEA) developed by Charnes et al. (1978) is usually applied to evaluate the performance of some choices with some attributes compared to each other (Wu and Blackhurst, 2009). A more detailed theoretical introduction of DEA can be found in Cooper (2004). In this study, we utilized this method to determine the efficiency of non-dominated solutions obtained by each method. To calculate this metric, each Pareto optimal solution can be considered as a DMU. Moreover, an objective (e.g., total costs) may be considered as input, while the other objective (e.g., environmental effects) is considered as output. Subsequently, all non-dominated solutions obtained by the algorithms are combined and the efficiency of these points is calculated by the DEA model, which introduced by Amin (2009); Amin and Toloo (2007). Clearly the algorithm with higher value has better efficiency than the other algorithms.




6.3. Parameter setting Parameter setting is identified as an important milestone to achieve robust metaheuristic algorithms by not producing functional variance under external environmental influence (Gholami et al., 2009). This process intends to find good values for the parameters before the run of the algorithm, which remains fixed during the run (Behnamian et al., 2010). Generally speaking, the efficiency of the algorithms is heavily based on miscellaneous factors, such as the assigned values to parameters (Vahdani and Zandieh, 2010).In this regard, this paper employs the response surface method (RSM) as a well-known technique that is widely applied in a variety of industrial settings and parameter optimizations, to determine the optimal values of the parameters (Rabiee et al., 2012). RSM is a collection of mathematical and statistical techniques developed by Box and Wilson in the early 1950s (Myers et al., 1971 and 2009). The response surface can be expressed as y = f (x 1, x 2 ,…, x k ) , while the variables are assumed measurable, k is the number of variables and the goal is to optimize the response variable Y. The variables should be independent, continuous, and controllable through experiments with negligible errors. Finding a suitable approximation for the true functional relationship between independent variables and the response surface is prerequisite. Neter et al. (1996) proposed a polynomial response surface function as: k

k

y = β0 + ¦β j x j + ¦¦βij xi x j + ¦β jj x2j + ε
j =1

i< j

j =1

(58)

where y is the predicted response; β0 is a constant; β j is the linear coefficient; βij is the crossproduct coefficient; and β jj is the squared coefficient. Eq.(58) should be considered when there is a curvature in the system. Considering k and two lower and upper bounds for each variable, the number of experiments is equal to (2k) or a fraction of it. The factors, their levels, and the k number of experiments for each algorithm ( nf = 2 , nax =number of axial points (2k) coded ±1 ,

nCp =number of central points) are presented in Table 5. It should be noticed that the effect of different strategies in parameter setting of MOHEV, i.e., (cd) and (md), is neutralized in our experiments by considering them as two different blocks. To compare the algorithms fairly, the population of non-dominated frontiers was set constant (i.e., nPop in NSGAII, nRep in MOPSO, and PopF in MOHEV are considered 50).


[Please insert Table 5 about here] The variables can be coded as Eq. (59) in which x i is the coded value of the i th test variable and X i is the uncoded value of it. xi =

X i − ((X upper + X Lower ) / 2) (59)

(X upper − X Lower ) / 2

As some performance assessment metrics are used together to comparethe multi-objective metaheuristics, a tuning problem is formedto set the input factor’s value based on represented performance metrics, simultaneously, which leads to the desired response in multi-objective decision making (MODM) environments. Accordingly, Derringer and Suich (1980) provided a utility function to optimize multipleresponses as follows: S

§ H −Y i · d i (Y i ) = ¨ i ¸ Li ≤ Y i ≤ H i © H i − Li ¹

(60)

where, response functions in the form of minimum Y(i) are converted to the utility function

d i (Y i ) . Li and H i are lower and upper bounds, respectively. S is the severity of the utility function by which the shape of the desirability can be changed for each goal according to severity value. Severities are used to put emphasis on the upper/lower bounds. If S =1, d i (Y i ) varies from 0 to 1 in a linear fashion. When S is between 1 and 10, it means that there is more emphasis on the goal, and conversely, when S is between 0.1 and 1, there is less emphasis on the goal as illustrated in Fig. 14. Parameters S are set to 1, 1, 2, and 2 for Diversity, SNS, POD and DEA, due to their considerable importance, respectively. POD and DEA were calculated by comparing the total experiments results of all algorithms. The desirability function is as follows:

( )

( )

D = m d1 y 1 × d 2 y 2 × … × d m

(y m )

(61)

where m is the number of objectives. The tuned values for parameters, R-Squared (R2) and desirability (D) are estimated as displayed in Table 6. [Please insert Figure 14 about here] [Please insert Table 6 about here]


6.4. Comparisons among algorithms In this section, the effectiveness of the algorithms is compared by using the proposed metrics as comparison measures. In order to determine whether there is a significant difference among the performance of algorithms after the calculation of Pareto optimal solutions for all four algorithms for each problem, the efficiency of each algorithm is assessed in terms of DM, SNS, POD, and DEA. According to Table 7, which presents results obtained for each problem, MOHEV(cd) clearly exceeds other approaches, especially in medium and large size problems in all metrics, while in small problems, although MOHEV(cd) had a minor preference over others, there was not any obvious priority among the algorithms in this size.

[Please insert Table 7 about here]


Fig. 15 demonstrates the non-dominated solutions of a single run obtained by algorithms in three scales. As it can be seen, along with the increase in the size of problems that lead to the complexity of problems, the increasing effectiveness of MOHEV(cd) can be observed, clearly. Also, while the size of the problems is increasing, the difference between MOHEV(cd) and MOHEV(md) increases compared to the other algorithms, i.e., NSGAII and MOPSO. [Please insert Figure 15 about here] In addition, to find the best algorithm decisively, a statistical comparison among the algorithms in each performance measure is done. Accordingly, first of all, the results obtained for each problem are transformed to Relative Percentage Deviation(RPD) obtained by the following formula: ALg sol − Min sol RPD = (62) Min sol
where, ALg sol is the objective value obtained by a given algorithm for a problem, and Minsol is the best solution obtained for that problem. The lower value of RPD is preferred. Then, a Kruskal–Wallis analysis of the variance (ANOVA) as a non-parametric method has been performed for testing whether the samples originate from the same distribution. The obtained results verify that at least two algorithms are not identical in all metrics ( DM, H = 17.67, P = 0.001 ; SNS, H = 21.78, P = 0.00; POD, H = 25.82, P = 0.00 , and

DEA, H = 25.82, P = 0.00 ).

In addition, Fig. 16 indicates that in SNS, POD, and DEA metrics, MOHEV is significantly better than the other algorithms. Also, according to SNS, POD, and DEA metrics, it can be declared that by growing the dimensions of the problems, the efficiency of the proposed algorithm outperformsthe other approaches. [Please insert Figure 16 about here]


7. Application of the proposed approach in a case study In this section, a real case study is provided to demonstrate the usefulness of the proposed model. The case is an Iranian automobile producer named Iran Khodro Co. (IK) as the largest automobile factory in the Middle East which is planning to redesign its SCN and to develop a network due to the increasing demand of its customers. This company aims to preserve its lion share in the internal market, with more than 50% of the total market share, and even improve it in the other countries of the Middle East. IK has formed a supply network based on its plants as the main part of the network and distributes the products to local retailers. The supplier management is the most important concern of this network. The suppliers try to exceed each other to gain a place in a better class to benefit from financial advantages, volume of purchase, etc. In this rat race, the IK supplier experts audit the suppliers by questionnaires designed based on economical and environmental aspects. The network has two active manufacturing centers in two plants that can be relocated and three potential ones considered to be opened as new plant(s). Also, there are two existing distribution centers and four potential locations which can be selected to serve thirteen retailers scattered throughout the Middle East and especially inside the country. Also, it is planned to use the cross-docks in the network to reduce the cost and environmental emission. Accordingly, three potential cross-docks are considered in this redesign process. Also, nine suppliers were audited and placed in four different classes A, B, C, and D, and are ready to sell their raw materials/automobile parts to the plants. The plants can use different production lines as technologies that each of them uses require different levels of investment and variable cost and the also have special impacts on the environment. The required data are presented in Tables 8 through 10. It is noteworthy that the fixed opening cost and emission of active plants and distribution centers are equal to zero. The maximum desired number of plants, distribution centers, and cross-docks are 4, 3, and 2, respectively. The suppliers are classified based on the resulting score of the auditing process to different classes as (0.75, 100] Class A, (0.6, 0.75] Class B, (0.5, 0.6] Class C, and [0, 0.5] Class D. The suppliers placed in class D are ignored by IK as a qualified supplier. Also, it should be mentioned that the Iranian rial is considered as the monetary unit. [Please insert Table 8 about here] [Please insert Table 9 about here]

The transportation cost, i.e., TC s g p , TC pd , TC dc , TC cr and TC dr , is estimated based on the distances between each two facilities. Furthermore, the per unit environmental impacts of transportation, i.e., EM pt , ET s g p , ET pd , ET dc , ET cr and ET dr , are calculated with a coefficient of this distance. According to Standard Unit Fuel Consumption (SUFC), one gram of CO2 released by gasoline diesel per kilometer can be estimated to liter per 100 kilometers as 1 gram/kilometer CO2 equals 0.043 liter/100 kilometer. In other words, each track produces 23.2 grams of CO2 for each kilometer driven. It should ne noted that due to space limitation, the value of some input parameters such as distances between facilities calculated based on the position of facilities is not provided and can be presented upon request. [Please insert Table 10 about here] To solve the case problem, the MOHEV(cd) algorithm, as the best proposed approach able to handle different sizes of the sustainable SCN design problem, is utilized. MOHEV(cd) is constructed such as to be able to find a good approximation of the Pareto front for the decisionmaker. The results are illustrated in Fig. 17 and reported in Table 11. According to the achieved results, when the decision maker is cost-oriented, the centers with lower opening and operating costs are selected. However, the environment-oriented decision-makers select centers which create less ramifications on environment. [Please insert Figure 17 about here] [Please insert Table 11 about here]
8. Conclusions and future work In this paper a novel hybrid optimization approach is applied to strategically design a sustainable forward supply chain network (SCN) with stochastic demand in order to minimize the total costs and environmental effects, simultaneously. In the proposed model, the impacts of suppliers on cost and environment objectives are considered by integrating sourcing decisions into the model. Due to the complexity of the problem, a new metaheuristic algorithm called MOHEV is used to generate Pareto-optimal solutions. This algorithm combines the adaptive multi-objective electromagnetism-like mechanism algorithm (AMOEMA) and the adaptive multi-objective variable neighborhood search (AMOVNS) as a local search is proposed. The performance and the reliability of the proposed algorithms were evaluated in comparison with available benchmark algorithms (i.e., NSGAII and MOPSO). Four different comparison metrics, i.e., Diversity of solutions, Spread of non-dominated solutions, Percent of domination, and Data envelopment analysis (DEA) were utilized to compare the efficiency of the algorithms. For a fair comparison, all the algorithms’ parameters were estimated by the response surface method (RSM) and the multi-objective decision making (MODM) approach based on four different comparison metrics. In this regard, for the MOHEV algorithm, the input parameters were estimated in their upper limit, i.e., coded parameters of the maximum number of iteration

(Maxit) and number of population (nPop) were tuned as 0.94 and 0.8, respectively. Also, the coded parameter of the local search iteration (LSiter) was estimated 1. Accordingly, the results obtained by the proposed algorithm almost dominate the results of NSGAII and MOPSO. Based on performance metrics resulting from DEA on all problems, the average performance of MOHEV(cd), MOHEV(md), MOPSO and NSGAII were calculated as 0.35, 0.295, 0.152, and 0.204, respectively. Moreover, the sensitivity analysis for the performance of algorithms versus the dimensions of the problem demonstrated that MOHEV yielded better solutions in comparison with NSGAII and MOPSO, and also MOHEV(cd) is superior to MOHEV(md) in most situations. For further research, according to achieved results, it is recommended to assess the performance of the proposed method in other optimization problems. In addition, regarding the proposed model, it is recommended to extend the model to a closed-loop supply chain through appending collecting, recycling, and disposal centers. In this model, the economic and environmental effects of suppliers are encountered in the supply chain effects; a possible future work is integrating the grading, classification and selection process explicitly in the model. Finally, the vehicle routing of can be added to achieve a more detailed optimization model.


References Ageron, B., Gunasekaran, A., Spalanzani, A., 2011. Sustainable Supply Management: An Empirical Study. International Journal of Production Economics. 140 (1), 168-182. Amin, G.R., 2009. Comments on finding the most efficient DMUs in DEA: An improved integrated model. Computers & Industrial Engineering. 56 (4), 1701-1702. Amin, G.R., Toloo, M., 2007. Finding the most efficient DMUs in DEA: An improved integrated model. Computers & Industrial Engineering. 52 (1), 71-77. Arroyo, J.E.C., dos Santos Ottoni, R., de Paiva Oliveira, A., 2011. Multi-objective Variable Neighborhood Search Algorithms for a Single Machine Scheduling Problem with Distinct due Windows. Electronic Notes in Theoretical Computer Science. 281 5-19. Bachlaus, M., Pandey, M., Mahajan, C., Shankar, R., Tiwari, M., 2008. Designing an integrated multiechelon agile supply chain network: a hybrid taguchi-particle swarm optimization approach. Journal of Intelligent Manufacturing. 19 (6), 747-761. Bai, C., Sarkis, J., 2010. Integrating sustainability into supplier selection with grey system and rough set methodologies. International Journal of Production Economics. 124 (1), 252-264. Basurko, O.C., Mesbahi, E., 2012. Methodology for the sustainability assessment of marine technologies. Journal of Cleaner Production. 68 (1), 155–164. Baumgartner, R.J., 2011. Critical perspectives of sustainable development research and practice. Journal of Cleaner Production. 19 (8), 783-786. Behnamian, J., Fatemi Ghomi, S., 2011. Hybrid flowshop scheduling with machine and resourcedependent processing times. Applied Mathematical Modelling. 35 (3), 1107-1123. Behnamian, J., Fatemi Ghomi, S., Zandieh, M., 2009. A multi-phase covering Pareto-optimal front method to multi-objective scheduling in a realistic hybrid flowshop using a hybrid metaheuristic. Expert Systems with Applications. 36 (8), 11057-11069. Behnamian, J., Ghomi, S., 2010. Hybrid flowshop scheduling with machine and resource-dependent processing times. Applied Mathematical Modelling. Behnamian, J., Zandieh, M., Fatemi Ghomi, S.M.T., 2010. Due windows group scheduling using an effective hybrid optimization approach. The International Journal of Advanced Manufacturing Technology. 46 (5), 721-735.

Bergenwall, A.L., Chen, C., White, R.E., 2012. TPS's process design in American automotive plants and its effects on the triple bottom line and sustainability. International Journal of Production Economics. 140 (1), 374-384. Bhattacharya, C.B., Sen, S., 2004. When, why, and how consumers respond to social initiatives. California Management Review. 47 (1), 9–24. Birbil, ù.ø., Fang, S.C., 2003. An electromagnetism-like mechanism for global optimization. Journal of Global Optimization. 25 (3), 263-282.

Brandenburg, M., Govindan, K., Sarkis, J., and S. Seuring. (2014). Quantitative models for sustainable supply chain management: developments and directions, European Journal of Operational Research, Vol. 233, No. 2, pp. 299-312
Brundtland, G.H., 1987. Our common future. Oxford paperbacks. (A/42/427). Carter, C.R., Rogers, D.S., 2008. A framework of sustainable supply chain management: moving toward new theory. International Journal of Physical Distribution & Logistics Management. 38 (5), 360-387. Chaabane, A., Ramudhin, A., Paquet, M., 2012. Design of sustainable supply chains under the emission trading scheme. International Journal of Production Economics. 135 (1), 37-49. Chang, P.C., Chen, S.H., Fan, C.Y., 2009. A hybrid electromagnetism-like algorithm for single machine scheduling problem. Expert Systems with Applications. 36 (2), 1259-1267. Charnes, A., Cooper, W.W., Rhodes, E., 1978. Measuring the efficiency of decision making units. European Journal of Operational Research. 2 (6), 429-444. Coello Coello, C.A., Lechuga, M.S., 2002. MOPSO: A proposal for multiple objective particle swarm optimization. IEEE, 1051-1056 Cooper, W.W., 2004. Handbook on Data Envelopment Analysis (International Series in Operations research & Management Science) POD. Cruz, J., Matsypura, D., 2009. Supply chain networks with corporate social responsibility through integrated environmental decision-making. International Journal of Production Research. 47 (3), 621648. Cruz, J.M., 2008. Dynamics of supply chain networks with corporate social responsibility through integrated environmental decision-making. European Journal of Operational Research. 184 (3), 10051031. Cruz, J.M., 2009. The impact of corporate social responsibility in supply chain management: Multicriteria decision-making approach. Decision Support Systems. 48 (1), 224-236. Dao, V., Langella, I., Carbo, J., 2011. From green to sustainability: Information Technology and an integrated sustainability framework. The Journal of Strategic Information Systems. 20 (1), 63-79. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: NSGA-II. Evolutionary Computation, IEEE Transactions on. 6 (2), 182-197. Derringer, G., Suich, R., 1980. Simultaneous optimization of several response variables. Journal of quality technology. 12 (4), 214-219. Du, F., Evans, G.W., 2008. A bi-objective reverse logistics network analysis for post-sale service. Computers & operations research. 35 (8), 2617-2634. Elhedhli, S., Merrick, R., 2012. Green supply chain network design to reduce carbon emissions. Transportation Research Part D: Transport and Environment. 17 (5), 370-379. Elkington, J., 1997. Cannibals with Forks: The Triple Bottom Line of Twenty First Century Business. Capstone, Mankato, MN. Elkington, J., 1998. Cannibals with Forks: The Triple Bottom Line of 21st Century Business (The Conscientious Commerce Series). Erol, I., Sencer, S., Sari, R., 2011. A new fuzzy multi-criteria framework for measuring sustainability performance of a supply chain. Ecological Economics. 70 (6), 1088–1100. Fleischmann, M., Bloemhof-Ruwaard, J.M., Beullens, P., Dekker, R., 2004. Reverse Logistics Network Design, in: Dekker, R., Fleischmann, M., Inderfurth, K., Van Wassenhove, L.N. (Eds.), Reverse logistics: quantitative models for closed-loop supply chains. Springer, Berlin, 65–94

Fleischmann, M., Kuik, R., Dekker, R., 2002. Controlling inventories with stochastic item returns: A basic model. European Journal of Operational Research. 138 (1), 63-75. Fleischmann, M., Van Wassenhove, L.N., van Nunen, J.A.E.E., van der Laan, E.A., Dekker, R., Bloemhof-Ruwaard, J.M., 1997. Quantitative models for reverse logistics: a review. European Journal of Operational Research. 103 (1), 1–17. Fombrun, C.J., 2005. The leadership challenge: Building resilient corporate reputations. Handbook on responsible leadership and governance in global business. 54 68. Frota Neto, J.Q., Bloemhof-Ruwaard, J.M., van Nunen, J.A.E.E., van Heck, E., 2008. Designing and evaluating sustainable logistics networks. International Journal of Production Economics. 111 (2), 195-208. Geiger, M.J., 2008. Randomised Variable Neighbourhood Search for Multi Objective Optimisation. Arxiv preprint arXiv:0809.0271. Gholami, M., Zandieh, M., Alem-Tabriz, A., 2009. Scheduling hybrid flow shop with sequencedependent setup times and machines with random breakdowns. The International Journal of Advanced Manufacturing Technology. 42 (1), 189-201. Gimenez, C., Sierra, V., Rodon, J., 2012. Sustainable operations: Their impact on the triple bottom line. International Journal of Production Economics. 140 (1), 149-159. Godfrey, R., 1998. Ethical purchasing: developing the supply chain beyond the environment. Greener purchasing: Opportunities and innovations. 244-251. Govindan, K., Khodaverdi, R., Jafarian, A., 2013. A Fuzzy Multi criteria approach for measuring sustainability performance of a Supplier based on triple bottom line approach. Journal of Cleaner Production. 47 345-354.

Govindan, K., & Popiuc, M. N. (2014). Reverse supply chain coordination by revenue sharing contract: A case for the personal computers industry.European Journal of Operational Research, 233(2), 326-336. Govindan, K., Jafarian, A., Khodaverdi, R., and K. Devika (2014a), Two-echelon multiplevehicle location-routing problem with time windows for optimization of sustainable supply chain network of perishable food, International Journal of Production Economics, Vol. 152, pp. 9-28 Govindan, K., Soleimani, H., & Kannan, D. (2015). Reverse logistics and closed-loop supply chain: A comprehensive review to explore the future. European Journal of Operational Research, 240(3), 603-626. Gunasekaran, A., Spalanzani, A., 2011. Sustainability of manufacturing and services: Investigations for research and applications. International Journal of Production Economics. Hansen, P., Mladenoviü, N., 2001. Variable neighborhood search: Principles and applications. European Journal of Operational Research. 130 (3), 449-467. Hassini, E., Surti, C., Searcy, C., 2012. A literature review and a case study of sustainable supply chains with a focus on metrics. International Journal of Production Economics. 140 (1), 69–82. Hsu, C.-W., Hu, A.H., 2009. Applying hazardous substance management to supplier selection using analytic network process. Journal of Cleaner Production. 17 (2), 255-264. Hsueh, C.F., Chang, M.S., 2008. Equilibrium analysis and corporate social responsibility for supply chain integration. European Journal of Operational Research. 190 (1), 116-129. Hu, C., Li, S., 2009. Two-phase interactive satisfying method of fuzzy multiple objective optimization with linguistic preference. International Journal of Information Technology & Decision Making. 8 (3), 427-443. Hyun, C.J., Kim, Y., Kim, Y.K., 1998. A genetic algorithm for multiple objective sequencing problems in mixed model assembly lines. Computers and Operations Research. 25 (7-8), 675-690.

Jamili, A., Shafia, M., Tavakkoli-Moghaddam, R., 2011. A hybridization of simulated annealing and electromagnetism-like mechanism for a periodic job shop scheduling problem. Expert Systems with Applications. 38 (5), 5895-5901. Jo, J.-B., Li, Y., Gen, M., 2007. Nonlinear fixed charge transportation problem by spanning tree-based genetic algorithm. Computers & Industrial Engineering. 53 (2), 290-298.

Kannan, D., Diabat, A., Alrefaei, M., Govindan, K., & Yong, G. (2012). A carbon footprint based reverse logistics network design model. Resources, Conservation and Recycling, 67, 7579. Kannan, G., Noorul Haq, A., & Devika, M. (2009). Analysis of closed loop supply chain using genetic algorithm and particle swarm optimisation.International Journal of Production Research, 47(5), 1175-1200. Khalili, M., Tavakkoli-Moghaddam, R., 2012. A multi-objective electromagnetism algorithm for a biobjective flowshop scheduling problem. Journal of Manufacturing Systems. 31 (2), 232–239. Ko, H.J., Evans, G.W., 2007. A genetic algorithm-based heuristic for the dynamic integrated forward/reverse logistics network for 3PLs. Computers & operations research. 34 (2), 346-366. Li, Z., Liao, H., Coit, D.W., 2009. A two-stage approach for multi-objective decision making with applications to system reliability optimization. Reliability engineering & system safety. 94 (10), 15851592. Liang, Y.C., Chen, A.H.L., Tien, C.Y., 2009. Variable Neighborhood Search for Multi-Objective Parallel Machine Scheduling Problems, 519-522 Lieckens, K., Vandaele, N., 2007. Reverse logistics network design with stochastic lead times. Computers and Operations Research. 34 (2), 395–416. Linton, J.D., Klassen, R., Jayaraman, V., 2007. Sustainable supply chains: An introduction. Journal of Operations Management. 25 (6), 1075-1082. Moslemi, H., Zandieh, M., 2011. Comparisons of some improving strategies on MOPSO for multiobjective (r, Q) inventory system. Expert Systems with Applications. 38 (10), 12051-12057. Mostaghim, S., Teich, J., 2003. The role of ε-dominance in multi objective particle swarm optimization methods. IEEE, 1764-1771 Vol. 1763 Myers, R.H., Montgomery, D.C., Anderson-Cook, C.M., 1971. Response surface methodology. Allyn and Bacon Boston. Myers, R.H., Montgomery, D.C., Anderson-Cook, C.M., 2009. Response surface methodology: process and product optimization using designed experiments. John Wiley & Sons Inc. Naderi, B., Tavakkoli-Moghaddam, R., Khalili, M., 2010. Electromagnetism-like mechanism and simulated annealing algorithms for flowshop scheduling problems minimizing the total weighted tardiness and makespan. Knowledge-Based Systems. 23 (2), 77-85. Neter, J., Wasserman, W., Kutner, M.H., 1996. Applied linear regression models. Irwin Homewood. Nydick, R.L., Hill, R.P., 1992. Using the analytic hierarchy process to structure the vendor selection procedure. International Journal of Purchasing and Materials Management. 28 (2), 31-36. Pagell, M., Krause, D., Klassen, R., 2008. Sustainable supply chain management: theory and practice. Journal of Supply Chain Management. 44 (1), 85-85. Pinto-Varela, T., Barbosa-Póvoa, A.P.F.D., Novais, A.Q., 2011. Bi-objective optimization approach to the design and planning of supply chains: Economic versus environmental performances. Computers & Chemical Engineering. 35 (8), 1454-1468. Pishvaee, M., Kianfar, K., Karimi, B., 2010a. Reverse logistics network design using simulated annealing. The International Journal of Advanced Manufacturing Technology. 47 (1), 269-281. Pishvaee, M.S., Farahani, R.Z., Dullaert, W., 2010b. A memetic algorithm for bi-objective integrated forward/reverse logistics network design. Computers & operations research. 37 (6), 1100-1112. Pishvaee, M.S., Razmi, J., 2012. Environmental supply chain network design using multi-objective fuzzy mathematical programming. Applied Mathematical Modelling. 36 (8), 3433-3446.

Pishvaee, M.S., Razmi, J., Torabi, S.A., 2012. Robust possibilistic programming for socially responsible supply chain network design: A new approach. Fuzzy Sets and Systems. 206 (0), 1-20. Polacek, M., Hartl, R.F., Doerner, K., Reimann, M., 2004. A variable neighborhood search for the multi depot vehicle routing problem with time windows. Journal of Heuristics. 10 (6), 613-627. Qiang, Q., Ke, K., Anderson, T., Dong, J., 2013. The closed-loop supply chain network with competition, distribution channel investment, and uncertainties. Omega. 41 (2), 186-194. Rabiee, M., Zandieh, M., Jafarian, A., 2012. Scheduling of a no-wait two-machine flow shop with sequence-dependent setup times and probable rework using robust meta-heuristics. International Journal of Production Research. 50 (24), 7428-7446. Robeson, J.F., Copacino, W.C., Howe, R.E., 1992. The Logistics Handbook. Macmillan, Inc., NewYork. Salema, M.I.G., Barbosa-Povoa, A.P., Novais, A.Q., 2007. An optimization model for the design of a capacitated multi-product reverse logistics network with uncertainty. European Journal of Operational Research. 179 (3), 1063-1077.

Sasikumar, P., and G. Kannan, (2009) “Issues in reverse supply chain, part III: Classification and simple analysis”, International Journal of Sustainable Engineering, Vol.2, No.1, pp.2 – 27 Sasikumar, P., and G. Kannan, (2008a) “Issues in reverse supply chains, Part-I: End-of-life product recovery and inventory management-An overview”, International Journal of Sustainable Engineering, Vol. 1, No. 3, pp. 154-172. Sasikumar, P., and G. Kannan, (2008b) “Issues in reverse supply chain, part II: reverse distribution issues – an overview”, International Journal of Sustainable Engineering, Vol. 1, No. 4, pp. 234-249. Soleimani, H and Govindan, K., (2014) Reverse Logistics Network Design and Planning Utilizing Conditional Value at Risk, European Journal of Operational Research , Vol. 237, No. 2, pp. 487-497 Sarkis, J., Talluri, S., 2002. A Model for Strategic Supplier Selection. Journal of Supply Chain Management. 38 (1), 18-28. Seuring, S., 2013. A review of modeling approaches for sustainable supply chain management. Decision Support Systems. 54 (4), 1513-1520. Seuring, S., Gold, S., 2013. Sustainability management beyond corporate boundaries: from stakeholders to performance. Journal of Cleaner Production. 56 (1), 1-6. Seuring, S., Müller, M., 2008. From a literature review to a conceptual framework for sustainable supply chain management. Journal of Cleaner Production. 16 (15), 1699-1710. Shen, L., Olfat, L., Govindan, K., Khodaverdi, R., Diabat, A., 2012. A fuzzy multi criteria approach for evaluating green supplier's performance in green supply chain with linguistic preferences. Resources, Conservation and Recycling. 74 (1), 170–179. Shen, Z., 2007. Integrated supply chain design models: a survey and future research directions. Journal of Industrial and Management Optimization. 3 (1), 1. Shen, Z.J.M., Daskin, M.S., 2005. Trade-offs between customer service and cost in integrated supply chain design. Manufacturing & Service Operations Management. 7 (3), 188-207. Smith, L., Ball, P., 2012. Steps towards sustainable manufacturing through modelling material, energy and waste flows. International Journal of Production Economics. 140 (1), 227-238. Spar, D.L., La Mure, L.T., 2003. The power of activism: Assessing the impact of NGOs on global business. California Management Review. 45 (3), 78-101. Srivastava, S.K., 2007. Green supply chain management: a state of the art literature review. International Journal of Management Reviews. 9 (1), 53-80. Srivastava, S.K., 2008. Network design for reverse logistics. Omega. 36 (4), 535-548. Su, C.T., Lin, H.C., 2011. Applying electromagnetism-like mechanism for feature selection. Information sciences. 181 (5), 972-986.

Syarif, A., Yun, Y., Gen, M., 2002. Study on multi-stage logistic chain network: a spanning tree-based genetic algorithm approach. Computers & Industrial Engineering. 43 (1–2), 299-314. Tang, C.S., Zhou, S., 2012. Research advances in environmentally and socially sustainable operations. European Journal of Operational Research. 223 (3), 585-594. Tavakkoli-Moghaddam, R., Khalili, M., Naderi, B., 2009. A hybridization of simulated annealing and electromagnetic-like mechanism for job shop problems with machine availability and sequencedependent setup times to minimize total weighted tardiness. Soft Computing-A Fusion of Foundations, Methodologies and Applications. 13 (10), 995-1006. Teuteberg, F., Wittstruck, D., 2010. A Systematic Review of Sustainable Supply Chain Management. Multikonferenz Wirtschaftsinformatik 2010. 203. Tseng, M.-L., Tan, R.R., Siriban-Manalang, A.B., 2013. Sustainable consumption and production for Asia: sustainability through green design and practice. Journal of Cleaner Production. 40 (1), 1-5. Tsou, C.S., 2009. Evolutionary Pareto optimizers for continuous review stochastic inventory systems. European Journal of Operational Research. 195 (2), 364-371. Tsou, C.S., Hsu, C.H., Yu, F.J., 2008. Using multi-objective electromagnetism-like optimization to analyze inventory tradeoffs under probabilistic demand. Journal of Scientific & Industrial Research. 67 (8), 569-573. Tsou, C.S., Kao, C.H., 2006. An electromagnetism-like meta-heuristic for multi-objective optimization. IEEE, 1172-1178 Tsou, C.S., Kao, C.H., 2008. Multi-objective inventory control using electromagnetism-like metaheuristic. International Journal of Production Research. 46 (14), 3859-3874. Vachon, S., Klassen, R.D., 2008. Environmental management and manufacturing performance: The role of collaboration in the supply chain. International Journal of Production Economics. 111 (2), 299-315. Vahdani, B., Zandieh, M., 2010. Scheduling trucks in cross-docking systems: Robust meta-heuristics. Computers & Industrial Engineering. 58 (1), 12-24. Van Hentenryck, P., Vergados, Y., 2007. Population-based simulated annealing for traveling tournaments. Menlo Park, CA; Cambridge, MA; London; AAAI Press; MIT Press, 267-272 Wang, F., Lai, X., Shi, N., 2011. A multi-objective optimization for green supply chain network design. Decision Support Systems. 51 (2), 262-269. Wang, H.F., Hsu, H.W., 2010. A closed-loop logistic model with a spanning-tree based genetic algorithm. Computers & operations research. 37 (2), 376-389. Weber, C.A., 1996. A data envelopment analysis approach to measuring vendor performance. Supply Chain Management: An International Journal. 1 (1), 28-39. Weber, C.A., Current, J.R., Benton, W.C., 1991. Vendor selection criteria and methods. European Journal of Operational Research. 50 (1), 2-18. Wu, T., Blackhurst, J., 2009. Supplier evaluation and selection: an augmented DEA approach. International Journal of Production Research. 47 (16), 4593-4608. Yao, M.J., Hsu, H.W., 2009. A new spanning tree-based genetic algorithm for the design of multi-stage supply chain networks with nonlinear transportation costs. Optimization and Engineering. 10 (2), 219237. Zhang, Y., Rao, Y., Zhou, M., 2007. GASA Hybrid Algorithm Applied in Airline Crew Rostering System. Tsinghua Science & Technology. 12 (1), 255-259. Zhu, Q., Sarkis, J., 2004. Relationships between operational practices and performance among early adopters of green supply chain management practices in Chinese manufacturing enterprises. Journal of Operations Management. 22 (3), 265-289. Zitzler, E., Thiele, L., 1998. An evolutionary algorithm for multiobjective optimization: The strength Pareto approach. Swiss Federal Institute of Technology, TIK-Report. 43.





 
 


 
 
Sets G

Description set of suppliers’ grades g ∈ G

Sg

set of potential suppliers in grade g s g ∈ S g

P

set of potential plants p ∈ P

D C

set of potential distribution centers d ∈ D set of potential cross-docks c ∈ C set of retailers r ∈ R set of production technologies t ∈T set of demand scenarios ξ ∈ Ξ

R T

ξ






 


 
 
Parameters TEC pt

Description acquisition cost of technology t ∈T at plant p ∈ P

OC p

opening cost of plant p ∈ P

OC d

opening cost of distribution center d ∈ D

OC c

opening cost of cross-dock c ∈ C

FC s g

fixed costs of long relationship with supplier s g ∈ S g

V Cd

variable cost for handling products at distribution center d ∈ D

V Cc

variable cost for handling products at cross-dock c ∈ C

TC s g p

per unit transportation cost from supplier s g ∈ S g to plant p ∈ P

TC pd

per unit transportation cost from plant p ∈ P to distribution center d ∈ D

TC dc

per unit transportation cost from distribution center d ∈ D to cross-dock c ∈ C

TC cr

per unit transportation cost from cross-dock c ∈ C to retailer r ∈ R

TC dr

per unit transportation cost from distribution center d ∈ D to retailer r ∈ R

MC pt

manufacturing cost of each unit product at plant p ∈ P with technology t ∈T

EM pt

manufacturing emissions of each unit product at plant p ∈ P with technology t ∈T

ET s g p

per unit environmental impacts of transporting from supplier s g ∈ S g to plant p ∈ P

ET pd

per unit environmental impacts of transporting from plant p ∈ P to distribution center d ∈ D

ET dc

per unit environmental impacts of transporting from distribution center d ∈ D to cross-dock c ∈ C

ET cr

per unit environmental impacts of transporting from cross-dock c ∈ C to retailer r ∈ R

ET dr

per unit environmental impacts of transporting from distribution center d ∈ D to retailer r ∈ R

ES s g

environmental impacts related to supplier s g ∈ S g

EO pt

environmental impacts of opening plant p ∈ P with technology t ∈T

EOd

environmental impacts of opening distribution center d ∈ D

EO c

environmental impacts of opening cross-dock c ∈ C

Caps g

capacity of supplier s g ∈ S g

Cap pt

production capacity of plant p ∈ P with technology t ∈T

Capd

handling capacity of distribution center d ∈ D

Capc

handling capacity of cross-dock c ∈ C

S Max

maximum desired number of suppliers

PMax

maximum desired number of plants

D Max

maximum desired number of distribution centers

C Max

maximum desired number of cross-docks

P (ξ )

probability of scenario ξ ∈ Ξ

Dem r (ξ )

demand of retailer r ∈ R under scenario ξ ∈ Ξ

























 
 


 
 
Decision variables

αs

g

Description 1 if supplier s g ∈ S g is selected, 0 otherwise

β pt

1 if plant p ∈ P with technology t ∈T is openned, 0 otherwise

γd

1 if distribution center d ∈ D is opened, 0 otherwise

λc

1 if cross-dock c ∈ C is opened, 0 otherwise

y rc (ξ )

1if retailer r ∈ R is assigned to cross-dock c ∈ C in scenario ξ ∈ Ξ,0 otherwise

x s g p (ξ )

flow of products from supplier s g ∈ S g to plant p ∈ P under scenario ξ ∈ Ξ

x pd (ξ )

flow of products from plant p ∈ P to distribution center d ∈ D under scenario ξ ∈ Ξ

x dc (ξ )

flow of products from distribution center d ∈ D to cross-dock c ∈ C under scenario ξ ∈ Ξ

x cr (ξ )

flow of products from cross-dock c ∈ C to retailer r ∈ R under scenario ξ ∈ Ξ

x dr (ξ )

flow of products from distribution center d ∈ D to retailer r ∈ R under scenario ξ ∈ Ξ

h pt (ξ )

amount of product manufactured at plant p ∈ P with technology t ∈T under scenario ξ ∈ Ξ


















Table 4. Factors and their levels



   
(5 × 4 × 3 × 3 × 2 × 4)











 
(6 × 5 × 4 × 4 × 3 × 5)  
(7 × 6 × 4 × 4 × 3 × 6)






(R ×C × D × P ×T × S )

 
(8 × 7 × 5 × 5 × 4 × 7)



 *
(9 × 8 × 5 × 6 × 4 × 8)






!  
"!#


 +
(9 × 8 × 6 × 7 × 5 × 9)

) '










$%&
"$#


 ,
(10 × 10 × 7 × 7 × 6 × 9)



 ' 
  
"#


 -
(11 × 12 × 8 × 8 × 7 × 10)




 ( 
"#


 /
(12 × 12 × 9 × 8 × 7 × 12)



 
"#



  (14 × 13 × 10 × 9 × 8 × 12)

.(






' 
"#









  (15 × 14 × 11 × 10 × 8 × 14)




  (18 × 15 × 12 × 8 × 12 × 16)




N (100,20)


 

  

" Demr (ξ ) #
N orm al distribution ( µ , σ )


' 0

s g


U (1,3)


1' 


(2


TECpt


U (100,400)





" VCc VCd VC s #


U (1,5)


g

T C cr TCdr


U (5,50 )

TCdc TC pd TCsg p


U (10,90)

 
'
  



ETcr ETdr


U ( 4,14)




ETdc ET pd ETsg p


U (10,20)




 
'
  












' '(



'
'

 EMpt


U ([10,40], P ,T )



' '(



'
'


MCpt

U ([1,5], P ,T )







*3
4( 
 
  
 
' 

5   
No. Experiments

Factors and their levels

Algorithm

total No . = ( n f , nax , ncp )

Lower limit(-), Upper limit(+)

NSGAII

MaxIt (-) , (+) (300) , (600)

Pc (-) , (+) (0.6,0.8)

Pm (-) , (+) (02,0.4)

20=(23,6,6)

MaxIt (-) , (+) (300,600) MaxIt

MOPSO

(-) , (+)

MOHEV

(300,600)

nPop (-) , (+) (40,100) nPop

nGrid (-) , (+) (5,15) LSiter

(-) , (+)

(-) , (+)

(40,100)

(8,12)

20=(23,6,6) 20=(23,6,6) Block1(cd)=(23,0,4) Block2(md)=(0,6,2)






Table 6. Tuned parameters, R-squared (R 2 ) , and Desirability (D) 


' 
 





488










!

!

!

!


 


6


74


D


  2






9+

9 ,/,



 9 





+:


-:


-+:


-:



 ,**


)66


9*-/



9





9 *





*-:


,*:


*:


,:



 ,,/


)6; 


9*/
 9/

9






+:


,:


+:


,*:



 -









,

  
 

( 
 
")
)8
6

74#
<










4( 
4( 
4( 
4( 

4( 
4( 








4( 
4( 

4( 
4( 
4( 
4( 

4( 
4( 
4( 
4( 




/
 +*+
-,+

--/
 ! //**
-+*,

 +






*,+
,
,-

 $ /
-,
-/



"#$  -
 -

 ,
 +-


"#

"#

"# "#  
"#
"%  +,
 +


"%



 % *,
/
,,+

,*+
**/


/,*
!
 *-
 
 +
"#!
 -


"$#  ,
 -




+
*,*
+
$#
*
-


! +,



"$  +-
 +-


"

 ,*
 +
 +

 ,/



*
*,,
$ +
-+,*

 +/
++
*-,

 /
 ,,
 
" $
"#  +,
 +
 +*
 &' 
+
 +*+
*+
-+,

 # //*-
,-,
-+-+


,
  -+/
--
,-*

% /+
//
/-,





"


-
 !$ */

+

!  /,/


*+
,



"!#  --
 *
 +

"  ,
 *-
 +



/
##! +
-+
/-*,

#  **/





*,,+



"!  //
 *+
 -

"  ,
 *,
 +




% +*+
+*
-,

 %#% *
,+,
-+--



"!#%  +*
 ,/
 ,

"$  ,*
 +
 ,




$ //**
+/
+-

!!! +/
**,
/*,



"!  +-
 *
 

"%  ,
 *
 */




  +
/-
,

$% ,-/

/



"  *,
 
 ,

"$  ,-
 **
 +


4( 3)6;7="#>
4(
3)6;7="#>
4( 3)66

4( 3 488





 ,
 ,
 *,

"$  ,+
 +
 ,


"%  *
 
 --

"%  ,
 +
 +/



-

 


 '2
" 
 ' 
  
%&#


Echelon Plants (p) (1) Tehran (IK central production unit) (Existing) (2) Mashhad (IK east production unit) (Existing) (3) Kermanshah (IK west production unit )

(4) Tehran Second unit (IK central production unit II ) (5) Yazd (IK south production unit )

Distribution center (d) (1) Tehran (Persian logistic Co.) (Existing) (2) Karaj (IK logistic Co.) (Existing) (3) Mashhad (4) Semnan (5) Isfehan (6) Zanjan

Cross-dock (c) (1) Kerman (2) Lorestan (3) Semnan

Parameter Production technology (t) Type 1 Type 2 Type 1 Type 2 Type 1 Type 2 Type 3 Type 1 Type 2 Type 3 Type 1 Type 2 Type 3

* p

* pt

(OC ,TEC , MC pt** , EO pt*** , EM pt**** ,Cap pt ) (0, 90, 20, 0, 1300, 52000) (0, 98, 18, 0, 900, 45000) (0, 76, 22, 0, 950, 65000) (0, 80, 17, 0, 1400, 57000) (245, 91, 24, 60, 1250, 68000) (165, 104, 19, 65, 1050, 50000) (188, 87, 18,52, 950, 56000) (156, 110, 18, 48, 1400, 64000) (194, 94, 19, 44, 1350, 56000) (173, 93, 23, 72, 980, 47000) (232, 85, 16, 35, 1080, 59000) (175, 79, 18, 42, 1050, 50500) (158, 95, 19, 50, 1150, 56500)

(OC d* ,VC d** , EOd*** ,Capd ) (74, 12, 18, 47500) (68, 14, 15, 55000) (66, 11, 14, 53500) (55, 9, 19, 48000) (49, 14, 14, 65000) (72, 13, 17, 62500)

(OC c* ,VC c** , EOc*** ,Capc ) (38, 4, 15.3, 38000) (37, 5, 14.8, 32000) (22, 4, 12.6, 35000)

*Thousand billion rials ** Million rials *** Million grams **** Thousand grams

Table 9. The data of case study (Suppliers) Evaluation Score Class (Out of 100) (1) Karaj (Mehvar Mashin Co.) 0.67 B (2) East of Tehran (Shetabkar Co.) 0.87 A (3) Rasht (Persian Khodro Co.) 0.45 C (4) West of Tehran (Mehvarsazan Co.) 0.62 B (5) West of Tehran (Mehr-kam pars Co.) 0.59 C (6) West of Tehran (Sapco Co.) 0.92 A (7) Tabriz (Pars Automobile part Co.) 0.82 A (8) Tabriz (Charkheshgar Co.) 0.62 B (9) Tehran (Crouz Co.) 0.74 B Suppliers (s)

* Thousand million rials ** Ten thousand grams

(FC s*g , ES s**g ,Cap s g ) (53, 21,17800) (74, 13, 14600) (78, 23, 18600) (52, 18, 16900) (54, 17, 15200) (56, 15, 14800) (71, 16, 13600) (72, 19, 17200) (59, 21, 15600)

Table 10.The data of case study (retailers) (×102 ) Scenarios ( ξ )

1

2

3

4

Scenario Probability ( P (ξ ) ) Retailers (r)

0.1

0.2

0.3

0.4

(1) Shiraz (Shiraz branch IK centre of sale) (2) Ilam (West Leasing Co.) (3) Hamehdan (Persian Leasing Co.) (4) Sari (Mazandaran IK centre of sale) (5) Karaj (IK centre of sale Alborz) (6) Center of Tehran (IK central sale unit.) (7) North of Tehran (IK sale centre Vanak) (8) Kermanshah (West Pars Leasing Co.) (9) Yazd (IK centre of sale Yazd) (10) Tabriz (Azarbayjan Leasing Co.) (11) Mashhad (IK centre of sale Khorasan) (12) Esfahan (Hamrah Leasing Co.) (13) East of Tehran (IK centre of sale Resalat)

Dem r (1)

Dem r (2)

Dem r (3)

Dem r (4)

~ Unit [50,63] ~ Unit [25,35] ~ Unit [13,25] ~ Unit [125,150] ~ Unit [115,138] ~ Unit [125,150] ~ Unit [70,80] ~ Unit [40,60] ~ Unit [105,120] ~ Unit [90,120] ~ Unit [30,40] ~ Unit [90,115] ~ Unit [135,155]

~ Unit [63,88] ~ Unit [35,45] ~ Unit [25,38] ~ Unit [100,130] ~ Unit [125,155] ~ Unit [115,140] ~ Unit [75,115] ~ Unit [55,75] ~ Unit [60,90] ~ Unit [115,135] ~ Unit [88,115] ~ Unit [115,140] ~ Unit [175,200]

~ Unit [80,120] ~ Unit [38,55] ~ Unit [30,63] ~ Unit [50,115] ~ Unit [150,180] ~ Unit [75,115] ~ Unit [90,145] ~ Unit [70,95] ~ Unit [90,110] ~ Unit [50,80] ~ Unit [60,100] ~ Unit [65,95] ~ Unit [150,180]

~ Unit [113,150] ~ Unit [50,80] ~ Unit [50,75] ~ Unit [75,100] ~ Unit [170,195] ~ Unit [50,100] ~ Unit [140,170] ~ Unit [30,45] ~ Unit [45,65] ~ Unit [65,95] ~ Unit [40,65] ~ Unit [50,70] ~ Unit [200,240]




'  
 ' 


 '2


Design

Total cost

Design1st

2805

Design4th

Design7th

Total Environmental

Type of design

Priority of suppliers in each class

Cost oriented

Class A=[ S6, S7, S2] Class B=[ S4, S8, S1, S9] Class C=[ S5, S3]

Cost oriented

Class A=[ S6, S7, S2] Class B=[ S4, S8, S9, S1] Class C=[ S5, S3]

Mediocre

Class A=[ S6, S2, S7] Class B=[ S8, S1, S4, S9] Class C=[ S5, S3]

9 (×1015 ) effect (× 10 )

3004

3254

Design10th 3318

Design13th 3609

1403

1282

1110

965

903

Environment oriented

Class A=[ S2, S6, S7] Class B=[ S1, S4, S9, S8] Class C=[ S3, S5]

Environment oriented

Class A=[ S2, S6, S7] Class B=[ S1, S4, S8, S9] Class C=[ S5, S3]

Optimal value of binary decision variables (first phase solution)

β pt = [1, 0, −;1,0, −;0, 0,0;1, 0,0;0, 0,1]

γ d = [0, 0,1,1,1, 0] λc = [0,1,1] β pt = [1, 0, −;1, 0, −; 0, 0, 0; 0, 0,1;0,1, 0]

γ d = [0, 0,1,1, 0,1] λc = [1, 0,1] β pt = [1, 0, −;0,1, −; 0,1, 0;0, 0,1;0, 0, 0]

γ d = [0,1, 0,1, 0,1] λc = [0,1, 0] β pt = [ 0,1, −; 0,1, −; 0, 0,1;0, 0, 0;0,1, 0]

γ d = [1, 0,1, 0, 0,1] λc = [1,0,0] β pt = [0,1, −;0,1, −;0, 0, 0;0,1, 0;1, 0, 0]

γ d = [1,1, 0, 0, 0,1] λc = [1,1, 0]

Figures

Figure 1. A schematic example of a supply network

Procedure MOHEV() Initialization () Set Parameters ( M axIt , PopSize

)

MaxIt : maximum number of iterations PopSize : number of charged particles Generate initial particles randomly Wh ile (hasn't met stop criterion (MaxIt )) do Calculate resultant force F () Move particle by F () S earch Locally () end while Figure 2. Pseudo-code of MOHEV main loop

1

     

 






 $








$




$









$










$













 






















 $


 


 












PMax = 3



































 !"# 






















  





 



































 



Figure 3. Schematic illustration of the structure of one sub-solution in the plant echelon with three parts: 1) five potential plants (p), 2) four selected distribution centers (d), and 3) four technology levels (t)

2

Initialization Procesdure() for i = 1to PopSize for n = 1to N λ ← U (0,1) x ni = l n + λ (u n − l n ) n = 1, 2,..., N end for find design by parsing method



f (x i ) := [f 1 (x i ),..., f k (x i )] A := Updating pareto archive set() end for Figure 4. Pseudo-code of initialization procedure Crowding distance procedure() L=S

/ /number of solutions in S

for each i , S [i ]distance = 0 / /initialize distance for each objective m S = sort (S , m ) / /sort using each objective function S [1]distance = S [L ]distance = ∞ / /boundarysolutions are alwaysselected for i = 2 to (L − 1) / /for all other points S [i ]distance = S [i ]distance + (S [i + 1].m − S [i − 1].m ) where S [i ]. m is m th objective function value of the i th solution in population S end for end for end for

Figure 5. Pseudo-code of CD procedure

Figure 6. Best particle selection in two strategies S1 (MD) and S2 (CD)

3


X1
X3
F 13


F 23


X2
F3 Figure 7. An example for the resultant of forces

Resultant force caculation procedure()
Select one of the two strategies S 1 or S 2 to determine x best if strategy == S 1 for i = 1 to Popsize




prox (x i ) = min( f ( x i ) − f (x ) )

x ∈A

end for


i x best = arg min

{ prox ( x )} // p = {pop1 ,pop 2 ,...,pop PopSize } x ∈p

else

x best =CD procedure() // x best ∈ A | CD x
best = max(CD i ) &CD i ≠ ∞ endif for i = 1 to PopSize § q i = exp ¨ −N ¨ ©
i F := 0 end for



· prox( x i ) − prox( x best ) ¸ m
j
best ¸ ¦ j =1 (prox(x ) − prox(x )) ¹

for i = 1 to PopSize for j = 1 to PopSize

if prox(x j ) < prox(x i ) then

qiq j

F i := F i + (x j − x i )

xi −x j

elseif prox(x j ) > prox(x i ) then

qiq j

F i := F i + (x i − x j )

xi −x j

2

//{Attraction}

2

//{Repulsion}

else



F i := F i // do nothing,

end if end for end for Figure 8. Pseudo-code of the Resultant Force procedure

4

Update Position procedure() for i = 1to PopSize

if x i ≠ x best then

Fi F i :=
i // normalizing the force vector F for n = 1to N do

λ ← U (0,1) if Fni > 0 then xni := xni + λ Fni (un − xni ) else xni := xni + λ Fni ( xni − ln ) end if end for end if



f ( x i ) := [ f1 ( x i ),..., f k ( x i )] // evaluation the new position A := Updating Pareto archive set () end for Figure 9. Pseudo-code of the Update Position procedure

5

Updating Pareto archive set ()
for each u ∈ P P : Initialized Population  A=φ nA = 0 n A :Number of members A for 1: PopSize




If x  u , ∃u ∈ A if x dominates u

A := A − {u} // Eliminate u from the set of non-dominated solutions nA = n A − nu



A := A ∪ x for all x  u n A = n A + 1 else do nothing to A endif end  if n A ≥ PopF // PopF :the capacity of front A Sort Pareto archive A based on CD A = A (1: PopF ) // Eliminate the excess solutions else do nothing to A end Figure 10. Pseudo-code of the Updating Pareto Archive Set procedure

6

Local search ( AMOVNS)procedure () recieve P and A // P = {1, 2,.., PopSize }, A = {non-dominated solutions} // initilalization P = [ ] / / Temporary archive LS

for LSitr = 1to Maxlsitr do Select strategy S 1 or S 2 based on Roulette Wheel Method () if S 1 is selected



s is chosen randomly among (x 1 ,.., x PopSize ) else
s is selected based on Crowding distance procedure() among x i ∈ A endif for N ns = 1 to ns max do determine a random solution s ′ from N ns //shaking perform a neighborhood search N ns on s ′ to find s ′′ Evaluate the solution s ′′ if strategy == S 1 P = P ∪ s ′′ LS

LS

else if s ′′  ε s Update A with s ′′ and break elseif s ′′  s / /no one dominates each other select s or s ′′ randomly and if s ′′ selected, Update A with s ′′ and break endif endif endfor A := Updating pareto archive set(PLS ) endfor Figure 11. Pseudo-code of the Local Search (AMOVNS) procedure

Figure 12. Neighborhood structures in the sub-solution selection process

7

Figure 13. NSs in the neighborhood structures in the transship rate determination process

Figure 14. Desirability Curves for Minimum Goal

8

x 10

MOHEV(cd)

MOHEV(md)

MOPSO

NSGAII

8.47

8.45

(a) Small scale 8.43

8.41

9.186

9.272

9.358

9.444

9.04

MOHEV(cd)

MOHEV(md)

MOPSO

NSGAII

9.025

(b) Medium scale

9.01

8.995

8.98 9.538

9.616

MOHEV(cd)

4.52

9.694

MOHEV(md)

9.772

MOPSO

9.85

NSGAII

4.465

(c) Large scale

4.41

4.355

5.526

5.652

5.778

5.904

6.03 5

Figure 15. Dispersion of non-dominated solutions obtained by different algorithms in three scales (×105 )

MOHEV(cd)

MOHEV(md)

MOPS O

MOHEV(cd)

NS GAII

MOHEV(md)

MOPS O

NS GAII

RPD for SNS

RPD for DM

0.8 0.7

0.4

0.45

0.1

0.1

Small MOHEV(cd)

Medium MOHEV(md)

Small

Large MOPS O

NS GAII

MOHEV(cd)

Medium MOHEV(md)

Large MOPS O

NS GAII

RPD for DEA

RPD for POD

1 0.7 0.4 0.1

.3 0.2 0.1 0

Small

Medium

Large

Small

Medium

Large

Figure 16. Plots of evaluation metric value for the interaction between the type of algorithm and size of problems 9

Figure 17. Dispersion of non-dominated solutions obtained by MOHEV(cd) for the case problem (×1015 )

10