Toward an integrated sustainable-resilient supply chain: A pharmaceutical case study

Transportation Research Part E 103 (2017) 109–142

Contents lists available at ScienceDirect

Transportation Research Part E journal homepage: www.elsevier.com/locate/tre

Toward an integrated sustainable-resilient supply chain: A pharmaceutical case study Behzad Zahiri a,⇑, Jun Zhuang b, Mehrdad Mohammadi c a

Department of Operations Management & Strategy, School of Management, State University of New York at Buffalo, Buffalo, NY, USA Department of Industrial and System Engineering, State University of New York at Buffalo, Buffalo, NY, USA c Ecole des Mines de Saint-Etienne, Department of Manufacturing Sciences and Logistics, CMP, CNRS UMR, 6158 LIMOS, 880 avenue de Mimet, 13541 Gardanne, France b

a r t i c l e

i n f o

Article history: Received 12 January 2017 Received in revised form 10 March 2017 Accepted 22 April 2017

Keywords: Sustainability Resiliency Uncertainty Pharmaceutical supply chain network design HIV LGBTQ

a b s t r a c t In this paper, a novel multi-objective integrated sustainable-resilient mixed integer linear programming model for designing a pharmaceutical supply chain network under uncertainty is presented. To cope with the uncertainty aspect of the model, a new fuzzy possibilistic-stochastic programming approach is developed. Additionally, due to NPhard nature of the problem, we propose a novel Pareto-based lower bound method as well as a new meta-heuristic algorithm. Several numerical examples, as well as a case study targeting TruvadaÓ supply chain for the LGBTQ community, as they account for majority of the market for such product, in France is proposed. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The pharmaceutical, as an immense global industry, is responsible for the manufacturing, development and marketing of medications. Based on the 2015 CMR Pharmaceutical R&D Factbook, global sales reached a milestone of $1 trillion in 2014, and will reach $1.4 trillion in 2020, driven by increased healthcare access in emerging markets and high-priced new drugs for cancer and other diseases (http:reuters.com). In general, a pharmaceutical supply chain (PSC) is a 5-tier supply chain which includes primary manufacturers, secondary manufacturers, main and local distribution centers (DCs), and destination zones/demand points (e.g., pharmacies, hospitals, clinics, etc.). The primary manufacturers are in charge of the production of required active ingredients (RAI), which normally include either several chemical synthesis and separation stages to build up the involved complex molecules, or purification and product recovery in case of biochemical processes (Shah, 2004). Secondary manufacturers are responsible for further production processes with different technology levels, packaging and finalizing the products that are usually in SKU form. Compared to a typical manufacturing supply chain, primary manufacturers can be referred to as suppliers of raw materials, and secondary manufacturers as manufacturing centers. Consequently, secondary manufacturers do not only play a major role in the production of final goods, but they can also store up a limited amount of products within the facility. Both main

⇑ Corresponding author. E-mail address: [email protected] (B. Zahiri). http://dx.doi.org/10.1016/j.tre.2017.04.009 1366-5545/Ó 2017 Elsevier Ltd. All rights reserved.

110

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

and local DCs are in charge of stocking goods to satisfy the market demand. Compared to main DCs, local ones are often more dispersed to cover more demand points, and more limited in terms of capacity levels. Despite coupling with sophisticated technologies and improvements in quantity and quality of associated products in PSCs, many companies are far from effectively satisfying market demands with respect to arisen concerns (Masoumi et al., 2012). With ever-increasing awareness of sustainability among government, industry, and the general public over the last two decades, policymakers worldwide have sought to incorporate sustainability considerations into urban and industrial development, and pharmaceutical industry has not been an exception (Fiksel, 2006, 2003). Sustainability is mainly referred to as a balance between social, economic, and environmental issues involved in human development. Social aspect, which has rarely been addressed in the literature, corresponds to forcing non-governmental organizations to take responsibilities for the social impacts of their actions. One of the main aspects of social responsibility is community involvement and development since there is a consensus that the communities around the workplace should be enhanced both socially and economically. More precisely, increasing employment opportunities as well as providing a balanced economic development for local communities are the main goals of social responsibility (Zhalechian et al., 2016). Based on the World Summit of Sustainable Development, environmental impacts are considered an important pillar of sustainability. There is an ever-increasing consensus that environmental impacts (EIs) left uncontrolled, may result in major changes to both the climate and the environmental systems. Governments are under immense pressure to enact legislation to curb these impacts. These restrictions, which include controlling greenhouse gas emissions, more specifically carbon dioxide (CO2), is becoming a growing interest and companies are being urged to incorporate these issues into supply chain management schemes (Pagell and Wu, 2009; Benjaafar et al., 2013; Ding et al., 2016). Among these mechanisms and legislations, carbon emission trading or cap-and-trade system is generally accepted as one of the most effective market-based mechanisms, which has recently been broadly adopted by the EU, UN, etc. More precisely, the cap-and-trade approach provides economic incentives to decrease the emissions of pollutants, and is considered as market-based government-mandated mechanism. To respond to this issue, and in addition to cap-and-trade approach, companies tend to adopt more energy efficient technologies, equipment, or vehicles. On the other hand, these goals can be approached with optimizing production decisions, facility locations, transportation and inventory (Hua et al., 2011; Stavins, 2003). Upon doing so, an optimized supply chain network design (SCND) toward sustainability can not only correspond to economically beneficiary, but also lead to a significant change in social and environmental impacts. However, with all increases in complexities in supply chains as in PSCs, vulnerabilities to disruptions have significantly been increased in such a way that any temporary stoppage out of disruption source either external (e.g., natural disasters, fire, theft, etc.) or internal (e.g., labor strike, technology malfunctions, etc.) leads to a significant loss. In order to improve the ability to respond rapidly and cost effectively to such unpredictable changes, SCs must adopt new strategies and incorporate the concept of resiliency in the designing process of the SC (Cardoso et al., 2015; Craighead et al., 2007). Supply chain resilience is defined as the capacity for a SC to survive and retain its basic function in the face of any turbulent change (Pettit et al., 2010). Despite the increasing efforts on design and management of a sustainable SC, the quantitative impact of sustainability interventions on the supply chain resiliency has remained unexplored. Following sustainability in an environment characterized by frequent inevitable disruptions necessitates modeling and analysis that can accommodate this dynamism and complexity. Additionally, static sustainability analysis is simplistic, since both the economic and non-economic sustainability performance of a SC can be affected by disruptive events. Accordingly, this requires further optimization techniques and management tools to develop resilient and sustainable SCs, in such a way that the effect of any disruptions on sustainability performance is minimized (Fahimnia and Jabbarzadeh, 2016). Given these two mainstream topics that were discussed (i.e., sustainability and resiliency), we aim to investigate how PSC can be coupled with these two goals. It should be mentioned that due to high volatility and dynamic nature of parameters, most especially supply and demand in such an industry, and knowing the fact that such change can greatly affect the structure of network designs, they are considered under uncertainty. In general, the techniques to cope with uncertainty can be divided into three main categories; Fuzzy programming, Robust optimization, and Stochastic programming. More details on how we coped with such uncertainty are provided in Section 4. With respect to the abovementioned issues, this paper develops a new integrated sustainable-resilient pharmaceutical supply chain network design (ISRPSCND) in a multi-period planning horizon under uncertainty. The considered objective functions try to minimize total network cost, maximize total sustainability in terms of minimizing environmental impact and maximizing social impacts, and minimize de-resiliency based on the developed measures. The main contributions of this paper that differentiate it from the existing literature are listed as follows:
Proposing a new multi-period pharmaceutical network design including both strategic decisions (e.g., location decisions) and tactical ones (i.e., network allocation).
Simultaneously considering both resiliency and sustainability as the objective functions and setting a trade-off between these two measures.
Both introducing and developing resilience measures for the considered network and applying them as optimization tools.

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

111


Developing a new technique to cope with the inherent uncertainty of the input parameters based on mean-absolute deviation and fuzzy possibilistic-stochastic programming.
Developing a new hybrid meta-heuristic algorithm to efficiently solve the large-sized instances.
Proposing a lower bound to evaluate the performance of the developed meta-heuristic algorithm.
Applying the presented model in a case study. The rest of paper is organized as follows. In Section 2, a comprehensive literature review is provided. Section 3 is dedicated to problem definition and mathematical formulation. The proposed solution approach is discussed in Section 4. In Section 5, the considered case study and some sensitivity analyses are provided. Finally, Section 6 discusses the conclusion and some new avenues for future research. 2. Literature review The literature in this field can be categorized into three main fields namely (i) PSC network design (ii) sustainable network design, and (iii) resilient supply chain. Table 1 summarizes the studied papers in these fields. Despite wide surveys on different challenges and managerial issues on pharmaceutical SCs (see Papageorgiou, 2009; Laínez et al., 2012; Yu et al., 2010; Narayana et al., 2012; Alnaji and Ridha, 2013; Rossetti et al., 2011), network optimization research remains far from adequately addressing industry challenges. In this domain, all efforts have been mainly devoted to optimizing certain parts of the chain, while only few articles have addressed the whole chain. In case of sustainable supply chains, research has tended to target empirical and conceptual studies with only a scant, but increasing number of articles published on quantitative analysis and analytical modeling (Fahimnia and Jabbarzadeh, 2016). Fiksel (2003, 2006) have investigated on how an industrial enterprise can work toward sustainability, and further resiliency by adopting a fresh perspective based on systems thinking. From the analytical perspective, extensive reviews have been done by Eskandarpour et al. (2015) and Seuring (2013). These reviews study published analytical papers in the field of sustainable supply chain network design, and is based on four aspects; (i) the included environmental and social objectives, (ii) the integration methods in the provided models, (iii) applied/developed methods, and (iv) the industrial applications and contexts. The reviews find some limitations to the current research in this topic; (1) limited and narrow-scope inclusion of environmental and social measures, and (2) scant inclusion of uncertainty and improved multi-objective approaches. While many studies have investigated sustainable SCND, the literature is scant on the sustainable network design of the PSC. The findings discovered from the research papers reveal the narrow scope of sustainability and resiliency in SC, the lack of uncertainty consideration, superficial utilization of sustainability and resiliency measures, and insufficient address of the tradeoff between sustainability and resiliency. Moreover, neither of these concerns has been considered on PSC. To address this gap, we propose a novel ISRPSCND introducing new measure of resiliency and considering issues like cap-and trade system in the network under uncertainty. 3. Problem framework Fig. 1 shows the schematic view of the considered PSC network. Since the RAIs are imported through different channels and countries and the corresponding strategy may vary based upon the importer’s policies, the studied chain starts from the secondary manufacturers to the demand points. As discussed earlier, RAIs after production are sent to secondary manufacturers (from now on manufacturers) for further production process with different technology levels, packaging and finalizing the products. Since there are different capacities of such technologies, managing their investment and their compatibility to social impacts are major concerns in the production aspect of such industry. To cope with this issue, both the purchasing and the selling of different technologies are allowable over time. Additionally, in case of any disruption resulted by external sources (e.g., natural disasters, etc.) or internal ones (e.g., unreliability of technologies, etc.), these technologies may become unavailable/fail. In such situations, backup technologies are needed to meet the faced demand. Therefore, at the production level, the productions of products are allocated to technologies level by level (i.e., backup technologies). By performing it like this, any case of any failure in a certain technology will be avoided due to the fact that production will be reallocated to the backup technology in a consequent level. This reallocation policy is continued until the production is allocated to a nonfailable technology. The final products are then shipped to DCs in order to be stocked and further shipped to demand points. In case of any demand, local DCs are responsible to meet the demand of regional markets. The presented model integrating both strategic decisions (i.e., optimizing location of facilities - manufacturing centers, main and local DCs) and tactical ones (i.e., production, purchasing and selling of set of technologies, decisions on carbon credits, inventory management and capacity planning, allocation and transportation, etc.) in addition to all abovementioned issues, focuses on two main goals in network design of a PSC. (i) sustainability, which itself includes three subclasses of (1) network’s total cost, (2) environmental impact in terms of omitted CO2, and (3) social impacts of business activities, and (ii) resiliency. For environmental impact, the main focus is made toward production and shipping of goods, in such a way that for production, the utilization of different technologies with respect to the buying and selling of carbon credits may end up having certain CO2 emission. For social impacts, establishment of new facilities of any kind (either manufacturing centers, or any

112

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Table 1 Reviewed papers. Reference

Sustainability Measures

Resiliency Measures

Planning Horizon

Type of Uncertainty

Time Period

Solution Method

Case Study/ Application

Papageorgiou et al. (2001) Levis and Papageorgiou (2004) Sousa et al. (2005) Susarla and Karimi (2012) Sousa et al. (2011) Grunow et al. (2003) Hansen and Grunow (2015) Mousazadeh et al. (2015) Weraikat et al. (2016) Varsei and Polyakovskiy (2017) Zhalechian et al. (2016) Pop et al. (2015) Balaman and Selim (2016) Weraikat et al. (2016) Klibi and Martel (2012) Carvalho et al. (2012) Torabi et al. (2015) Cardoso et al. (2015)

– – – – – – C C – C, GE, ED C, GE, JO, ED C, GE C C C C C C

Str, Tc Str, Tc Str, Tc Tc Tc Tc Tc St, Tc Tc St, Tc St, Tc St, Tc St, Tc Tc St, Tc Tc Tc Tc

– – – – – – S R – – F – F – S – S S

M M M M M M M M Sg Sg M Sg M Sg M M Sg M

Ex He He Ex He He Ex Ex Ex Ex He He Ex He Ex Sm He Ex

Ph Ph Ph Ph Ph Ph Ph Ph Ph FB E G G Ph G Am G G

Fahimnia and Jabbarzadeh (2016) Levalle and Nof (2015a) Levalle and Nof (2015b) Falasca et al. (2008) Our paper

C, GE, RC, ED, PR, LR – – – C, GE, CT, JO, ED

– – – – – – – CDL – – – – CDL – MS, CDL LT SBD NCr, FC, NC, D, CF, CD, ENPV, CDL SBD TC, TQ, SBD TC, TQ D, FC, NC, NCr NCr, MS, FC, NC, CDL

St Tc Tc – St, Tc

S – – – F

Sg Sg Sg – M

Ex Sm Sm Sm He

Cl G G G Ph

Am: Automotive. C: Cost. CD: Coutnd. CDL: Customer de-service Level/Coverage. CF: Coutflow. Cl: Clothing. CT: Cap and Trade. D: Density. E: Electronics. ED: Economic development. ENPV: Expected net present value. Ex: Exact. F: Fuzzy. FB: Food and beverage. FC: Flow complexity. G: General. GE: GHG emission/Environmental performance score. He: Heuristic/meta-heuristic. JO: Job opportunities. LR: Labor and human rights. LT: Lead time. M: Multi period. MS: Multiple sourcing/Reassignment. NC: Node complexity. NCr: Node criticality. Op: Operational. Ph: Pharmaceutical and medication. PR: Product responsibility. R: Robust optimization. RC: Resource consumption. S: Stochastic programming. SBD: Scenario-based disruption. Sg: Single period. Sm: Simulation. St: Strategic. Tc: Tactical. TCF: Total cost of flow. TQS: Total quality of service.

kind of DCs) are not only taking place with respect to potential created job opportunities and unemployment rate, but also the regional development rate. This will give more focus toward the less developed areas. For resiliency, even though the main focus has been made towards the production aspect of the chain, different measures of resilience have been developed and inserted into the proposed model as optimization tools. This makes the whole network

113

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Fig. 1. The schematic view of a typical pharmaceutical supply chain.

more reliable toward any kind of disruptions/disasters, and guaranteed to meet the faced demand in any case. More details about these measures are provided in Section 3.1. The following notations are introduced for the proposed mathematical model: Indices i j k h l r, s p m t

index of candidate locations for manufacturers i ¼ f1; . . . ; Ig index of candidate locations for main DCs j ¼ f1; . . . ; Jg index of candidate locations for local DCs k ¼ f1; . . . ; Kg index of demand zones h ¼ f1; . . . ; Hg index of technology levels l ¼ f1; . . . ; L þ 1g, in such way that l ¼ f1; . . . ; Og indicates the old technology levels and l ¼ fO þ 1; . . . ; L þ 1g ¼ N denotes the new technology levels index of allocation levels r; s ¼ f1; . . . ; Rg index of product types p ¼ f1; . . . ; Pg index of transportation modes m ¼ f1; . . . ; Mg index of time periods t ¼ f1; . . . ; Tg

Parameters ~f fixed establishment cost of a manufacturing center in candidate location i i ~f 0 fixed establishment cost of a main DC in candidate location j j ~f 00 fixed establishment cost of a local DC in candidate location k k f cq ilp unit production cost of product type p using technology l in manufacturing center i ~cijm unit transportation cost between nodes i and j using transportation mode m ~c0ikm unit transportation cost between nodes i and k using transportation mode m ~c00jkm unit transportation cost between nodes j and k using transportation mode m ~c000 khm f ch ip

f0 ch jp f ch 00

kp

cf vi cf v 0i f cp il f cp 0il ~j i ~j0 j ~j00 k ur i ur 0j ur 00k ef vi ef v 0j f ev 00k rdi 0 rdj 00 rdk t ~r il r 0tilp ð:Þ e ei ijm e ei 0ikm e 00 ei jkm

unit transportation cost between nodes k and h using transportation mode m unit holding cost of product type p in manufacturing center i unit holding cost of product type p in main DC j unit holding cost of product type p in local DC k unit selling price of carbon credit unit purchasing price of carbon credit cost of purchasing new technology level l in manufacturing center i cost of selling new/old technology level l in manufacturing center i number of created job opportunities out of opening manufacturing center i number of created job opportunities out of opening main DC j number of created job opportunities out of opening local DC k unemployment rate in region i unemployment rate in region j unemployment rate in region k regional economic value of manufacturing center i regional economic value of main DC j regional economic value of local DC k regional development factor in region i (0 6 rdi 6 1) 0 regional development factor in region j (0 6 rdj 6 1) 00 regional development factor in region k (0 6 rdk 6 1) fixed setup emission of technology level l in manufacturing center i in period t unit production emission function for product type p with technology level l in manufacturing center i in period t environmental impact of shipping one unit product from node i to j using transportation mode m environmental impact of shipping one unit product from node i to k using transportation mode m environmental impact of shipping one unit product from node j to k using transportation mode m (continued on next page)

114

e 000 ei khm ~ n 0 ne f n00

a~ ae0 af00 af000 ~ b be0

f b00

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

environmental impact of shipping one unit product from node k to h using transportation mode m penalty coefficient for critical manufacturing centers penalty coefficient for critical main DCs penalty penalty penalty penalty penalty penalty penalty

coefficient coefficient coefficient coefficient coefficient coefficient coefficient

for for for for for for for

critical local DCs flow complexity between nodes i and j flow complexity between nodes i and k flow complexity between nodes j and k flow complexity between nodes k and h node complexity of manufacturing centers node complexity of main DCs

gilp

penalty coefficient for node complexity of local DCs penalty coefficient for demand dissatisfaction level penalty coefficient for production with old technologies maximum capacity level of manufacturing center i maximum storage capacity level of main DC j maximum storage capacity level of local DC k maximum production capacity of manufacturing center i for product type p using technology l

et D hp

demand of product type p in demand zone h in period t

~ tl u atil

probability of incident for technology l in period t disruption probability of technology level l in manufacturing center i in period t

Gti Smax t Sþmax t M

maximum allowable emission for manufacturing center i in period t maximum amount of carbon credit allowed to be sold at period t maximum amount of carbon credit allowed to be bought at period t A reasonably large number

c~ s~ Capi Cap0j Cap00k

Variables wi 1; yj 1; zk 1; v tij 1;

v 0tik v 00tjk v 000t kh

if if if if

a manufacturing center is established in candidate location i, and 0; otherwise a main DC is established in candidate location j, and 0; otherwise a local DC is established in candidate location k, and 0; otherwise node i is assigned to node j in period t, and 0; otherwise

1; if node i is assigned to node k in period t, and 0; otherwise 1; if node j is assigned to node k in period t, and 0; otherwise

etil util Itip

1; if node k is assigned to node h in period t, and 0; otherwise 1; if product type p is produced using technology l at level r in manufacturing center i in period t, and 0; otherwise 1; if new technology l 2 N is purchased in period t in manufacturing center i, and 0; otherwise 1; if a technology is sold in period t in manufacturing center i, and 0; otherwise inventory level of product type p in manufacturing center i in period t

I0tjp I00t kp qtilp xtijpm

inventory level of product type p in main DC j in period t inventory level of product type p in local DC k in period t quantity of produced product type p with technology level l in period t in manufacturing center i flow of product type p transferred from node i to j using transportation mode m in period t

x0tikpm

flow of product type p transferred from node i to k using transportation mode m in period t

x00t jkpm

flow of product type p transferred from node j to k using transportation mode m in period t

x000t khpm st i sþt i

flow of product type p transferred from node k to h using transportation mode m in period t

mrt ilp

quantity of carbon units to sell in period t quantity of carbon units to buy in period t

3.1. Measures of resilience Even though measuring resiliency of a network is still debatable, the existing literature provides some indicators, according to which such issue can be addressed. In this section, five measures according to the structure of the network along with some newly introduced ones are updated. 3.1.1. Node criticality According to this measure, a node is considered to be critical if the total inflows and outflows exceed a certain threshold. Eqs. (1)–(5) show the critical nodes for manufacturing centers, main DCs, and local DCs, respectively.

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

w0i ¼ 1j

XX X XX X XX qtilp þ xtijpm þ x0tikpm > U ti ; p2P l2L

y0j ¼ 1j

X XX

m2M p2P j2J

m2M p2P i2I

z0k ¼ 1j

ð1Þ

m2M p2P k2K

X XX 0t x00t jkpm > U j ;

xtijpm þ

8i; t;

115

8j; t;

ð2Þ

m2M p2P k2K

X XX X XX X XX 00t x0tikpm þ x00t x000t jkpm þ khpm > U k ;

m2M p2P i2I

m2M p2P j2J

000t qtilp ; xtijpm ; x0tikpm ; x00t jkpm ; xkhpm P 0;

w0i ; y0j ; z0k 2 f0; 1g;

8k; t;

ð3Þ

m2M p2P h2H

8i; j; k; h; l; p; m; t;

ð4Þ

8i; j; k:

ð5Þ

3.1.2. New technology and the (L + 1)th level/multiple assignment As discussed before, manufacturing process, most especially in pharmaceutical industry, are done using different sets of technologies (i.e., old and new), which even though purchasing/utilization of new ones are more costly, they are more reliable. On the other hand, in case of any disruption resulted by external or internal issues, utilization of some technology sets may become unavailable. In this situation, backup technology levels are needed to meet the demand. Therefore, in the setting of production, products are assigned to different technologies level by level, which means that if the corresponding technology level is unavailable/disrupted for any reason, the production taking place will use the backup technology level in a consequent level. This reassignment policy is continued until the production is assigned and the possibility of the technology level failing is avoided entirely (i.e., level ðL þ 1Þ). In doing so, the technology levels are categorized into three main levels: (i) old technology with lowest cost and highest disruption probability, (ii) new technology with medium costs and medium disruption probability, and (iii) the non-failable level with the highest cost and zero failure. 3.1.3. Flow complexity This measure calculates the total interaction between nodes, where the first term corresponds the assignments of main DCs to manufacturing centers, the second term shows the interaction between main DCs and local ones, and the last term refers to allocation of demand zones to local DCs. According to this measure, total flow in a network is considered to be complex if the total number of corresponding links is high. Eqs. (6) and (7) calculate total number of links in the network.

X XX

v tij þ

t

j2J i2I

XX

XX

XX

k2K i2I

k2K j2J

h2H k2K

v 0tik þ

v 00tjk þ

v 000t kh

!

;

v tij ; v 0tik ; v 00tjk ; v 000t 8i; j; k; h; t: kh 2 f0; 1g;

ð6Þ ð7Þ

3.1.4. Node complexity Similar to flow complexity, a network is called to be node complex, if the total number of active nodes is high. Eqs. (8) and (9) show the total number of established facilities in the network.

X X X wi þ yj þ zk ; i2I

j2J

ð8Þ

k2K

wi ; yj ; zk 2 f0; 1g;

8i; j; k:

ð9Þ

3.1.5. Customer de-service level We calculate customer de-service level as the unmet demand over each product type, demand zone and time period with a weight factor wthp (see Eq. (10)). More precisely, Eq. (10) tries to build a balance between unmet demand fractions by minimizing the maximum unmet demand ratio (i.e., equals to wthp  1) over different product types, time periods and demand zones.

min max h;p;t

wthp

et  D hp

X X

!

x000t khpm

;

ð10Þ

k2K m2M

Conversion of Eq. (10) into its linear counterpart is achieved using Eqs. (11)–(13);

min U et  U P wthp D hp

ð11Þ XX

! x000t khpm ;

ð12Þ

k2K m2M

U P 0:

ð13Þ

116

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

3.2. Mathematical formulation This section proposes a new mathematical model for designing an ISRPSCN. In doing so, the considered assumptions are listed as follows:
Location of demand zones are known, while only candidate location of manufacturing centers, main and local DCs are predetermined.
The tactical decisions are considered to be dynamic, in such a way that a certain decision may vary from period to period.
All facilities (i.e., manufacturing centers, main and local DCs) are capacitated.
There is a reliable technology production level (i.e., (L þ 1)th level), called non-failable technology level, that all production of products should be assigned to at a specific level.
Intermodal freight transport is allowable. In other word, different transportation modes are available to ship the products from different nodes to the end demand point.
Establishment of manufacturing centers are accompanied with built-in primary technology levels.
Because of the high volatility of demand and supply in such industry, the relevant parameters of supply and demand are considered to be uncertain (see Section 4).

min Z1 ¼

X X X XX X X X X XXX ~f w þ ~f 0 y þ ~f 00 z þ ~ ilp qtilp þ ~cijm xtijpm cq i i j j k k i2I

j2J

t2T p2P l2Lþ1 i2I

k2K

t2T m2M p2P j2J i2I

X X XXX X X XXX X X XXX XXX 000t ~ It ~c0ikm x0tikpm þ ~c00jkm x00t ~c000 ch þ ip ip jkpm þ khm xkhpm þ t2T m2M p2P k2K i2I

t2T m2M p2P k2K j2J

t2T m2M p2P h2H k2K

XXX XXX XX XX XXX ~ 0 I0t þ ~ 00 I00t  ~ il etil ch ch c~v i st c~v 0i sþt cp þ jp jp kp kp i þ i þ t2T p2P j2J

XXX ~ 0il util  cp

t2T p2P k2K

t2T i2I

t2T i2I

t2T p2P i2I

t2T l2Lþ1 i2I

ð14Þ

t2T l2Lþ1 i2I

0

"

X X X ~ji ur i wi þ ~j0 ur 0 y þ ~j00 ur 00 zk j j j k k

#max

1

X X X  ~ji ur i wi þ ~j0j ur 0j yj þ ~j00k ur 00k zk

B C B C B C i2I j2J i2I j2J k2K k2K C min Z2 ¼ xB " # " # B max min C X X X X X B X C 0 0 00 00 @ A ~ji ur i wi þ ~j0 ur 0 y þ ~j00 ur 00 zk ~ ~ ~ ji ur i wi þ jj ur j yj þ jk ur k zk  j j j k k i2I

j2J

0

"

k2K

i2I

j2J

k2K

#max X X X X X X 0 00 0 00 e~v i wi ð1  rdi Þ þ e~v 0j yj ð1  rdj Þ þ e~v 00k zk ð1  rdk Þ e~v i wi ð1  rdi Þ þ e~v 0j yj ð1  rdj Þ þ e~v 00k zk ð1  rdk Þ 

1

B C B C B C i2I j2J i2I j2J k2K k2K 0B þ x B" #max " #min C C X X X X X B X C 0 00 0 00 @ A e~v i wi ð1  rdi Þ þ e~v 0 y ð1  rd Þ þ e~v 00 zk ð1  rd Þ e~v i wi ð1  rdi Þ þ e~v 0 y ð1  rd Þ þ e~v 00 zk ð1  rd Þ  j j

i2I

j

j j

k

k

j2J

i2I

k2K

j

k

j2J

k

k2K

ð15Þ

min Z3 ¼

XXX XX X X X X XXX X X XXX ~ ijm xt þ ~ 0 x0t ~r til wi þ r 0tilp ðqtilp Þqtilp þ ei ei ijpm ikm ikpm t2T l2Lþ1 i2I

t2T p2P l2Lþ1 i2I

t2T m2M p2P j2J i2I

t2T m2M p2P k2K i2I

X X XXX X X XXX ~ 00 x00t þ ~ 000 x000t þ ei ei jkm jkpm khm khpm t2T m2M p2P k2K j2J

min Z4 ¼

ð16Þ

t2T m2M p2P h2H k2K

X

X 0 X f ~nw0 þ ne y0j þ n00 z0k i j2J k2K |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} i2I

þ

X

Nodecriticality

X X j2J

i2I

a~ v tij þ

X

X

X

X

t2T

k2K

i2I

X X f ~ iþ bw þ be0 yj þ b00 zk þ i2I j2J k2K |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Node complexity

ae0 v 0tik þ

X

X

X

t2T

k2K

Flow complexity

c~U |{z} Demand dissatisfaction level

þ

j2J

af00 v 00tjk þ

X

X

X

af000 v 000t kh k2K |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} t2T

X

X

t2T

X

h2H

X

s~qtilp i2I |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} t2T

p2P

l2O

ð17Þ

Production with old technologies

Objective function (14) tries to minimize total cost in the network consisting establishment cost, transportation cost, inventory holding cost, purchasing and selling carbon credits, and buying and selling technologies. Objective function (15) is referred to the social impact of opening facilities (i.e., manufacturing centers, main and local DCs) in the network. The first weighted term tries to maximize total created job opportunities out of opening centers with respect to the unemployment rate, and the second term maximizes the balanced economic development. Since the function is in the form of f max f , f max f min

minimization of this ratio would result in achieving the maximum value for f . It should be mentioned that both

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

f

max

and f

min

117

are parameters and in order to obtain their respective values, the objective function (i.e. f ) should be separately max

min

optimized under maximization (to achieve f ) and minimization (to achieve f ). Objective function (16) tries to minimize total environmental impacts in the network. The unit production emission of CO2 is considered to be a function of both quantity of products and time periods. Further details are provided in Section 3.3. Objective function (17) aims to minimize nonresiliency of the network with respect to the discussed measures in Section 3.1.1.

X XX XX qtilp  xtijpm  x0tikpm ;

t1 Itip ¼ Iip þ

m2M j2J

l 0t1 þ I0tjp ¼ Ijp

XX XX xtijpm  x00t jkpm ; m2M i2I

00t1 I00t þ kp ¼ Ikp

8i; p; t;

ð18Þ

m2M k2K

8j; p; t;

ð19Þ

m2M k2K

XX XX XX x0tikpm þ x00t x000t jkpm  khpm ; m2M i2I

m2M j2J

8k; p; t;

ð20Þ

m2M h2H

Constraints (18)–(20) are inventory balance constraints, for manufacturing center, main DC, and local DC, respectively.

X t1 X ðIip þ qtilp Þ 6 Capi wi ; p2P

8i; t;

X t1 X X ðIjp þ xtijpm Þ 6 Cap0j yj ; p2P

ð21Þ

l

8j; t;

ð22Þ

m2M i2I

XX X 00t1 X X 00 ðIkp þ x0tikpm þ x00t jkpm Þ 6 Capk zk ; p2P

m2M i2I

8k; t;

ð23Þ

m2M j2J

Constraints (21)–(23) are capacity constraints for manufacturing center, main DC, and local DC, respectively, and guarantee that total inflow in each period does not exceed the predefined maximum capacity.

XX

xtijpm 6 M v tij ;

8i; j; t;

ð24Þ

x0tikpm 6 M v 0tik ;

8i; k; t;

ð25Þ

00t x00t jkpm 6 M v jk ;

8j; k; t;

ð26Þ

p2P m2M

XX p2P m2M

XX p2P m2M

XX

000t x000t khpm 6 M v kh ;

8k; h; t;

ð27Þ

p2P m2M

Constraints (24)–(27) are assignment constraints. They guarantee that the link between each pair of nodes is only established when a flow exists in the associated link.

8i; l 2 fL þ 1g; p; r ¼ 1; t;

Prtilp ¼ 1  atil ;

0

Prtil0 p ¼ ð1  atil0 Þ

l X l

atil Pr1;t mr1;t ilp ; 1  atil ilp

8i; l 2 fL þ 1g; p; r > 1; t;

ð28Þ

ð29Þ

Constraints (28) and (29) are the transitional probability equations. For r ¼ 1, the probability of production of product type p with technology level l in manufacturing center i equals to technology l in manufacturing center i in period t is not disrupted. For 1 < r 6 R, the corresponding probability equals to atil ð1  atil0 Þ=ð1  atil ÞP r1;t when product type p is produced using techilp nology l in manufacturing center i in period t.

X rt Pilp qtilp 6 gilp ;

8i; l 2 fL þ 1g; p; t;

ð30Þ

r2R r X X mrtilp þ mstiðLþ1Þp ¼ 1; l2L

X mrtilp ¼ 1; r2R

8i; r; p; t;

ð31Þ

s

8i; p; l 2 fL þ 1g; t;

ð32Þ

118

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Constraint (30) is the maximum production capacity of manufacturing centers. Constraint (31) guarantees that for the production of a certain product type p in manufacturing center i in period t, either p is assigned to a technology type at level P r or is assigned to the non-failable ðL þ 1Þ at a certain level s 6 r (considering rs mst iðLþ1Þp ¼ 0 if r ¼ 1). Constraint (32) ensures that each technology is only assigned to one level.

qtilp 6 M

X mrtilp ;

8i; l 2 fL þ 1g; p; t;

ð33Þ

r2R

8i; l 2 fL þ 1g; p; r; t;

mrtilp 6 wi ;

ð34Þ

Constraint (33) links the production quantity to assignment variable, in such a way that in case of any production of product type p with technology l in manufacturing center i in period t, the respective assignment variable takes value 1. Constraint (34) assures that the production can only take place if the corresponding manufacturing center is open. 0

etil 6 1  etil ;

8i; l 2 N; t; t0 < t

ð35Þ

8i; l 2 fL þ 1g; t; t 0 < t

ð36Þ

0

util 6 1  util ; T t X X 0 0 util 6 etil ; t 0 ¼t

8i; l 2 N;

ð37Þ

t 0 ¼1 0

etil 6 1  util ; mrtilp 6

t X 0 etil ;

8i; l 2 N; t; t 0 < t

ð38Þ

8i; l 2 N; r; t;

ð39Þ

t0 ¼1 0

8i; 2 fL þ 1g; r; t; t 0 < t;

mrtilp 6 1  util ; etil 6

T X X mrtilp ;

ð40Þ

8i; l 2 N; t;

ð41Þ

t 0 ¼t r2R

8i; l 2 fL þ 1g; t;

util 6 wi ;

ð42Þ

Constraint (35) guarantees that in case of purchasing a new technology, purchasing the same technology in subsequent periods is not allowable. Similarly, constraint (36) ensures that in case of selling a technology, the same technology cannot be sold in subsequent periods. Constraints (37) and (38) assure that a technology can only be sold if it has been purchased in the previous periods, and after purchasing, it cannot be sold again. Constraints (39)–(41) link the assignment production variable to the purchasing and selling of technology variables. Constraint (42) guarantees that the selling of technologies is only possible if the corresponding manufacturing center is open.

X X X X þt ½wi~r til þ r 0tilp ðqtilp Þqtilp  Gti wi  þ ½wi~r til þ r0tilp ðqtilp Þqtilp  Gti wi  þ st i  si 6 0; p2P

l2N

8i; t;

ð43Þ

p2P

l2O

max st ; i 6 St

8i; t;

ð44Þ

þmax ; sþt i 6 St

8i; t;

ð45Þ

þt st i si ¼ 0;

8i; t;

ð46Þ

Constraint (43) enforces that the amount of carbon emission at any period, considering both buying and selling credits, incurred by any technology type to be lower or equal than the maximum allowable environmental impact (i.e., Gti ). Constraints (44) and (45) restrict the maximum carbon credit to be sold and bought in each period, respectively. Eq. (46) implies and sþt cannot take positive values concurrently. that st i i

U ti

XX XX X X XX  qtilp þ xtijpm þ x0tikpm p2P l2L

U 0tj 

X XX m2M p2P i2I

j2J p2P m2M

xtijpm þ

ð1  w0i Þ P eð1  w0i Þ;

8i; t;

ð47Þ

m2M p2P k2K

X XX x00t jkpm m2M p2P k2K

!!

!! ð1  y0j Þ P eð1  y0j Þ;

8j; t;

ð48Þ

119

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

U 00t k 

X XX

x0tikpm þ

m2M p2P i2I

X XX X XX x00t x000t jkpm þ khpm m2M p2P j2J

!! ð1  z0k Þ P eð1  z0k Þ;

8k; t;

ð49Þ

m2M p2P h2H

Constraints (47)–(49) are the non-criticality condition for manufacturing centers, main DCs, and local DCs, respectively. The transformation is achieved using Eqs. (1)–(3). t t 0t 00t 000t t þt Itip ; I0tjp ; I00t kp ; qilp ; xijpm ; xikpm ; xjkpm ; xkhpm ; si ; si ; U P 0; rt t t wi ; yj ; zk ; v tij ; v 0tjk ; v 00t kh ; milp ; eil ; uil 2 f0; 1g;

8i; j; k; h; l; p; m; t;

ð50Þ

8i; j; k; h; l; t:

ð51Þ

Finally, constraints (50) and (51) determine the type of decision variables. 3.3. Calculation of disruption probability (atil ) In order to derive the technology disruption/failure probability atil , we propose a model which not only considers this failure out of natural disasters, but also out of technology malfunction itself (here we call internal). We assume that a technology type can be disrupted by multiple types of causes, whose expected probability for each technology type in each period and manufacturing center can be obtained by integrating the impact of each resource. In general, when joint occurrence of different types of resources are considered, they can be divided into natural and internal sources as follows: Let us consider the set of all natural disaster sources by N , and the joint occurrence of a subset of F denotes a scenario. We define Cs ¼ fF 0 jF 0 # F ; jF 0 j ¼ sg; 8s ¼ 0; 1; . . . ; jF j as the collection of all subsets of natural disaster scenarios of size s, S j where jF s¼0 Cs accounts for all the natural disaster scenarios. The expected disruption probability of technology l out of all possible natural disaster scenarios (atN il ) can be derived through enumerating all the scenarios as follows:

atN il ¼

jF j X X s¼0 F 0 2C

ptil ðF 0 Þvtil ðF 0 Þ;

8i; l; t;

ð52Þ

s

where ptil ðF 0 Þ is the frequency of occurring scenario F 0 as a matter of using technology level l in manufacturing center i in period t, and vtil ðF 0 Þ is the average failure time of technology level l in manufacturing center i in period t when natural disaster scenario F 0 happens. On the other hand, the expected disruption probability of technology l out of all possible internal sources scenarios (atI il ) can be derived as the maximum expected probabilities of all productions with technology l in manufacturing center i in period t;

~ t rt rt atI il ¼ maxful P ilp milp g; ðp;rÞ

8i; l; t;

ð53Þ

Finally, the disruption probability of technology level l in manufacturing center i in period t (atil ) induced by all possible scenarios can be calculated as follows:

T

l∈O

T-1 ξlp(n-1)

r'tlp

T-2 ξlp(n-2)

t=2 t=1 ξlp(1)

l∈N ξ'lp(m-1)

ξlp(2) ξ'lp(2)

ξ'lp(1) t q Fig. 2. Multiple breakpoint unit production emission function.

120

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142 tI atil ¼ maxfatN il ; ail g;

8i; l; t;

ð54Þ

3.4. Unit production gas emission As discussed earlier, reducing emissions from industry requires a sustained and focused effort. One of the major challenges toward this aim is better consideration of technology utilization and production management. In order to significantly reduce the industrial emissions and to maximize the energy efficiency, older and inefficient technologies should be replaced with the current best available technologies and best practice technologies (Brown et al., 2012). In this regard, the unit production emission is considered as a multiple break point function of total production each product type in each period, as opposed to be linear. This allows the ability to capture the increasing behavior of emission with machines’ partial degradation over time and the volume of batch in the production system. In other words, we have r0tlp ðqtilp Þ, which means r0tlp is a function of qtilp with an increasing slope of nlp , where nlp > nlp , r 0tlpð1Þ > r 0tlpð1Þ and m < n (see Fig. 2). l2O

l2N

l2O

l2N

According to above-mentioned descriptions, the unit production emission can be stated as:

r 0tlp ðqtilp Þ ¼

8" !# X X X r;t1 X r;t1 > > 0t1 rt 0t rt > > r lpð1Þ ð1 þ dÞ milp milp þ r lpð1Þ 1  þ nlpð2Þ ðqtilp  Q 1 Þ milp milp > > > r2R r2R r2R r2R > > " !# > > > X X X r;t1 X r;t1 > > 0t1 rt 0t rt < r lp ð1 þ dÞ milp milp þ r lpð2Þ 1  milp milp þ nlpð2Þ ðqtilp  Q 2 Þ ð2Þ r2R

r2R

r2R

r2R

r2R

r2R

r2R

r2R

if Q 1 6 qtilp 6 Q 2 if Q 2 6 qtilp 6 Q 3

> > > .. > > > . > " !# > > > X X X r;t1 X r;t1 > > 0t1 rt 0t rt > þ nlpð2Þ ðqtilp  Q n1 Þ if Q n1 6 qtilp 6 Q n milp milp > : r lpð2Þ ð1 þ dÞ milp milp þ r lpð2Þ 1  ð55Þ

By introducing a positive variable qtilp ¼

r0tlp ðqtilp Þ ¼

Pn1

t n0 ¼1 qilp ðn0 Þ ,

formulation (55) can be transformed as follows:

8 ! ! X r;t1 X r;t1 X X > > 0t1 rt 0t rt > > rlpð1Þ ð1 þ dÞ milp milp þ rlpð1Þ 1  milp milp  nlpð1Þ Q 1 þ nlpð1Þ qtilp ð1Þ > > > r2R r2R r2R r2R > > ! ! > > > X r;t1 X r;t1 X X > > 0t1 rt 0t rt < rlp ð1 þ dÞ milp milp þ rlpð2Þ 1  milp milp  nlpð2Þ Q 2 þ nlpð2Þ qtilp ð2Þ ð2Þ r2R

r2R

r2R

r2R

r2R

r2R

r2R

r2R

if Q 1 6 qtilp ð1Þ 6 Q 2 if Q 2 6 qtilp ð2Þ 6 Q 3

> > > .. > > > . > ! ! > > > X r;t1 X r;t1 X X > > 0t1 rt 0t rt > r ð1 þ dÞ m m þ r 1  m m Q  n þ nlpðn1Þ qtilpðn1Þ if Q n1 6 qtilpðn1Þ 6 Q n > lpðn1Þ n1 ilp lpð2Þ ilp ilp ilp : lpð2Þ ð56Þ

then, we would have:

r0tlp ðqtilp Þqtilp ¼

8 ! ! X X X X r;t1 > 2 > > > r0t1 ð1 þ dÞ mrtilp mr;t1 qtilp ð1Þ þ r0tlpð1Þ qtilp ð1Þ  mrtilp milp qtilp ð1Þ  nlpð1Þ Q 1 qtilp ð1Þ þ nlpð1Þ qtilp ð1Þ > lp ilp ð1Þ > > r2R r2R r2R r2R > > ! ! > > > X X X X r;t1 > 2 > t 0t t < r0t1 mrtilp mr;t1 mrtilp milp qtilp ð2Þ  nlpð2Þ Q 2 qtilp ð2Þ þ nlpð2Þ qtilp ð2Þ lpð2Þ ð1 þ dÞ ilp qilp ð2Þ þ r lpð2Þ qilp ð2Þ  r2R

r2R

r2R

if Q 1 6 qtilp ð1Þ 6 Q 2 if Q 2 6 qtilp ð2Þ 6 Q 3

r2R

> > > . > > > .. > ! ! > > > X X X r;t1 X r;t1 > 2 > 0t1 rt t 0t t rt t t > ð1 þ dÞ m m q þ r q  m m q Q q r  n þ nlpðn1Þ qtilpðn1Þ if Q n1 6 qtilpðn1Þ 6 Q n > lpðn1Þ n1 ilpðn1Þ ilp ilpðn1Þ lpð2Þ ilpðn1Þ ilp ilpðn1Þ ilp ilp : lpðn1Þ r2R

r2R

r2R

r2R

ð57Þ

The mathematical model presented so far shows some non-linear expressions, which significantly increases the complexity of the model and does not guarantee that a global solution will be obtained. To convert the model into its linear counterpart, a linearization section is provided in Appendix A.

4. Solution approach In this section, we are going to address how the considered uncertain parameters should be treated, how to deal with the multi-objectivity aspect of the model, and how to cope with the complexity of the problem.

121

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

4.1. Proposed mean-absolute deviation possibilistic-stochastic programming (MADPSP) model The proposed model is a possibilistic-stochastic model with mixed uncertainty in costs, environmental and social impacts, incident probability and demand, some of which like demand and incident probability are tainted with deeper uncertainty. In such environments, a synergy of fuzzy and stochastic tools can effectively resolve the existing issue (Mohammadi et al., 2014; Liu et al., 2003). The classical possibilistic-stochastic programming method, in a nutshell, is an extension of the fuzzy-robust optimization and chance constraint programming tools. Despite its superiority in being capable of dealing with different kinds of uncertainties (e.g., fuzzy, stochastic, interval numbers, etc.) in terms of random-fuzzy variables, it cannot adequately address the risk factor of the objective function, as in many papers, only the average state of uncertain parameters have been addressed. However, in real world applications, the decision makers (DMs) mainly want to impose the risk factor in the decision making process. In addition, considering the risk factor in resilient logistics is of a concern in many industries (see Waters, 2011). Towards this, we develop the previous classical method to a mean-absolute deviation fuzzy possibilistic-stochastic programming (MADFPSP) model. In doing so, the model should be converted into its crisp counterpart. To start with, consider the compact model as follows:

~ min Z ¼ CX s.t.

~ 6 B; ~ AX xj P 0;

xj 2 X;

j2J

ð58Þ

~ is a fuzzy vector with elements ~cj , and B ~¼ where C

ðpnn Þ ðp1 Þ ðp2 Þ ðb1 ; b2 ; . . . ; bn1 ; bn1þ1 ; bn1þ2 ; . . . ; bn 1 Þ

is an uncertain vector consisting ðp Þ

ðp Þ

ðpnn Þ

1 2 ; bn1þ2 ; . . . ; bn 1 ). both fuzzy (i.e., b1 ; b2 ; . . . ; bn1 ) and random fuzzy elements with deeper uncertainty (i.e., bn1þ1 According to Carlsson and Fullér (2001), the lower and upper possibilistic mean values of fuzzy number ~cj with n-level set

½cj n ¼ ½f cj ðnÞ; g cj ðnÞ; ðn 2 ½0; 1Þ can be defined as Eqs. (59) and (60), respectively:

R1 M ð~cj Þ ¼

0

Z 1 f cj ðnÞPosðcj 6 f cj ðnÞÞdn ¼2 nf cj ðnÞdn R1 0 Posðcj 6 f cj ðnÞÞdn 0

ð59Þ

Z 1 g cj ðnÞPosðcj P g cj ðnÞÞdn ng cj dn ¼2 R1 0 Posðcj P g cj ðnÞÞdn 0

ð60Þ

and

R1 

M ð~cj Þ ¼

0

where Pos denotes possibility and can be defined as:

Posðcj 6 f cj ðnÞÞ ¼ Pð1; f cj ðnÞÞ ¼ sup cj ðuÞ ¼ n;

ð61Þ

u6f c ðnÞ j

and

Posðcj P g cj ðnÞÞ ¼ Pðg cj ðnÞ; þ1Þ ¼ sup cj ðuÞ ¼ n: uPg c ðnÞ

ð62Þ

j

Definition 1 Carlsson and Fullér (2001). For a fuzzy number ~cj with ½cj n ¼ ½f cj ðnÞ; g cj ðnÞ; ðn 2 ½0; 1Þ, the crisp possibilistic mean can be stated as:

Mðcj Þ ¼ ðM  ðcj Þ þ M ðcj ÞÞ=2 ¼

Z

1 0

nðf cj ðnÞ þ g cj ðnÞÞdn

ð63Þ

~ with Definition 2 Zhang and Zhang (2014). Given two fuzzy numbers ~cj with ½cj n ¼ ½f cj ðnÞ; g cj ðnÞ; ðn 2 ½0; 1Þ and d j n

½dj  ¼ ½hcj ðnÞ; kcj ðnÞ; ðn 2 ½0; 1Þ, the absolute deviation between cj and dj can be defined as:

 j þ dj  Mðc  j Þ  Mðd  j ÞÞ=2 rðcj ; dj Þ ¼ ðMjc Considering a trapezoidal membership function for the fuzzy number cj , we would have:

ð64Þ

122

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

lc~j ðxÞ ¼

8 xc1 j > > > c2j c1j > > > > <1

if c1j 6 x 6 c2j

> > > > > > > :

c4j x c4j c3j

if c3j 6 x 6 c4j

0

if x 6 c1j or x P c4j

if c2j 6 x 6 c3j

ð65Þ

Obviously, in case where c2j ¼ c3j , the trapezoidal membership function will be converted to a triangular one. Since the membership function is continuous, all n-cuts are bounded and can be expressed as: cjn ¼ ½c1j þ nðc2j  c1j Þ; c4j  nðc4j  c3j Þ. Accordingly, the lower and upper possibilistic means and the possibilistic mean value can be respectively stated as:

M ð~cj Þ ¼ ð2c2j þ c1j Þ=3;

ð66Þ

M ð~cj Þ ¼ ð2c3j þ c4j Þ=3;

ð67Þ

Mð~cj Þ ¼ ðc1j þ 2c2j þ 2c3j þ c4j Þ=4:

ð68Þ

Theorem 1 Zhang and Zhang (2014). For a trapezoidal fuzzy number ~cj ¼ ðc1j ; c2j ; c3j ; c4j Þ, the possibilistic absolute deviation can be expressed as:

rð~cj Þ ¼ c3j  c2j þ ððc4j  c3j Þ þ ðc2j  c1j ÞÞ=3

ð69Þ

Now, since the right hand side of constraint of model (58) contains an uncertain parameter, including elements of random ðp Þ

ðp Þ

ðpnn Þ

1 2 fuzzy ones (i.e., bn1þ1 ; bn1þ2 ; . . . ; bn 1 ), if condition (70) holds, the possibilistic-stochastic constraint of model (58) can be replaced by the following 2k precise inequalities in which k indicates k levels of n–cut Liu et al., 2003.

flAij ðnij Þjnij 2 ½0; 1g ¼ fni1 ; ni2 ; . . . ; nik g 0 6 ni1 6 ni2 6    6 nik 6 1;

i2I

ð70Þ

We can thus have:

lX 6 B l ; A

8l ¼ 1; 2; . . . ; k

ð71Þ

Al X P Bl ;

8l ¼ 1; 2; . . . ; k

ð72Þ

 l ¼ supðAl Þ, B  l ¼ supðBl Þ, Al ¼ infðAl Þ and Bl ¼ infðBl Þ. where A With respect to the discussed formulations, the MADFPSP counterpart of model (58) can be expressed as:

~ þ ð1  cÞrðCÞX ~ min f ¼ cMðCÞX s.t. J X s sij XÞ 6 B ða i j¼1

( s ¼ with B i

s b i sðp Þ  b i i

J X

when i ¼ 1; 2; . . . ; m1 ; s ¼ 1; 2; . . . ; k1 when i ¼ m1 þ 1; m1 þ 2; . . . ; m; s ¼ k1 þ 1; k1 þ 2; . . . ; k

ðasij XÞ P Bsi

j¼1

( with

Bsi

xj P 0;

¼

bsi sðpi Þ

bi

when i ¼ 1; 2; . . . ; m1 ; s ¼ 1; 2; . . . ; k1 when i ¼ m1 þ 1; m1 þ 2; . . . ; m; s ¼ k1 þ 1; k1 þ 2; . . . ; k

ð73Þ

j2J

where c is a weight factor (0 6 c 6 1) imposing a tradeoff between mean and standard deviation of the uncertain parameters. Noteworthy, the value of this factor is purely subjective and should be determined based on the DM’s preference.

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

123

~ under any n–cut level have random features Additionally, since boundaries intervals of some fuzzy elements of vector B (see Mohammadi et al., 2014), hence, they can be presented as normal distributions as following Eqs.:

( ) 2 ½b2 ðsÞ  l 1 ; p½b2 ðsÞ ¼ pffiffiffiffiffiffiffi exp  2r2 2pr

ð74Þ

( 2)  1 ½b2 ðsÞ  l p½b2 ðsÞ ¼ pffiffiffiffiffiffiffi exp  ; 2 2r  2pr

ð75Þ

and

where

l and r2 are the expected value and variance of b2 ðsÞ, respectively. The same holds for l and r 2 for b2 ðsÞ.

4.2. Approximated lower bound The primary model (‘‘PM” in short) in Section 3.2, even though linearized, is NP-hard in nature and difficult to solve in large scale instances. In order to face that, a two-step multi-objective lower bound procedure based on the relaxed version of the model (‘‘RM” in short) and augmented e-constraint method is presented in this section, followed by a developed evolutionary algorithm in Section 4.3. In the RM, substituting Prtilp with a given value provides a lower bound of objective function value for the PM (Aboolian et al., 2012). Let ati½1 6 ati½2 6    6 ati½L1 6 ati½L be an ordering of failure probabilities in L. We have:

Hrtil ¼

8 ! r Y > > t > a ð1  atil Þ; 8i; l 2 L; r; t; > > < 1¼1 i½1 r1 > Y > > t > > : ai½1 ;

ð76Þ

8i; l ¼ L þ 1; r; t:

1¼1

where Hrtil is an optimistic version of P rtilp 8p 2 P, and provides a lower bound for that (i.e., Hrtil 6 P rtilp - see Aboolian et al. (2012) for proof). The RM can then be stated as:

min ½Z1 ; Z2 ; Z3 ; Z4  s.t.

X

Hrtil qtilp 6 gilp ;

ð14-17Þ

8i; l 2 fL þ 1g; p; t;

ð77Þ

r2R

Constraints (12), (18)–(27), (31)–(45), (50), (51), (A.1)–(A.10), (A.17)–(A.28), (A.34)–(A.36). Theorem 2. Consider xPM and ZPM ðxPM Þ to be the optimal solution and objective function of the PM, respectively, and xRM and ZRM ðxRM Þ to be the optimal solution and objective function of the RM, respectively. ZRM ðxRM Þ is a lower bound for ZPM ðxPM Þ in such a way that:

ZRM ðxRM Þ 6 ZPM ðxPM Þ:

ð78Þ

Proof. Clearly, the optimal solution of PM (i.e., xPM ) is a feasible solution for RM, and since Hrtil 6 P rtilp , we have ZRM ðxPM Þ 6 ZPM ðxPM Þ. Additionally, ZRM ðxRM Þ 6 ZRM ðxPM Þ, since xRM is the optimal solution to the RM. Consequently, we can derive that ZRM ðxRM Þ 6 ZPM ðxPM Þ. h

Now, to provide a Pareto lower bound frontier, first, consider the following multi-objective problem (MOP) with all minimized objectives:

min Z j ðxÞ s.t.

Z 1 ðxÞ 6 e1 ; Z 2 ðxÞ 6 e2 ; . . . ; Z j1 ðxÞ 6 ej1 ; x 2 X:

ð79Þ

where e1 ; e2 ; . . . ; ej1 are the satisfaction levels for the objective functions, and impose a maximum perturbation level for the constrained objectives. By parametrical variations of the e values step-by-step, Pareto optimal solutions can be achieved (see Mousazadeh et al. (2015) for further information). However, the derived solutions from these series of models are usually not efficient (Zhang and Reimann, 2014). In order to face such inefficiency, a couple of modified versions, namely AUGMECON

124

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

proposed by Mavrotas (2009), AUGMECON2 presented by Mavrotas and Florios (2013), and SAUGMECON introduced by Zhang and Reimann (2014) have been developed. In step two of our proposed lower bound approximation method, we utilize the SAUGMECON model. The superiority of this approach over the two other developed ones is aimed at its computational efficiency. Unlike AUGMECON and AUGMECON2 models, SAUGMECON tries to marginally optimize of all the objective functions beside the main one. Accordingly, model (79) can be transformed to:

min Z j ðxÞ þ gðZ 1 ðxÞ=.1 þ Z 2 ðxÞ=.2 þ    þ Z j1 ðxÞ=.j1 Þ s.t.

Z 1 ðxÞ 6 e1 ; Z 2 ðxÞ 6 e2 ; . . . ; Z j1 ðxÞ 6 ej1 ; x 2 X:

ð80Þ 3

where g is an adequately small number usually between 10

6

and 10 , and

.j is the range of the jth objective function.

4.3. Developed evolutionary algorithm To solve NP-hard SCND problems, different techniques in the literature have been proposed including Lagrangian relaxation, Benders decomposition algorithm, meta-heuristics and evolutionary algorithms (see Zhalechian et al., 2016; Zahiri et al., 2014; Mohammadi et al., 2014; Niakan et al., 2015; Shaw et al., 2016; Diabat et al., 2015; Tong et al., 2015; Jeihoonian et al., 2016; Keyvanshokooh et al., 2016). In this paper, we propose a new multi-objective evolutionary algorithm called DVG, based on Differential Evolution (DE) algorithm, Variable Neighborhood Search (VNS) algorithm and game theory. For more information regarding the application of game theory in MOPs, see Mohammadi et al. (2016). 4.3.1. Game-based EA and rewarding procedure In this algorithm, in order to optimize the objective functions of MOPs, an evolutionary game is utilized. In order to reward each solution in the population, a procedure based on the percentage of victory during the game is utilized. This is further determined due to the number of solutions dominated by each solution. The first step in designing a game-based EA is generating a population corresponding to each objective function. Next, each solution (i.e. player) of each population plays with the players of the other populations and vice versa are gathered. The reward value ðRÞ of each player in competition with other populations can be derived from Eqs. (81)–(83). Finally, based on the obtained reward values, each player gets a ‘‘fitness”, which is calculated as Eq. (84). Fig. 3 shows the population architecture for a game-based evolutionary algorithm.



w1

OFV 2 ðxj ÞOFV 2 ðxi Þ

max OFV 1

max OFV 2





OFV 3 ðxj ÞOFV 3 ðxi Þ

OFV 4 ðxj ÞOFV 4 ðxi Þ

max OFV 3

max OFV 4

þw4

ð81Þ ð82Þ



OFV 1 ðxl ÞOFV 1 ðxi Þ OFV 2 ðxl ÞOFV 2 ðxi Þ OFV 3 ðxl ÞOFV 3 ðxi Þ OFV 4 ðxl ÞOFV 4 ðxi Þ þw2 þw3 þw4 max OFV 1 max OFV 2 max OFV 3 max OFV 4

...

... ...

ð83Þ

...

...

...

w1

OFV 1 ðxk ÞOFV 1 ðxi Þ OFV 2 ðxk ÞOFV 2 ðxi Þ OFV 3 ðxk ÞOFV 3 ðxi Þ OFV 4 ðxk ÞOFV 4 ðxi Þ þw2 þw3 þw4 max OFV 1 max OFV 2 max OFV 3 max OFV 4

...

 Rði; lÞ ¼ e

þw3



w1

Rði; kÞ ¼ e

þw2

...

Rði; jÞ ¼ e

OFV 1 ðxj ÞOFV 1 ðxi Þ

Fig. 3. Population architecture for a game-based evolutionary algorithm.

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

PN

j¼1 Rði; jÞ

FðiÞ ¼

N

PN þ

k¼1 Rði; kÞ

N

125

PN þ

l¼1 Rði; lÞ

N

ð84Þ

where xi , xj , xk and xl are the solutions in the first to forth populations, respectively, w1 to w4 are the importance coefficients of the objective functions (i.e. w1 þ w2 þ w3 þ w4 ¼ 1), and N is the population size. 4.3.2. Differential evolutionary algorithm Differential evolution (DE), introduced by Storn and Price (1997), is an efficient evolutionary heuristic, which belongs to the category of evolutionary algorithms ðEAsÞ and is used to obtain near-optimal solutions using operators (i.e., crossover and mutation). DE, similar to other EAs, starts with a randomly chosen initial population over the variable domain. In each generation, mutation operator is deployed to create a new solution vector as Eq. (85).

MV i ¼ xv 1 þ Fðxv 2  xv 3 Þ

ð85Þ

where MV i is a mutant vector, xv 1 , xv 2 and xv 3 are different random vectors selected from ði ¼ 1; . . . ; N pop Þ, which are not equal to i, and F is the amplification coefficient (i.e. F 2 ½0; 2). The diversity of the algorithm is increased using a binomial crossover operator (Storn and Price, 1997). In doing so, each mutated vector shares its information with another predetermined (target) vector to generate trial vector li ¼ fli1 ; . . . ; lik ; . . . ; lin g as follows:



lik ¼

MV ik ; if randðkÞ 6 CR and k ¼ rnbrðiÞ

xik ;

if randðkÞ > CR and k–rnbrðiÞ

ð86Þ

where randðkÞ denotes the k-th element of an N-dimensional uniform random number 2 ½0; 1, CR shows the crossover rate and rnbrðiÞ is a randomly chosen index 2 f1; . . . ; Ng, which ensures that the new solution vector gets at least one of the dimensional values of mutated vector. Finally, if the value of the fitness function for the solution vector is lower, the target vector will be replaced by the new solution in the following generation. 4.3.3. Proposed DVG algorithm As previously mentioned, the proposed DVG algorithm is based on Differential Evolution (DE) algorithm, Variable Neighborhood Search (VNS) algorithm and game theory. In this algorithm, the optimization process starts with generating four random populations. Afterwards, the solutions are ranked based on their fitness values obtained from the rewarding procedure (see Section 4.3.1). In the next step, the ranked solutions are evolved using DE algorithm explained in Section 4.3.2. Next, in order to have local search, a neighborhood search structure (NSS) procedure is applied. In the proposed algorithm, NSS is applied on the solutions by utilizing VNS to achieve better solutions. With proven merits and ever-increasing successful application of VNS, the importance of such algorithm has been highlighted in recent years (Hansen and Mladenovic´, 2001). The superiority of VNS over other local search heuristics is to use two or more neighborhoods, instead of one in its structure. More precisely, during the search process, the neighborhood undergo a systematic change (Behnamian et al., 2009). Additionally, to avoid spending overwhelming computational time, the best number of neighborhoods is often three (Mohammadi et al., 2013), as we followed in our algorithm. The three neighborhoods employed in our algorithm are as follows: I. Swap: In the swap operator, places of two random selected parts are exchanged. II. Reversion: In the reversion operator, a random part of a solution is selected and its permutation is reversed. III. Inversion: In the inversion operator, one bit is chosen randomly and its value is replaced by a new random value. Furthermore, the number of outer and inner iterations of the proposed VNS are equal to ItVNSout and ItVNSin , respectively. The termination criteria is the number of function calls (NFC). Finally, the Pseudo code of the proposed DVG algorithm is provided as Fig. 4. After applying VNS on the Pareto solutions, the new solution is again rewarded using Eqs. (81)–(83) to be compared with the current solution. 5. Computational experiments In this section, first, in order to show the efficiency of the proposed DVG algorithm, several small, medium and large scale instances are presented in Benchmarking subsection. The considered case study is presented in Case study section followed by sensitivity analyses and managerial insights. 5.1. Benchmarking In this section, the performance of the propose DVG algorithm in terms of solution quality and computational time is compared with two other widely used algorithms, NSGA-II and MOICA, under different randomly generated data sets (for

126

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Fig. 4. The pseudo code of the proposed DVG.

further information about these two algorithms, see Mohammadi et al., 2013). These data sets contain 60 different problems ranging from small, medium and large scale, 20 of each. In order to tune the parameters of the algorithms, the response surface methodology (RSM) has been utilized (see Calégari et al., 1999) and the values are provided in Table A1 in Appendix A. In addition, parameter ranges are provided in Table A2 in the same file. It should be noted that these parameter settings, before being used in the comparative study, were pre-tested over randomly selected problem instances. In order to evaluate the solutions out of the DVG, NSGA-II and MOICA, five commonly used metrics in the literature, namely (i) quality metric (QM), (ii) convergence metric (CM), (iii) divergence metric (DM), (iv) spacing metric (SM), and (v) mean ideal distance (MID) metric have been used. With aid of these metrics, both qualitative and quantitative comparisons between multi-objective EAs can be addressed. The presented model and its relaxed version were coded in GAMS software. We compared the results of these two models in small scale and some medium size instances, since with increasing the size of the problem, due to its NP-hard nature, obtaining the optimal solution from the PM was impossible. To provide a notation for each, the results out of the primary model without any relaxation and lower bound approximation model are referred to PM and RM, respectively. The results are provided in Tables 2 (below) and A3 and A4 in Appendix A. As can be seen in Tables 2 and A3, the CPU time for the PM significantly rises from 118 to 19,982 as the size of problem instances increases. This range for the RM is between 31 and 2330 for the corresponding problem instances. Due to low memory of GAMS software, the larger instances could not be solved under the PM. The mean tightness of the derived solutions (see CM column), which is the gap between the results out of the PM and RM, is equal to 1.068%. The closeness to the optimal solution, in terms of low gap, and low CPU time guarantee both the quality of the solutions as well as their efficiency. Finally, for the large instances, which are provided in Table A4, the maximum CPU time for the RM is equal to 2.77 h (9972 s).

127

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142 Table 2 The metric results out of small scale instances. Problem no.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

PM

RM

DM

SM

MID

CPU (s)

CM

DM

SM

MID

CPU (s)

0.636 0.42 0.895 0.319 0.805 0.624 0.688 0.779 0.656 0.406 0.582 0.367 0.667 0.733 0.78 0.606 0.587 0.724 0.86 0.564

0.535 0.442 0.767 0.847 0.761 0.477 0.897 0.697 0.902 0.419 0.658 0.462 0.814 0.739 0.667 0.872 0.732 0.834 0.31 0.65

0.605 0.852 0.934 0.846 0.936 0.903 0.894 0.546 0.903 0.77 0.712 0.874 0.583 0.928 0.616 0.608 0.936 0.839 0.768 0.591

118 148 226 333 391 559 603 611 830 1391 1662 1833 2212 4492 5511 6329 7836 8238 8882 9103

0.553 0.641 0.783 0.741 0.777 0.999 1.149 1.086 1.092 1.144 1.129 1.147 0.666 1.369 1.097 1.09 1.43 1.138 1.119 1.426

0.929 1.464 1.418 1.411 1.053 1.193 1.273 0.911 0.984 1.243 1.352 1.336 0.812 1.181 1.326 0.985 1.153 1.308 0.922 0.959

0.474 0.275 0.716 0.79 0.565 0.253 0.824 0.655 0.89 0.277 0.55 0.36 0.767 0.535 0.441 0.763 0.668 0.757 0.261 0.483

0.449 0.703 0.542 0.771 0.813 0.702 0.877 0.418 0.816 0.573 0.634 0.518 0.542 0.645 0.605 0.436 0.781 0.737 0.707 0.452

31 63 71 83 99 121 196 277 301 333 386 411 557 599 654 724 803 872 994 1048

Now, in order to validate the performance and efficiency of the proposed DVG algorithm compared to NSGA-II, MOICA and RM in terms of CM, DM, SM and MID metrics, the experiments were re-conducted. The solutions out of these algorithms for the small, medium and large test problems are provided in Tables 3–5. Noteworthy, the values of metrics are calculated compared to the RM. The results in Tables 3–5 show the higher performance of DVG algorithm compared to NSGA-II and MOICA considering lower values for CM, SM and MID and higher values for DM. Additionally, lower CPU time for DVG compared to the two other metaheuristics guarantees the efficiency of the algorithm. Considering CM metric, the mean gap between the proposed DVG and RM is 2.07%, while this value for NSGA-II and MOICA is 4.03% and 4.022, respectively. With increasing the size of the problem, the CPU time for NSGA-II starts from 15 s and for the largest problem, equaling to 3113. The range for MOICA is between 17 to 2881 s, while for DVG, the CPU time varies from 11 to 1249 for all tested problems. The efficiency of the algorithms under the tested problem instances are illustrated in Fig. 5. As can be seen, comparing NSGA-II and MOICA, for small and medium scale instances, NSGA-II shows higher efficiency, while for large scale problems, MOICA outperform NSGA-II in terms of CPU time. Now, in order to schematically evaluate the quality performance of the solutions out of the proposed DVG algorithm, Pareto fronts (PF) of all three algorithms along with the RM model for objectives Z 1 and Z 2 (a), Z 1 and Z 4 (b), Z 2 and Z 3 (c), Z 2 and Z 4 (d), and Z 3 and Z 4 (e) for problem #35 are provided in Fig. 6. As the results show, the proposed DVG algorithm provides high quality non-dominated solutions, which proves its superiority over the two other meta-heuristic algorithms.

5.2. Case study Now, with the proposed DVG algorithm and its proved efficiency, the model can be implemented on a real size problem. Ranked sixth in the world, in terms of medical concentration of healthcare providers, France is among the biggest consumers of pharmaceutical products in the world and one of the most complex markets to distribute pharmaceutical goods. France has the world’s highest consumption of medicines per capita at 300 USD. Additionally, with a population of 65 million, France is the second largest European market for medical care (http:finesco-medica.com). According to National Institute of Statistics and Economic Studies (INSEE), the French Ministry of Industry (FMI) and the French Statistical Office (FSO), the annual turnover of French pharmaceutical sector is equal to €36 billion. In terms of expenditure, pharmaceutical share has risen rapidly in France, as in other European countries, more than doubling from 1990s to 2000s and 2010s (Sermet et al., 2010). This further signifies the inefficient pharmaceutical network in France. HIV exposure continues to be a major global public health issue, having claimed >35 million infections so far. According to World Health Organization (WHO), in 2015, 940,000–1.3 million people died from HIV-related causes globally, 34.0– 39.8 million people living with HIV and 1.8–2.4 million people becoming newly infected (http:who.int). This is while Western Europe account for a further quarter of this number - France (8%), Spain (6%), United Kingdom (5%) and Italy (5%) http: avert.org. Newly reported statistics by WHO reveals a 35% decrease in new HIV infections, and HIV-related deaths by 28% between 2000 and 2015, as a result of international efforts that led the global achievement of the HIV targets of the Millennium Development Goals. However, the number of new infections in France has remained stable despite a decline in other countries (http:who.int; http:reuters.com). Among different categories, rates of HIV infection among men who have sex with

128

Table 3 Comparison of algorithms under small scale instances.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

DVG

NSGA-II

MOICA

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

0.8 1 1 0.9 1 1 0.8 1 1 1 0.7 1 1 0.95 1 1 1 1 1 1

0.662 0.751 0.762 0.817 0.838 0.885 0.967 0.982 1.015 1.047 1.1 1.058 1.29 1.29 1.314 1.225 1.545 1.613 1.621 1.699

1.504 1.33 1.637 1.268 1.171 1.194 0.993 1.204 1.649 1.382 1.37 1.473 0.857 1.317 1.42 1.412 1.721 0.975 1.363 1.301

0.417 0.248 0.378 0.347 0.24 0.288 0.251 0.312 0.314 0.241 0.234 0.327 0.118 0.201 0.31 0.442 0.266 0.392 0.21 0.246

0.655 0.548 0.407 0.428 0.497 0.549 0.515 0.611 0.32 0.426 0.233 0.452 0.433 0.442 0.222 0.286 0.518 0.433 0.493 0.566

11 23 83 90 92 94 117 124 126 135 136 261 278 302 329 339 349 356 369 393

0.2 0 0.2 0 0 0.2 0 0 0.1 0 0 0.3 0 0 0.05 0 0 0 0 0

1.123 1.149 1.24 1.243 1.398 1.536 1.823 1.657 1.869 1.937 2.077 2.12 2.125 2.27 2.091 2.392 2.432 2.154 2.602 2.632

0.958 0.76 1.275 0.796 0.543 1.127 0.888 1.114 1.297 1.159 0.995 0.946 0.705 0.768 1.226 1.066 1.076 0.659 0.921 0.948

0.439 0.409 0.952 0.717 0.869 0.563 0.513 0.532 0.451 0.804 0.604 0.712 0.61 0.856 0.697 0.967 0.405 0.623 0.855 0.328

0.787 0.901 0.662 0.666 0.657 0.778 0.93 0.746 0.508 0.893 0.585 0.867 0.83 0.803 0.529 0.502 1.099 0.703 0.82 0.863

15 48 111 304 341 352 394 531 543 560 618 637 656 663 703 746 756 762 765 810

0.2 0 0.05 0 0 0 0 0 0.1 0 0 0 0.1 0 0 0 0 0 0 0

1.078 1.333 1.464 1.467 1.517 1.543 1.334 1.723 1.725 1.837 2.042 2.045 2.063 2.302 2.41 2.424 2.62 2.621 2.445 2.647

0.501 1.177 0.732 1.265 0.569 1.063 0.57 0.533 0.871 0.61 0.548 0.782 0.849 0.609 0.7 0.632 0.534 0.784 1.116 1.285

0.635 0.996 0.612 0.736 0.596 0.596 0.252 0.56 0.405 0.807 0.359 0.469 0.218 0.213 0.517 0.608 0.516 0.619 0.281 0.881

0.855 0.868 0.525 0.703 0.572 0.732 0.881 0.728 0.542 0.877 0.692 0.947 0.94 0.738 0.555 0.565 0.878 0.609 0.81 0.908

17 56 165 449 453 565 621 720 776 790 840 874 888 899 951 1002 1057 1118 1152 1158

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Problem no.

Table 4 Comparison of algorithms under medium scale instances.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

DVG

NSGA-II

MOICA

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

1.717 1.848 1.88 1.809 1.94 1.952 1.97 1.971 1.987 1.792 2.02 2.033 2.093 2.11 2.116 2.169 2.137 2.303 2.386 2.39

0.758 1.699 1.426 1.123 1.258 1.647 1.402 1.379 1.123 1.511 1.305 1.37 1.578 1.321 1.129 1.595 1.37 1.305 1.443 1.324

0.24 0.219 0.201 0.525 0.365 0.228 0.414 0.535 0.331 0.236 0.42 0.236 0.51 0.417 0.222 0.552 0.347 0.354 0.326 0.414

0.545 0.441 0.598 0.645 0.621 0.565 0.468 0.404 0.445 0.354 0.404 0.365 0.382 0.527 0.663 0.429 0.438 0.412 0.337 0.542

399 410 416 428 430 459 494 499 508 529 548 577 587 604 676 728 728 742 792 798

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

2.647 2.74 2.89 3.182 3.419 3.501 3.47 3.602 3.604 3.617 3.628 4.001 4.051 4.295 4.155 4.631 4.658 4.879 4.931 5.01

0.682 1.179 1.033 0.505 1.244 1.007 0.837 1.085 0.763 0.611 0.659 0.916 1.053 0.966 1.127 0.661 0.969 1.032 0.906 0.698

0.592 0.511 0.326 0.774 0.714 0.395 0.443 0.901 0.414 0.788 0.99 0.306 0.61 0.72 0.777 0.595 0.782 0.94 0.452 0.434

0.772 0.775 0.706 1.048 0.917 0.819 1.057 0.604 0.748 0.523 0.94 0.86 0.694 0.937 1.013 0.911 0.604 0.848 0.539 0.778

836 862 962 1022 1083 1096 1113 1128 1143 1175 1180 1262 1275 1455 1478 1702 1949 1999 2128 2160

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

2.74 2.752 2.84 2.847 2.868 2.908 3.004 3.306 3.121 3.332 3.59 3.731 4.106 4.448 4.651 4.572 4.873 4.94 5.035 5.177

0.739 1.098 0.688 0.569 0.802 1.041 0.768 0.551 0.664 0.526 0.893 0.951 0.512 1.265 0.788 0.903 0.678 1.287 0.938 0.77

0.786 0.427 0.227 0.641 0.853 0.864 0.802 0.735 0.578 0.358 0.648 0.313 0.963 0.747 0.293 0.673 0.348 0.506 0.463 0.723

0.617 0.781 0.884 1.021 1.011 0.744 0.814 0.813 0.601 0.545 0.803 0.796 0.609 1.045 0.951 0.908 0.685 0.954 0.664 0.771

1255 1259 1311 1366 1377 1424 1450 1472 1570 1573 1590 1598 1615 1658 1713 1756 1834 1844 1856 1866

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Problem no.

129

130

Table 5 Comparison of algorithms under large scale instances.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

DVG

NSGA-II

MOICA

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

2.441 2.467 2.384 2.622 2.661 2.726 2.754 2.798 2.82 2.865 2.98 3.045 3.118 3.191 3.511 3.601 3.438 3.851 3.904 3.966

1.212 1.406 1.625 1.508 1.518 1.372 1.419 1.294 1.48 1.464 1.078 1.335 1.108 1.072 1.523 1.513 1.22 1.461 1.532 1.614

0.402 0.339 0.333 0.278 0.398 0.228 0.234 0.367 0.322 0.395 0.201 0.119 0.222 0.225 0.115 0.378 0.318 0.201 0.311 0.474

0.387 0.423 0.625 0.465 0.323 0.227 0.434 0.542 0.443 0.611 0.42 0.625 0.34 0.332 0.467 0.423 0.504 0.574 0.442 0.423

799 800 809 830 849 850 865 879 887 935 940 980 984 1083 1084 1100 1105 1130 1149 1249

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

5.042 5.56 5.574 5.702 5.72 5.904 6.004 6.079 6.096 6.182 6.035 6.32 6.512 6.537 6.89 7.051 7.333 7.235 7.714 7.928

1.127 0.926 1.151 0.663 1.195 1.209 0.544 0.901 0.707 0.709 0.679 1.052 0.683 0.809 1.202 0.826 0.556 1.212 0.751 0.991

0.623 0.604 0.591 0.313 0.497 0.389 0.729 0.537 0.348 0.74 0.424 0.441 0.27 0.267 0.203 0.476 0.388 0.706 0.425 0.858

0.645 0.719 0.99 0.888 0.644 0.527 0.899 0.936 1.096 0.872 0.827 1.054 0.639 0.591 0.757 0.824 0.715 0.814 0.639 0.831

2182 2190 2288 2316 2344 2377 2507 2568 2571 2575 2739 2747 2812 2843 2916 2947 3014 3020 3027 3113

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

5.185 5.185 5.39 5.493 5.541 5.646 5.702 5.941 6.12 6.291 6.362 6.374 6.747 6.886 7.017 7.326 7.583 7.47 7.702 7.895

0.779 0.802 0.938 0.707 0.591 1.039 1.239 1.232 1.19 0.945 0.772 0.51 0.878 1.047 1.086 0.993 1.154 0.617 0.51 0.765

0.61 0.784 0.403 0.98 0.892 0.899 0.298 0.826 0.552 0.458 0.224 0.219 0.86 0.901 0.202 0.634 0.392 0.294 0.802 0.603

0.642 0.742 0.841 0.897 0.69 0.53 0.963 1.041 0.835 0.82 0.611 0.94 0.509 0.572 0.884 0.719 0.716 0.855 0.691 0.78

2085 2123 2149 2157 2178 2225 2261 2290 2310 2312 2400 2418 2465 2508 2540 2586 2753 2770 2868 2881

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Problem no.

131

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

3500

DVG NSGA-II MOICA

3000

CPU time (s)

2500 2000 1500 1000 500 0 1 5 9 13 17 21 25 29 33 37 41 45 49 53 57

Problem instance

RM DVG NSGA-II MOICA

0.7 0.6

Z4

Z2

0.3

RM DVG NSGA-II MOICA

200 150

30

100

20

0.2

50

10

0.1

0

0

0

Z1 10000

77

RM DVG NSGA-II MOICA

67 57

4000 6000 Z1

60

8000

0 0.00

0.50

Z2

1.00

RM DVG NSGA-II MOICA

50

30

37

20

27

10

17 0.00

2000

40 Z4

47

x 10000

0

Z4

RM DVG NSGA-II MOICA

50 40

0.5 0.4

60

Z3

0.8

x 10000

Fig. 5. Efficiency of the algorithms under different problem instances.

0

0.20 Z2

0.40

0

20

Z3

40

60

Fig. 6. Pareto front solutions for problem #35.

men (‘‘MSM” in short) remain high almost everywhere - account for 63% of new HIV infections annually - and new prevention options are urgently needed (http:who.int; http:hrc.org). According to French National Institute for Public Health Surveillance, nearly half of newly infected people with HIV in France were gay men, and the incidence among homosexual men is 200 times higher than in the heterosexual population (http:reuters.com). Currently, there is not an absolute cure for HIV. Nevertheless, effective antiretroviral (ARV) drugs can control the virus and help prevent transmission (http:who.int). One of the most widely used ARV drugs is TruvadaÓ. The use of TruvadaÓ as a pre-exposure prophylaxis (PrEP) has become the topic of conversation for sexually active people in the LGBTQ community and beyond (http:lgbtweekly.com). PrEP is a strategy for people who do not have HIV, but who are at risk of getting it, to prevent them from getting HIV by taking a single pill (usually a combination of two antiretrovirals – as in TruvadaÓ) daily. Results show that people who are at high risk of infection (most especially in MSM), if are under PrEP, can reduce the risk of HIV infection up to 92% (http:who.int).

132

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

To address these issues, we consider the TruvadaÓ chain as our case study, focusing on the LGBTQ community (most especially gay and bisexual men) demand, as they account for the majority of this population. Additionally, as TruvadaÓ is completely imported to France, having such a chain would significantly help France to be self-dependent on meeting demands. The French health and medical distribution network (FDN) involves (i) manufacturers/laboratories (secondary manufacturers in our model), (ii) wholesalers (main DCs in the model), (iii) distributors (local DCs), and (iv) hospitals/pharmacies (demand points). To as potential locations to establish manufacturing centers, main and local DCs, as well as market points, resulting in a j48j  j48j  j48j  j48j  j5j  j5j  j1j  j3j  j24j problem size (see Fig. 7 for marked locations). It should be noted that the planning horizon is considered to be T ¼ 24, with a period of one month per unit. The value/range of the input parameters for the case study are provided in Table A5 in Appendix A. It should be noted that for uncertain parameter, the value/range shows the most likely points (c2j and c3j ) – see Section 4.1 for more details – and the prominent points are achieved with 5% perturbation over the most likely ones (i.e., ðc1j ; c2j ; c3j ; c4j Þ ¼ ð0:95c2j ; c2j ; c3j ; 1:05c3j Þ). 5.2.1. Sensitivity analysis and final results In this section, some sensitivity analyses are conducted to show the behavior of the objective functions, with respect to a change in storage capacity levels of manufacturing center ðCapÞ, main DC ðCap0 Þ, local DC ðCap00 Þ, production capacity ðgÞ, probability of incident ðuÞ, disruption probability ðaÞ, allowable emission threshold ðGÞ, maximum selling carbon credit ðS Þ, and maximum purchasing carbon credit ðSþ Þ. The results are provided in Tables 6–9, for objective functions Z 1 , Z 2 , Z 3 , and Z 4 , respectively. The illustrative version of the results is given in Fig. 8(a)–(d). It should be noted that the uncertainty

Fig. 7. Cities of France.

Table 6 Changes in Z 1 (%) vs. changes in parameters values. Changes in parameter value (%)

20 15 10 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100

Changes in Z 1 (%) Cap

Cap0

Cap00

g

u

a

G

S

Sþ

9.6 6.3 1.2 1.4 5.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4 7.4

15.7 15.1 9 1.2 4.8 8.5 19.6 25.2 26.1 27.2 27.3 32.5 32.5 32.5 32.5 32.5 32.5 32.5 32.5 32.5 32.5 32.5

16.4 10.2 4.3 5.7 8.4 9.1 10.7 14 21.3 22.4 23.7 25.2 34.3 36.1 36.1 36.1 36.1 36.1 36.1 36.1 36.1 36.1

10 8.4 6.3 0.5 0.9 2.8 4 8.6 12.3 33 33 33 33 33 33 33 33 33 33 33 33 33

16.6 15.3 15.1 12 13.6 15.8 17.1 17.8 23.2 32.8 35.1 45.6 45.6 45.6 45.6 45.6 45.6 45.6 45.6 45.6 45.6 45.6

15.8 15.5 11.4 13.4 14.8 15.7 17.9 18.3 22.1 23.4 25.5 26.2 27.9 31.9 37.7 37.7 37.7 37.7 37.7 37.7 37.7 37.7

6.4 5.9 5.5 0.7 1.8 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3 4.3

8.6 7.5 7.2 3.3 3.9 5.1 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9

8.2 7.1 4.4 3.1 4.3 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

133

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142 Table 7 Changes in Z 2 (%) vs. changes in parameters values. Changes in parameter value (%)

20 15 10 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100

Changes in Z 2 (%) Cap

Cap0

Cap00

g

u

a

G

S

Sþ

10.1 8.4 7.3 11.4 11.7 17.3 18.5 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7 19.7

31.2 26.8 17.3 13.8 16 27.2 28.7 33.1 33.4 33.4 33.4 33.4 33.4 33.4 33.4 33.4 33.4 33.4 33.4 33.4 33.4 33.4

28.3 21.5 18.8 13.3 17 18.6 21 29.5 32.7 32.7 32.7 32.7 32.7 32.7 32.7 32.7 32.7 32.7 32.7 32.7 32.7 32.7

33.2 29.9 27.1 14 20.3 21.7 22.2 22.6 26 27.8 27.8 30.2 34.3 41.7 44.6 61.4 61.4 61.4 61.4 61.4 61.4 61.4

9.2 7.4 3 4.2 5.4 6.5 8.3 9.1 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2

13.1 12.6 7 3.3 5 9.1 14 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3 16.3

3.3 2.5 1.9 1.7 2.5 3.3 4.4 4.9 5.5 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9 5.9

7.3 4.9 2.1 2.3 3.1 3.9 5.1 7.7 7.9 9.3 9.7 10.2 11.1 12.7 12.7 12.7 12.7 12.7 12.7 12.7 12.7 12.7

6.5 4 2.2 5 5.4 6.5 7.4 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5 9.5

Table 8 Changes in Z 3 (%) vs. changes in parameters values. Changes in parameter value (%)

20 15 10 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100

Changes in Z 3 (%) Cap

Cap0

Cap00

g

u

a

G

S

Sþ

23.3 15.4 11.8 13.2 15 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3 17.3

30 26.1 19.9 13 13.7 26.2 27.7 28.3 28.7 32.5 39.1 39.2 41 49.1 49.1 49.1 49.1 49.1 49.1 49.1 49.1 49.1

29.1 22.3 20.2 29.2 34.1 36 36.1 36.4 39 39.5 40.1 44.4 44.8 45.7 50 52 52.8 52.8 52.8 52.8 52.8 52.8

35.8 21.6 17.3 17.8 19 20.3 21.6 23.4 26.7 35.8 39.3 41.1 42.7 44.4 47.7 47.7 47.7 47.7 47.7 47.7 47.7 47.7

34 23.3 11.1 10 12.8 14.2 18.7 20.2 21.7 22.3 38.3 42.1 42.8 46.5 47.1 49 49 49 49 49 49 49

23.1 13.3 9.4 13.7 15.5 23.5 28.5 31 45.2 46.8 51.6 51.6 51.6 51.6 51.6 51.6 51.6 51.6 51.6 51.6 51.6 51.6

22.3 19 11.1 20 22.6 25.5 32.8 34.4 35.6 43.4 51 52 52.9 52.9 52.9 52.9 52.9 52.9 52.9 52.9 52.9 52.9

20.1 17.8 13.3 16.5 19.8 20.3 27.4 27.8 29 35 38.7 39.5 39.5 39.5 39.5 39.5 39.5 39.5 39.5 39.5 39.5 39.5

19.3 15.5 11 10.2 19.2 20.6 22.4 25.6 30.2 37.9 39.9 40.1 42.3 42.3 42.3 42.3 42.3 42.3 42.3 842.3 42.3 42.3

level for all analyses is set to n ¼ 0:9. The overall scheme shows that with increasing the values of u and a, there is an increasing behavior for all objective functions, while with increasing Cap, Cap0 , Cap00 , g, G, S , and Sþ , we see a decrease in the objective values. According to Table 6 and Fig. 8(a), among capacity levels of nodes (i.e., Cap, Cap0 , Cap00 , g), change in capacity level of local DCs, and manufacturing centers has the highest and lowest effect on Z 1 , respectively. The high change in the objective function as a result of change in u and a shows the high effect of these two parameters in the network’s overall cost. For Z 2 , the highest and lowest effect belongs to production capacity (g) and maximum allowable emission (G), respectively. The rationale behind this is that with increasing the production capacity, the total number of established manufacturing centers would decrease and accordingly, the job opportunity declines significantly. That further justifies the high impact of Cap0 and Cap00 on Z 2 . For the third objective function, we see highest impact for G, S , and Sþ compared to the other objectives, as they directly deal with environmental issues. Finally, for Z 4 , we see significant effect of capacity levels

134

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Table 9 Changes in Z 4 (%) vs. changes in parameters values. Changes in Z 4 (%)

20 15 10 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100

60

Changes in Z1 (%)

40

Cap Cap" φ G S+

Cap

Cap0

Cap00

g

u

a

G

S

Sþ

13.7 12.3 9.1 5.5 7.8 9.3 10.1 11.8 12.2 12.3 815.5 17.7 17.9 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2

22.7 21.3 11 10.5 17.6 21.1 22.6 22.6 25.7 27.6 29.2 31.3 33.4 33.6 34.2 35 43.7 43.7 43.7 43.7 43.7 43.7

21.2 18.2 13.1 10.4 11.1 23.5 28 32.8 35.4 42.5 43.2 47.8 50.6 58.3 62.8 63.4 66.8 68.7 68.7 68.7 68.7 68.7

23 11.3 7.7 12.1 22.7 26.1 27.9 28.7 35.3 36.1 39.6 42.1 48.5 49.8 51.2 56.5 56.8 59.6 59.6 59.6 59.6 59.6

21.2 14.3 10.3 12.7 12.9 22.2 23.5 23.6 26.3 28.6 28.8 30.5 35.3 36.8 37.2 39.6 40.8 41.1 41.8 42 47 47

39.9 23.7 11.2 10 21.6 25.2 26.4 27.8 27.9 32 39.5 42.7 47.6 49.3 50.5 50.5 50.5 50.5 50.5 50.5 50.5 50.5

9.9 7.4 5.9 10.7 11.1 11.3 15.1 15.8 15.9 16.3 16.7 18.2 18.4 18.7 20.7 21 21.8 24.3 25.2 25.2 27.5 27.5

5.3 4.9 4.5 6.9 7.3 10.7 11.6 11.8 15 15.8 16.1 18.6 19.5 22.3 22.3 22.3 22.3 22.3 22.3 22.3 22.3 822.3

6.9 6.1 3.7 7.5 10.6 11.4 13.3 13.8 14.7 14.7 15.8 15.9 16.4 18.5 18.5 18.5 18.5 18.5 18.5 18.5 18.5 18.5

40

Cap' η a S-

20

20 0 -20 -10 15 25 35 45 55 65 75 85 95 -20

Changes in Z2 (%)

Changes in parameter value (%)

-40

-80

Changes in parameter value (%)

60

60

40

40

20 0

Cap Cap" φ G S+

Cap' η a S-

-20 -10 15 25 35 45 55 65 75 85 95 -20 -40 -60

Changes in Z4 (%)

Changes in Z3 (%)

-20 -10 15 25 35 45 55 65 75 85 95 -20

-60

-40 -60

0

Cap' η a S-

Changes in parameter value (%)

20 0 -20 -10 15 25 35 45 55 65 75 85 95 -20 -40 -60

Changes in parameter value (%)

Cap Cap" φ G S+

-80

Cap Cap" φ G

Cap' η a S-

Changes in parameter value (%)

Fig. 8. Changes in objective functions (%) vs. changes in parameters values (%).

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

135

Fig. 9. Optimum TruvadaÓ network in France.

(Cap, Cap0 , Cap00 , and g) as well as disruption and incident probabilities (a and u). With increasing the capacity levels, the values for node criticality, node complexity and flow complexity measures reduce significantly, while with increasing the disruption and incident probabilities, the old technology measure increases, and the production system has to be assigned to different technology levels in order to meet the demand. The Optimum TruvadaÓ network in France is provided in Fig. 9. With such network, the total cost is equal to Z 1 ¼ €2228 million, with €1195 million for the fixed cost of facilities establishment, and €1033 million for variable cost (e.g., transportation, inventory, etc.), which further signifies the fitting balance between fixed and variable costs. With establishment of such facilities, up to 1497 job opportunities at a first step will be provided, which is quite significant for production of a single type product. It should be noted that if the manufacturing center/DCs are capable of production/storage of different product types (which is the case in real application), the job opportunities will increase accordingly. The third objective function is equal to Z 3 ¼ 875; 327, which according to the on hand cap-and-trade strategy and predefined threshold, is an acceptable amount for the whole chain. With applying node and flow complexity measures, the network’s structure has been simplified, in such a way that with the establishment of one manufacturing center, two main DCs, and five local DCs, the network can fully respond to all demands (U ¼ 0). It should also be mentioned that even though the direct shipment from manufacturing centers to local DCs is allowable in our model, the final result finds it optimal to ship the products mostly through the main DCs. This strategy significantly reduces excessive links throughout the network (=low flow complexity) as well as allows demand consolidation in main DCs, which ultimately leads to achieve highest demand satisfaction level.

6. Conclusion In this paper, we proposed a new mathematical model for designing a pharmaceutical supply chain network in a multiperiod planning horizon. High vulnerability of pharmaceutical industry to disasters (either internal or external), which ultimately results in production stoppage or failing to meet demand, along with arising consciousness toward environmental issues led to considering the chain as an integrated sustainable-resilient one. In doing so, new measures of resiliency and sustainability have been introduced and developed and further applied into the model. Minimization of total cost, maximization of overall social impact out of opening facilities, minimization of environmental and predefined resilience measures are the four targeted objective functions for the problem. In order to cope with the inherent uncertainty of the input data and their high variation over time, they have been considered under uncertainty and a new approach based on mean-absolute deviation and fuzzy possibilistic-stochastic programming. Furthermore, due to NP-hard nature of the problem, a novel Pareto-based lower bound as well as a new meta-heuristic algorithm called DVG has been developed. The results show the superiority of the developed algorithm over the two widely used ones in the literature (i.e., NSGA-II and MOICA).

136

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Several numerical examples along with a case study focusing on TruvadaÓ medicine on LGBTQ community in France have been provided to validate both efficiency of the introduced approaches and applicability of the model. Several sensitivity analyses are provided to address the behavior of objective functions with respect to variations of these data. Even though TruvadaÓ itself is not considered as a highly perishable product, perishability is becoming a rising issue in pharmaceutical industry. Accordingly, focusing on this issue is important as a future research in this field. Additionally, introducing other new approaches to cope with both uncertain nature of input parameters and complexity of the problem is recommended. To do so, one can consider designing a scenario tree and model it through stochastic programming and compare the results in terms of efficiency. For complexity, developing other meta-heuristic or exact algorithms is recommended. Most importantly, even though the resiliency has been considered throughout the network in this paper, the disasters only affect the production system (or equivalently manufacturing centers). It is recommended that disruption of DCs as well as routes between parties are taken into considerations. The analyses and results lead us to the following conclusions:
In the case where the fixed cost exceeds the variable cost (as in ours), increasing the number of DCs (either main or local) would result in increasing the total cost, since the reduction in transportation cost out of adding an extra warehouse is less than the corresponding establishment cost.
In the case where the variable cost is equal or greater than  115% of the fixed cost, adding an extra local DC decreases the total cost, as it provides more connection possibilities and flow options. For the values equal or greater than  153%, adding an extra main DC is recommended as it reduces the total cost of the network.
The capacity level of manufacturers and DCs play a significant role in network resiliency. Increase in the capacity levels will directly decrease the required number of manufacturers and/or DCs in the network and ultimately reduces the flow and node complexities (defined as resiliency measures). Furthermore, higher capacity levels translate to higher predefined threshold level in node criticality measure, which further lead to decrease in the number of corresponding nodes in the network.
As previously mentioned, the uncertainty level is set to 0.9. For the level of uncertainty 6 0:6, both the network configuration (in terms of locations and links) and the objective function values change, in such a way that Z 1 ; Z 2 ; Z 3 and Z 4 would have 18.1%, 9.3%, 12.4% and 7.2% decrease in value, respectively. However, when the uncertainty level varies in the interval of 0:6 6 n 6 0:9, only the objective function values change, while the network structure remains the same. These findings highlight the importance of considering uncertainty in our network as well as the vulnerability of the network structure, when the uncertainty level goes beyond 0.6.
Considering the established nodes, only one plant can satisfy all the existing demand throughout the country in our studied case. Nevertheless, as previously mentioned, this plant is equipped with different levels of backup technologies, in such a way that in case of any disruptions, the production process would switch to another level of technology, without any stoppage in production. Among the established local DCs, #35 is the least, and #26 and #11 are the most resilient local DCs according to the flow complexity measure. Accordingly, DC #35 is expected to be the most susceptible to disruption among all existing nodes. For established main DCs, DC #36 is the most and DC #15 is the least resilient main DCs.
As discussed in Section 3.1.2, in case of unavailability of any technology level resulted by external or internal issues, products are assigned to backup technologies. These reassignments are often costly, in such a way that in our case study, reassignment costs dominate other variable costs. Therefore, it is recommended to assign the products to the technology levels with the lower probability of failure for the cases similar to ours. This further avoids possible additional high cost of reassignment.

Acknowledgments This research was supported by the United States National Science Foundation (NSF) under award number 1334930. However, any opinions, findings, and conclusions or recommendations in this document are those of the authors and do not necessarily reflect views of the NSF. Appendix A. Linearization A.1. Linearization of constraint sets (47)–(49) Since the linearization structure of constraint sets (47)–(49) are quite similar to each other, we only go through constraint (47) and skip (48) and (49). Starting with inequality (47), we have three nonlinear terms qtilp  w0i and xtijpm  w0i , and x0tikpm  w0i . To convert them into the linear form, we introduce a nonnegative decision variable for each in such a way that t

0t

t

kilp ¼ qtilp  w0i , jijpm ¼ xtijpm  w0i , and jikpm ¼ x0tikpm  w0i , and following set of constraints:

kilp 6 Mw0i ; t

8i; l; p; t;

ðA:1Þ

137

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142 t

kilp 6 qtilp ;

8i; l; p; t;

ðA:2Þ

t

kilp P qtilp  Mð1  w0i Þ; jijpm 6 Mw0i ;

8i; l; p; t;

ðA:3Þ

t

8i; j; p; m; t;

ðA:4Þ

t

8i; j; p; m; t;

ðA:5Þ

jijpm 6 xtijpm ; t

jijpm P xtijpm  Mð1  w0i Þ;

8i; j; p; m; t;

ðA:6Þ

0t

8i; k; p; m; t;

ðA:7Þ

0t

8i; k; p; m; t;

ðA:8Þ

jikpm 6 Mw0i ; jikpm 6 x0tikpm ; 0t

jikpm P x0tikpm  Mð1  w0i Þ; t

t

0t

kilp ; jijpm ; jikpm P 0;

8i; k; p; m; t;

ðA:9Þ

8i; j; k; l; p; m; t;

ðA:10Þ

The linearized form of Eq. (47) can then be stated as:

ðU ti  eÞð1  w0i Þ þ

XX t XX X t X XX 0t kilp þ jijpm þ jikpm p2P l2L

j2J p2P m2M

m2M p2P k2K

XX XX X X XX P qtilp þ xtijpm þ x0tikpm ; p2P l2L

j2J p2P m2M

8i; t;

ðA:11Þ

m2M p2P k2K

The linearized forms of Eqs. (48) and (49) can be achieved using the same procedure. A.2. Linearization of Eq. (29) As probability Prtil0 p always lies between 0 and 1, with defining a new non-negative continuous variable r1;t bilp

¼ P r1;t  mr1;t and the following constraint sets, ilp ilp r1;t

6 mr1;t ilp ;

r1;t

P mr1;t þ Pr1;t  1; ilp ilp

r1;t

6 Pr1;t ilp ;

r1;t

P 0;

bilp bilp bilp bilp

8i; l 2 fL þ 1g; p; r > 1; t;

ðA:12Þ

8i; l 2 fL þ 1g; p; r > 1; t;

ðA:13Þ

8i; l 2 fL þ 1g; p; r > 1; t;

ðA:14Þ

8i; l 2 fL þ 1g; p; r > 1; t;

ðA:15Þ

the linear counterpart of Eq. (29) can be stated as: 0

Prtil0 p ¼ ð1  atil0 Þ

l X l

atil r1;t b ; 1  atil ilp

8i; l 2 fL þ 1g; p; r > 1; t:

ðA:16Þ

A.3. Linearization of Eq. (30) Since the Eq. includes a product of two continuous variables, we introduce two nonnegative continuous variables krtilp and 2

2

2

2

rt t rt rt rt 0rt rt 0rt rt 0rt t t k0rt ilp in such a way that P ilp qilp ¼ kilp  kilp , kilp ¼ 0:5ðP ilp þ qilp Þ and kilp ¼ 0:5ðP ilp  qilp Þ. Now the nonlinear terms kilp and kilp

should be linearized, whose linearization structure is discussed in detail in Appendix A.5. A.4. Linearization of Eq. (46) and sþt cannot take positive values at simultaneously. As it was mentioned before, Eq. (46) indicates that variables st i i Therefore, we define a new binary variable in such a way that:

138

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

( oti ¼

1 if st i ¼ 0 0

if sþt i ¼ 0

8i; t;

ðA:17Þ

Constraint (46) can then be transformed into its linear counterpart as follows; t sþt i 6 Moi ;

8i; t;

ðA:18Þ

t st i 6 Mð1  oi Þ;

8i; t:

ðA:19Þ

A.5. Linearization of Eq. (57) For linearizing Eq. (57), we should first define a binary variable ntil , where ntil ¼ ables in such a way that

n0t;t1 6 ntil ; il

n0t;t1 il

¼ ntil  nt1 and il

ðA:21Þ

8i; l; t;

6 Mn0t;t1 ; il

t;t1

6 qtilp ;

t;t1

P qtilp  ð1  xÞn0t;t1 ; il

t;t1

hilp

P 0;

ðA:22Þ

8i; l; p; t;

ðA:23Þ

8i; l; p; t;

2 f0; 1g; n0t;t1 il

Then, we define two new vari-

 qtilp , and the following constraint sets:

ðA:20Þ

t;t1

hilp

rt r2R milp .

8i; l; t;

P ntil þ nt1  1; n0t;t1 il il

hilp

¼

n0t;t1 il

8i; l; t;

6 nilt1 ; n0t;t1 il

hilp

t;t1 hilp

P

ðA:24Þ

8i; l; p; t;

ðA:25Þ

8i; l; t;

ðA:26Þ

8i; l; p; t;

ðA:27Þ

Eq. (57) would then be transformed as:

f(x)

c-1

3

2 1

1

2

3

x

Fig. 10. Piecewise linear estimation.

c

139

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

r 0tlp ðqtilp Þqtilp ¼

8 2 t;t1 t;t1 0t1 > ð1 þ dÞhilp þ r0tlpð1Þ ðqtilp ð1Þ  hilpð1Þ Þ  nlpð1Þ Q 1 qtilp ð1Þ Þ þ nlpð1Þ qtilp ð1Þ ðr lp > > ð1Þ > > > 2 t;t1 t;t1 > 0t1 < ðr lp ð1 þ dÞhilp þ r0tlp ðqtilp ð2Þ  hilp Þ  nlpð2Þ Q 2 qtilp ð2Þ Þ þ nlpð2Þ qtilp ð2Þ ð2Þ

if Q 1 6 qtilp ð1Þ 6 Q 2 if Q 2 6 qtilp ð2Þ 6 Q 3

ð2Þ

ð2Þ

> .. > > . > > > 0t1 2 > t;t1 t;t1 : ðr þ r 0tlpðn1Þ ðqtilpðn1Þ  hilpðn1Þ Þ  nlpðn1Þ Q n1 qtilpðn1Þ Þ þ nlpðn1Þ qtilpðn1Þ lpðn1Þ ð1 þ dÞhilp

if Q n1 6 qtilpðn1Þ 6 Q n ðA:28Þ

however, Eq. (A.28) is still nonlinear due to the term

qtilp2 .

For linearizing that, different linearization techniques can be uti-

lized. One of the mostly used and efficient approaches to deal with such nonlinearities is piecewise linear programming. Based on this approach, a non-convex function f ðxÞ can be linearized using Eq. (A.29) as follows:

f ðxÞ ¼ f ða1 Þ þ #1 ðx  a1 Þ þ

C1 X #c  #c1

2

c¼2

ðjx  ac j þ x  ac Þ:

ðA:29Þ

Table A1 Parameter setting used in the algorithms. DVG

Parameters Values

N 125

ItVNSout 20

ItVNSin 8

NFC max 20,000

G 1.5

w1 ; w2 ; w3 ; w4 0.25

MOICA

Parameters Values Parameters Values Parameters Values

n-Pop 300 b 2.15 n-Pop 300

N-imp 8 NFC max 20,000 PC 0.6

PA 0.64

PC 0.6

PR 0.32

n 0.125

Pm 0.32

NFC max 20,000

NSGA-II

Table A2 Random data generation range. Parameters

Problem size

Parameters

Small

Medium

Large

I J K H L R P M T ~f i ~f 0 j ~f 00

3 6 I 6 10 3 6 J 6 10 5 6 K 6 10 5 6 H 6 15 26L63 26R63 26P64 16M62 26T65 Uð500; 700Þ Uð40; 60Þ

15 6 I 6 25 15 6 J 6 20 15 6 K 6 25 20 6 H 6 30 56L67 56R67 56P69 36M64 10 6 T 6 20 Uð500; 700Þ Uð50; 80Þ

30 6 I 6 50 30 6 J 6 50 40 6 K 6 50 50 6 H 6 70 8 6 L 6 10 8 6 R 6 10 10 6 P 6 15 5 6 M 6 10 25 6 T 6 30 Uð600; 900Þ Uð50; 90Þ

f n00 a~ ae0 af00 af000 ~ b be0 f00 b

Uð10; 30Þ

Uð20; 30Þ

Uð20; 40Þ

f cq ilp a ~cijm a ~c0ikm a ~c00jkm a a ~c000 khm a f ch

Uð30; 50Þ Uð1; 5Þ Uð1; 5Þ Uð1; 5Þ Uð1; 5Þ

Uð30; 50Þ Uð1; 5Þ Uð1; 5Þ Uð1; 5Þ Uð1; 5Þ

Uð30; 50Þ Uð1; 5Þ Uð1; 5Þ Uð1; 5Þ Uð1; 5Þ

Uð10; 15Þ

Uð10; 15Þ

Uð10; 15Þ

f ch 0jp a f 00 a ch

Uð10; 15Þ

Uð10; 15Þ

Uð10; 15Þ

Uð10; 15Þ

Uð10; 15Þ

Uð10; 15Þ

cf v ia cf v 0i a f cp il f cp 0il ~j i ~j0 j ~j00

Uð2; 4Þ

Uð3; 5Þ

Uð4; 7Þ

Uð4; 7Þ Uð0:5; 1Þ Uð0:2; 0:4Þ Uð400; 900Þ Uð100; 150Þ

Uð5; 7Þ Uð0:5; 1Þ Uð0:2; 0:4Þ Uð400; 900Þ Uð100; 150Þ

Uð7; 10Þ Uð1; 2Þ Uð0:5; 0:9Þ Uð500; 1000Þ Uð150; 200Þ

Uð50; 150Þ Uð0:3; 0:6Þ Uð0:3; 0:6Þ Uð0:3; 0:6Þ Uð40; 60Þ Uð30; 50Þ

Uð50; 150Þ Uð0:3; 0:6Þ Uð0:3; 0:6Þ Uð0:3; 0:6Þ Uð40; 60Þ Uð30; 50Þ

Uð100; 150Þ Uð0:3; 0:6Þ Uð0:3; 0:6Þ Uð0:3; 0:6Þ Uð40; 80Þ Uð30; 50Þ

k

ip

kp

k

uri ur0j ur00k ~ n ne0 a b

Parameters are multiplied by 106 . Parameters are multiplied by 105 .

Problem size Small

Medium

Large

ef vi

Uð30; 50Þ Uð4; 7Þ Uð2; 5Þ Uð2; 5Þ Uð2; 5Þ Uð100; 150Þ Uð50; 70Þ Uð50; 70Þ Uð3; 9Þ Uð0:2; 0:8Þ Uð0:5; 1:5Þ

Uð30; 50Þ Uð4; 7Þ Uð2; 5Þ Uð2; 5Þ Uð2; 5Þ Uð100; 150Þ Uð50; 70Þ Uð50; 70Þ Uð3; 9Þ Uð0:2; 0:8Þ Uð0:5; 1:5Þ

Uð30; 50Þ Uð5; 9Þ Uð2; 5Þ Uð2; 5Þ Uð2; 5Þ Uð150; 170Þ Uð50; 100Þ Uð50; 100Þ Uð3; 11Þ Uð0:2; 0:8Þ Uð1; 2Þ

ef v 0j

Uð0:3; 0:4Þ

Uð0:3; 0:4Þ

Uð0:4; 0:6Þ

ef v 00k rdi 0 rdj 00 rdk ~r til

Uð0:1; 0:2Þ Uð0:4; 0:7Þ Uð0:4; 0:7Þ Uð0:4; 0:7Þ Uð5; 7Þ

Uð0:1; 0:2Þ Uð0:4; 0:7Þ Uð0:4; 0:7Þ Uð0:4; 0:7Þ Uð5; 7Þ

Uð0:2; 0:3Þ Uð0:4; 0:7Þ Uð0:4; 0:7Þ Uð0:4; 0:7Þ Uð10; 15Þ

r 0tilp ð:Þ e ei

Uð0:5; 1Þ

Uð0:5; 1Þ

Uð0:5; 1Þ

Uð0:1; 0:3Þ

Uð0:1; 0:3Þ

Uð0:1; 0:3Þ

e0 ei ikm e 00 ei

Uð0:1; 0:3Þ

Uð0:1; 0:3Þ

Uð0:1; 0:3Þ

jkm

Uð0:1; 0:3Þ

Uð0:1; 0:3Þ

Uð0:1; 0:3Þ

e 000 ei khm Capi b Cap0j b Cap00k b

gilp b

Uð0:1; 0:3Þ Uð0:1; 0:2Þ Uð3; 5Þ Uð1; 2Þ Uð4; 10Þ

Uð0:1; 0:3Þ Uð0:1; 0:2Þ Uð3; 5Þ Uð1; 2Þ Uð4; 10Þ

Uð0:1; 0:3Þ Uð0:1; 0:2Þ Uð3; 5Þ Uð1; 2Þ Uð4; 10Þ

~t b D hp

Uð1:5; 3:5Þ Uð0:1; 0:8Þ Uð0:2; 0:8Þ Uð500; 1500Þ Uð200; 300Þ Uð200; 300Þ

Uð1:5; 3:5Þ Uð0:1; 0:8Þ Uð0:2; 0:8Þ Uð500; 1500Þ Uð200; 300Þ Uð200; 300Þ

Uð2; 4Þ Uð0:1; 0:8Þ Uð0:2; 0:9Þ Uð500; 1500Þ Uð200; 300Þ Uð200; 300Þ

c~ s~

ijm

~ tl u

atil Gti Smax t Sþmax t

140

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

Table A3 The metric results out of medium scale instances. Problem no.

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

PM

RM

DM

SM

MID

CPU (s)

CM

DM

SM

MID

CPU (s)

0.64 0.707 – – – – – – – – – – – – – – – – – –

0.892 0.935 – – – – – – – – – – – – – – – – – –

0.738 0.666 – – – – – – – – – – – – – – – – – –

11276 19982 – – – – – – – – – – – – – – – – – –

1.435 1.493 – – – – – – – – – – – – – – – – – –

1.112 1.159 1.385 1.523 1.156 1.45 1.205 1.429 0.963 0.888 1.425 1.274 1.277 1.459 1.115 1.245 1.022 0.947 0.895 0.85

0.61 0.81 0.251 0.258 0.292 0.409 0.805 0.231 0.738 0.214 0.747 0.614 0.893 0.467 0.393 0.463 0.634 0.883 0.546 0.413

0.597 0.417 0.897 0.42 0.757 0.621 0.682 0.502 0.888 0.877 0.505 0.62 0.694 0.459 0.488 0.7 0.444 0.527 0.754 0.43

2009 2330 2678 2772 2808 2881 2952 2997 3112 3201 3299 3303 3320 3389 3409 3480 3505 3595 3661 3721

Table A4 The metric results out of large scale instances. Problem no.

PM

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

RM

DM

SM

MID

CPU (s)

CM

DM

SM

MID

CPU (s)

– – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – –

1.468 0.915 1.136 1.242 1.468 0.937 1.345 1.166 0.932 1.397 1.21 1.012 1.453 0.824 1.198 1.465 0.93 1.229 1.414 0.905

0.746 0.626 0.468 0.743 0.764 0.479 0.544 0.4 0.71 0.617 0.833 0.538 0.451 0.728 0.709 0.807 0.513 0.55 0.71 0.572

0.691 0.501 0.643 0.695 0.771 0.678 0.885 0.533 0.66 0.541 0.421 0.737 0.861 0.47 0.678 0.882 0.656 0.47 0.436 0.77

3882 4020 4329 4566 4625 4711 4857 4901 5082 5179 5350 5412 5512 5662 5790 6077 6619 7022 8224 9972

where f ðxÞ is the piecewise linear form of f ðxÞ, C is the number of divisions, and #c is the slope of the lines, which can be derived as follows (see Fig. 10):

#c ¼

f ðalþ1 Þ  f ðal Þ ; alþ1  al

ðA:30Þ

According to Chang (2002), Li and Yu (1999), Mirzapour Al-e-Hashem et al. (2013), if #c > #c1 , the linear counterpart for the absolute term could be written as:

f ðxÞ ¼ f ða1 Þ þ #1 ðx  a1 Þ þ

C1 X ð#c  #c1 Þðd1c þ x  ac Þ;

ðA:31Þ

c¼2

x  ac þ d1c P 0; d1c P 0;

8c:

8c;

ðA:32Þ ðA:33Þ

141

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142 Table A5 Values of case study input parameters.

a b

Parameters

Value/range

Parameters

Value/range

Parameters

Value/range

Parameters

Value/range

I J K H L R P M T ~f i ~f 0 j ~f 00 k f cq ilp a

48 48 48 48 5 5 1 3 24 [650, 700] [100, 120] [60, 70] [35, 40]

~cijm a ~c0ikm a ~c00jkm a a ~c000 khm f a ch ip f0 a ch jp f 00 a ch

[2, 3] [2, 3] [2, 3] [2, 3] [11, 12] [10, 11] [10, 11] [3, 4] [6, 7] [0.7, 0.9] [0.3, 0.6] [600, 700] [170, 190]

~j00 k uri ur0j ur00k ef vi ef v 0j ef v 00k rdi 0 rdj 00 rdk ~r til r 0tilp ð:Þ e ei ijm

[100, 120] [0.3, 0.4] [0.3, 0.4] [0.3, 0.4] [1. 1.7] [0.3, 0.35] [0.1, 0.15] [0.5, 0.6] [0.5, 0.6] [0.5, 0.6] [5, 7] [0.6, 0.9] [0.2, 0.3]

e0 ei ikm e 00 ei jkm e ei 000

[0.2, 0.3] [0.2, 0.3] [0.2, 0.3] [0.7, 1] [4, 5] [1, 1.5] [8, 9.5] [300, 350] [0.3, 0.35] [0.25, 0.4] [1200, 1300] [250, 270] [250, 270]

kp

cf v ia cf v 0i a f cp il f cp 0il ~j i ~j0 j

khm

Capi b Cap0j b Cap00k b

gilp b

~t b D hp ~ tl u t ail Gti Smax t Sþmax t

Parameters are multiplied by 106 . Parameters are multiplied by 105 .

Now, with the aid of Eqs. (A.31)–(A.33), qtilp2 can be transformed into its equivalent linear form as:

qtilp2 ¼ f ða1 Þ þ #1 ðqtilpðnÞ  a1 Þ þ

C1 X ð#c  #c1 Þðdt1ilpnc þ qtilpðnÞ  ac Þ;

ðA:34Þ

c¼2

qtilpðnÞ  ac þ dt1ilpnc P 0; dt1ilpnc P 0;

8i; l; p; n; c; t;

8i; l; p; n; c; t:

ðA:35Þ ðA:36Þ

See Tables A1–A5. References Aboolian, R., Cui, T., Shen, Z.-J.M., 2012. An efficient approach for solving reliable facility location models. INFORMS J. Comput. 25, 720–729. Alnaji, L., Ridha, M., 2013. The role of supply chain applications in jordanian pharmacies: a case study on pharmacies in the capital city Amman. Ind. Eng. Lett. 3, 65–71. Balaman, Sß.Y., Selim, H., 2016. Sustainable design of renewable energy supply chains integrated with district heating systems: a fuzzy optimization approach. J. Clean. Product. 133, 863–885. Behnamian, J., Zandieh, M., Ghomi, S.F., 2009. Parallel-machine scheduling problems with sequence-dependent setup times using an ACO, SA and VNS hybrid algorithm. Expert Syst. Appl. 36, 9637–9644. Benjaafar, S., Li, Y., Daskin, M., 2013. Carbon footprint and the management of supply chains: Insights from simple models. IEEE Trans. Autom. Sci. Eng. 10, 99–116. Brown, T., Gambhir, A., Florin, N., Fennell, P., 2012. Reducing CO2 emissions from heavy industry: a review of technologies and considerations for policy makers. Grantham Institute Briefing Paper 7. Calégari, P., Coray, G., Hertz, A., Kobler, D., Kuonen, P., 1999. A taxonomy of evolutionary algorithms in combinatorial optimization. J. Heurist. 5, 145–158. Cardoso, S.R., Barbosa-Póvoa, A.P., Relvas, S., Novais, A.Q., 2015. Resilience metrics in the assessment of complex supply-chains performance operating under demand uncertainty. Omega 56, 53–73. Carlsson, C., Fullér, R., 2001. On possibilistic mean value and variance of fuzzy numbers. Fuzzy Sets Syst. 122, 315–326. Carvalho, H., Barroso, A.P., Machado, V.H., Azevedo, S., Cruz-Machado, V., 2012. Supply chain redesign for resilience using simulation. Comput. Ind. Eng. 62, 329–341. Chang, C.-T., 2002. A modified goal programming model for piecewise linear functions. Eur. J. Oper. Res. 139, 62–67. Craighead, C.W., Blackhurst, J., Rungtusanatham, M.J., Handfield, R.B., 2007. The severity of supply chain disruptions: design characteristics and mitigation capabilities. Dec. Sci. 38, 131–156. Diabat, A., Battaïa, O., Nazzal, D., 2015. An improved lagrangian relaxation-based heuristic for a joint location-inventory problem. Comput. Oper. Res. 61, 170–178. Ding, H., Liu, Q., Zheng, L., 2016. Assessing the economic performance of an environmental sustainable supply chain in reducing environmental externalities. Eur. J. Oper. Res. 255, 463–480. Eskandarpour, M., Dejax, P., Miemczyk, J., Péton, O., 2015. Sustainable supply chain network design: an optimization-oriented review. Omega 54, 11–32. Fahimnia, B., Jabbarzadeh, A., 2016. Marrying supply chain sustainability and resilience: a match made in heaven. Transport. Res. Part E: Logist. Transport. Rev. 91, 306–324. Falasca, M., Zobel, C.W., Cook, D., 2008. ‘‘A decision support framework to assess supply chain resilience”. In: Proceedings of the 5th International ISCRAM Conference, pp. 596–605. Fiksel, J., 2003. Designing resilient, sustainable systems. Environ. Sci. Technol. 37, 5330–5339. Fiksel, J., 2006. Sustainability and resilience: toward a systems approach. Sustain.: Sci., Pract., Pol. 2. Grunow, M., Günther, H.-O., Yang, G., 2003. Plant co-ordination in pharmaceutics supply networks. OR Spectrum 25, 109–141. Hansen, K.R.N., Grunow, M., 2015. Planning operations before market launch for balancing time-to-market and risks in pharmaceutical supply chains. Int. J. Prod. Econ. 161, 129–139. Hansen, P., Mladenovic´, N., 2001. Variable neighborhood search: principles and applications. Eur. J. Oper. Res. 130, 449–467. (date visited: 03/05/2017). (date visited: 03/05/2017). (date visited: 03/05/2017). (date visited: 03/05/2017). , Date visited: 03/05/2017.

142

B. Zahiri et al. / Transportation Research Part E 103 (2017) 109–142

(date visited: 03/05/2017). (date visited: 03/05/2017). (date visited: 03/05/2017). Hua, G., Cheng, T., Wang, S., 2011. Managing carbon footprints in inventory management. Int. J. Prod. Econ. 132, 178–185. Jeihoonian, M., Zanjani, M.K., Gendreau, M., 2016. Accelerating Benders decomposition for closed-loop supply chain network design: case of used durable products with different quality levels. Eur. J. Oper. Res. 251, 830–845. Keyvanshokooh, E., Ryan, S.M., Kabir, E., 2016. Hybrid robust and stochastic optimization for closed-loop supply chain network design using accelerated Benders decomposition. Eur. J. Oper. Res. 249, 76–92. Klibi, W., Martel, A., 2012. Modeling approaches for the design of resilient supply networks under disruptions. Int. J. Prod. Econ. 135, 882–898. Laínez, J.M., Schaefer, E., Reklaitis, G.V., 2012. Challenges and opportunities in enterprise-wide optimization in the pharmaceutical industry. Comput. Chem. Eng. 47, 19–28. Levalle, R.R., Nof, S.Y., 2015a. Resilience by teaming in supply network formation and re-configuration. Int. J. Prod. Econ. 160, 80–93. Levalle, R.R., Nof, S.Y., 2015b. A resilience by teaming framework for collaborative supply networks. Comput. Ind. Eng. 90, 67–85. Levis, A.A., Papageorgiou, L.G., 2004. A hierarchical solution approach for multi-site capacity planning under uncertainty in the pharmaceutical industry. Comput. Chem. Eng. 28, 707–725. Li, H.-L., Yu, C.-S., 1999. A global optimization method for nonconvex separable programming problems. Eur. J. Oper. Res. 117, 275–292. Liu, L., Huang, G., Liu, Y., Fuller, G., Zeng, G., 2003. A fuzzy-stochastic robust programming model for regional air quality management under uncertainty. Eng. Optimiz. 35, 177–199. Masoumi, A.H., Yu, M., Nagurney, A., 2012. A supply chain generalized network oligopoly model for pharmaceuticals under brand differentiation and perishability. Transport. Res. Part E: Logist. Transport. Rev. 48, 762–780. Mavrotas, G., 2009. Effective implementation of the e-constraint method in multi-objective mathematical programming problems. Appl. Math. Comput. 213, 455–465. Mavrotas, G., Florios, K., 2013. An improved version of the augmented e-constraint method (AUGMECON2) for finding the exact pareto set in multi-objective integer programming problems. Appl. Math. Comput. 219, 9652–9669. Mirzapour Al-e-Hashem, S., Baboli, A., Sazvar, Z., 2013. A stochastic aggregate production planning model in a green supply chain: considering flexible lead times, nonlinear purchase and shortage cost functions. Eur. J. Oper. Res. 230, 26–41. Mohammadi, M., Jolai, F., Tavakkoli-Moghaddam, R., 2013. Solving a new stochastic multi-mode p-hub covering location problem considering risk by a novel multi-objective algorithm. Appl. Math. Model. 37, 10053–10073. Mohammadi, M., Torabi, S., Tavakkoli-Moghaddam, R., 2014. Sustainable hub location under mixed uncertainty. Transport. Res. Part E: Logist. Transport. Rev. 62, 89–115. Mohammadi, M., Dehbari, S., Vahdani, B., 2014. Design of a bi-objective reliable healthcare network with finite capacity queue under service covering uncertainty. Transport. Res. Part E: Logist. Transport. Rev. 72, 15–41. Mohammadi, M., Tavakkoli-Moghaddam, R., Siadat, A., Rahimi, Y., 2016. A game-based meta-heuristic for a fuzzy bi-objective reliable hub location problem. Eng. Appl. Artif. Intell. 50, 1–19. Mousazadeh, M., Torabi, S.A., Zahiri, B., 2015. Robust possibilistic programming approach for pharmaceutical supply chain network design (with a case study). Comput. Chem. Eng. 82, 115–128. Mousazadeh, M., Torabi, S., Zahiri, B., 2015. A robust possibilistic programming approach for pharmaceutical supply chain network design. Comput. Chem. Eng. 82, 115–128. Narayana, S.A., Pati, R.K., Vrat, P., 2012. Research on management issues in the pharmaceutical industry: a literature review. Int. J. Pharmaceut. Healthc. Market. 6, 351–375. Niakan, F., Vahdani, B., Mohammadi, M., 2015. A multi-objective optimization model for hub network design under uncertainty: an inexact rough-interval fuzzy approach. Eng. Optim. 47, 1670–1688. Pagell, M., Wu, Z., 2009. Building a more complete theory of sustainable supply chain management using case studies of 10 exemplars. J. Supply Chain Manage. 45, 37–56. Papageorgiou, L.G., 2009. Supply chain optimisation for the process industries: advances and opportunities. Comput. Chem. Eng. 33, 1931–1938. Papageorgiou, L.G., Rotstein, G.E., Shah, N., 2001. Strategic supply chain optimization for the pharmaceutical industries. Ind. Eng. Chem. Res. 40, 275–286. Pettit, T.J., Fiksel, J., Croxton, K.L., 2010. Ensuring supply chain resilience: development of a conceptual framework. J. Bus. Logist. 31, 1–21. Pop, P.C., Pintea, C.-M., Sitar, C.P., Hajdu-Ma˘celaru, M., 2015. An efficient reverse distribution system for solving sustainable supply chain network design problem. J. Appl. Logic 13, 105–113. Rossetti, C.L., Handfield, R., Dooley, K.J., 2011. Forces, trends, and decisions in pharmaceutical supply chain management. Int. J. Phys. Distrib. Logist. Manage. 41, 601–622. Sermet, C., Andrieu, V., Godman, M.B., Van Ganse, E., Haycox, A., Reynier, J.-P., 2010. Ongoing pharmaceutical reforms in France. Appl. Health Econ. Health Policy 8, 7–24. Seuring, S., 2013. A review of modeling approaches for sustainable supply chain management. Decis. Support Syst. 54, 1513–1520. Shah, N., 2004. Pharmaceutical supply chains: key issues and strategies for optimisation. Comput. Chem. Eng. 28, 929–941. Shaw, K., Irfan, M., Shankar, R., Yadav, S.S., 2016. Low carbon chance constrained supply chain network design problem: a benders decomposition based approach. Comput. Ind. Eng. 98, 483–497. Sousa, R.T., Shah, N., Papageorgiou, L.G., 2005. Global supply chain network optimisation for pharmaceuticals. Comp. Aided Chem. Eng. 20, 1189–1194. Sousa, R.T., Liu, S., Papageorgiou, L.G., Shah, N., 2011. Global supply chain planning for pharmaceuticals. Chem. Eng. Res. Des. 89, 2396–2409. Stavins, R.N., 2003. Experience with market-based environmental policy instruments. Handb. Environ. Econ. 1, 355–435. Storn, R., Price, K., 1997. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. J. Global Optim. 11, 341–359. Susarla, N., Karimi, I.A., 2012. Integrated supply chain planning for multinational pharmaceutical enterprises. Comput. Chem. Eng. 42, 168–177. Tong, L., Zhou, X., Miller, H.J., 2015. Transportation network design for maximizing space–time accessibility. Transport. Res. Part B: Methodol. 81, 555–576. Torabi, S., Baghersad, M., Mansouri, S., 2015. Resilient supplier selection and order allocation under operational and disruption risks. Transport. Res. Part E: Logist. Transport. Rev. 79, 22–48. Varsei, M., Polyakovskiy, S., 2017. Sustainable supply chain network design: a case of the wine industry in Australia. Omega 66, 236–247. Waters, D., 2011. Supply Chain Risk Management: Vulnerability and Resilience in Logistics. Kogan Page Publishers. Weraikat, D., Zanjani, M.K., Lehoux, N., 2016. Two-echelon pharmaceutical reverse supply chain coordination with customers incentives. Int. J. Prod. Econ. 176, 41–52. Weraikat, D., Zanjani, M.K., Lehoux, N., 2016. Coordinating a green reverse supply chain in pharmaceutical sector by negotiation. Comput. Ind. Eng. 93, 67–77. Yu, X., Li, C., Shi, Y., Yu, M., 2010. Pharmaceutical supply chain in China: current issues and implications for health system reform. Health Policy 97, 8–15. Zahiri, B., Tavakkoli-Moghaddam, R., Mohammadi, M., Jula, P., 2014. Multi-objective design of an organ transplant network under uncertainty. Transport. Res. Part E: Logist. Transport. Rev. 72, 101–124. Zhalechian, M., Tavakkoli-Moghaddam, R., Zahiri, B., Mohammadi, M., 2016. Sustainable design of a closed-loop location-routing-inventory supply chain network under mixed uncertainty. Transport. Res. Part E: Logist. Transport. Rev. 89, 182–214. Zhang, W., Reimann, M., 2014. A simple augmentede-constraint method for multi-objective mathematical integer programming problems. Eur. J. Oper. Res. 234, 15–24. Zhang, P., Zhang, W.-G., 2014. Multiperiod mean absolute deviation fuzzy portfolio selection model with risk control and cardinality constraints. Fuzzy Sets Syst. 255, 74–91.