A multi-stage stochastic program for the sustainable design of biofuel supply chain networks under biomass supply uncertainty and disruption risk: A real-life case study

Transportation Research Part E 118 (2018) 534–567

Contents lists available at ScienceDirect

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

A multi-stage stochastic program for the sustainable design of biofuel supply chain networks under biomass supply uncertainty and disruption risk: A real-life case study Mohammad Fattahia, Kannan Govindanb,

T

⁎

a

Department of Industrial Engineering and Management, Shahrood University of Technology, Shahrood, Iran Center for Sustainable Supply Chain Engineering, Department of Technology and Innovation, University of Southern Denmark, Campusvej 55, Odense, Denmark b

A R T IC LE I N F O

ABS TRA CT

Keywords: Biofuel supply chain Sustainable development Multi-stage stochastic programming Disruption risk Scenario tree construction Rolling horizon simulation

We focus on design and planning of a biofuel supply chain (SC) network from biomass to demand centers where biomass supply is stochastic and seasonal, and facilities’ capacity varies randomly because of possible disruptions. We propose a cost-efficient multi-stage stochastic program in which the greenhouse gas emissions are mitigated and the social impact of the SC is considered. A rolling horizon procedure is presented to implement and evaluate the stochastic model solution. Computational results on a real-life case study demonstrate the proposed optimization tools’ applicability as well as the effect of disruption risk and sustainability dimensions on biofuel SC planning.

1. Introduction In order to achieve future energy security and sustainability, countries all over the world are diversifying their energy portfolio and focusing on renewable energy (Sharma et al., 2013). Production of biofuels, as a renewable energy, has obtained the significant attention in both academia and practice over recent years. In 2010, the worldwide biofuel production increased 17% from 2009 and reached to 105 billion liters, and biofuels provided 2.7% of the world's fuels for road transportation (Shalaby, 2013). Biofuels can be derived from biomass including agricultural, domestic, commercial, and industrial wastes, agricultural crops, and naturally grown biomass (Sharma et al., 2013). Bioethanol and biodiesel as liquid biofuels can be used as a fuel for vehicles or additives to petroleum-based fuels. The first generation biofuels are produced from food crops such as corn, sugar beets, sugarcane, and any sugar or starch, and their production is not sustainable due to the competition with food resources (Awudu & Zhang, 2012). However, their conversion technologies are well-developed and commercialized (Parker et al., 2008). On the other hand, the second and third generation biofuels, advanced biofuels, are obtained from lignocellulosic biomass such as non-food energy crops, herbaceous crops, agricultural residues, forest residues, industrial wastes, and urban woody wastes (Awudu & Zhang, 2012). Although using many types of lignocellulosic biomass can reduce the used areas for food cultivation due to the land use changes, lignocellulosic biomass has higher per-acre ethanol yield, lower impact on the land use and agriculture, better efficiency in terms of life-cycle environmental performance, and more variety of resources in comparison with the feedstock of first generation biofuels (Chen & Fan, 2012; Farrell et al., 2006; Naik et al., 2010). Many countries have set mandates and targets for future biofuel production in response to increasing trend of energy consumption. For example, in accordance with the U.S. Environmental Protection Agency (2007), 36

⁎

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

https://doi.org/10.1016/j.tre.2018.08.008 Received 25 February 2018; Received in revised form 12 August 2018; Accepted 14 August 2018 1366-5545/ © 2018 Elsevier Ltd. All rights reserved.

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

billion gallons of biofuels should be produced by 2022 which about 15 billion gallons can be corn-based ethanol and the remaining 21 billion gallons will be from advanced biofuels. The success of a biofuel industry relies on the efficiency and reliability of its supply chain (SC) from feedstock sites to energy demand centers (Huang et al., 2014). Main important features of biofuel SCs include (Lim & Ouyang, 2016; Huang et al., 2014; Marufuzzaman & Ekşioğlu, 2016): (1) high transportation cost because of the bulky volume of biomass, (2) uncertainty and seasonality of biomass supply, and (3) vast geographically dispersed supply network. The goal of this paper is to design sustainable cellulosic biofuel SCs in which the feedstock supply is stochastic and highly seasonal, and facilities’ capacity varies randomly because of possible disruptions. To guarantee efficient and reliable practices through biofuel SC networks, and capture seasonal nature of the feedstock supply, the proposed model in this paper will integrate tactical and strategic planning decisions over a multi-period planning horizon. In accordance with Govindan et al. (2017), the uncertainty sources in designing SC networks include: (1) inherent uncertainties, and (2) man-made or natural disruptions. In biofuel SCs, one of the most important inherent uncertainties that results in operational SC risk is biomass supply (Awudu & Zhang, 2012). Further, the disruption risk caused by SC disruptions can affect functionality of SC’s elements either completely or partially for uncertain time duration (Tang, 2006). In this study, the proposed model will serve managers of biofuel companies to design their SCs under operational and disruption risk. Further, biofuel SC design with stochastic parameters and multi-period setting can result in a multi-stage stochastic program (MSSP). This paper is the first attempt to present an MSSP for planning of biofuel SCs seeking mitigation strategies against the feedstock supply uncertainty, feedstock seasonality, and disruption risk. To make implementable the planning decisions of the model’s optimal policy, we also contribute to the adjustment of strategic and tactical decisions in a rolling horizon framework. Sustainable development is one of the main concerns of countries all over the globe in which the social impacts and the environmental performance of their industries should be guaranteed in addition to the economic development (Eskandarpour et al., 2015; Kannan, 2018). Based on Awudu & Zhang (2012) and our literature survey, a few papers such as Hombach et al. (2018) and You et al. (2012) designed a sustainable biofuel SC network consistent with economic objective, environmental performance, and social responsibility impact. In this paper, to have a sustainable biofuel SC, a lower bound and upper bound are set for the social impact and environmental cost of the SC, respectively. Further, using the ε-constraint method, we show the existing trade-off between sustainability dimensions. To examine and justify the applicability of the proposed stochastic program, a real-life case study from Iran is studied. The aim of this paper is to determine planning decisions related to a cellulosic bioethanol SC that can be dispersed in Iran’s geographical area. Iran has a high potential for using biomass resources because of its good climate condition and waste lands for successful agricultural activities. Because of increasing trend in energy consumption in Iran, there exists a great tendency for using renewable energy. In this paper, to analyze and investigate the applicability, justification and commercial feasibility of bioethanol SCs in Iran, we have made a great effort in collecting, preparing, and aggregating large data sets. In summary, the key contributions of this study are as follows:

• The development of a novel MSSP that is linearized as a mixed integer linear programming model for integrated strategic and tactical sustainable planning of biofuel SC networks under operational and disruption risk as well as feedstock seasonality. • The development of a framework for scenario tree generation to model the corresponding stochastic parameters. • The development of a rolling horizon procedure to make planning decisions implementable and adjust strategic and tactical decisions over a multi-period planning horizon under uncertainty. • The development of a real-life case study based on data from Iran. The remainder of this paper is organized as follows: In Section 2, we briefly review the related literature. Characteristics of our optimization problem are explained in Section 3. In Section 4, the MSSP is represented. Scenario generation approach is explained in Section 5. Benders decomposition (BD) algorithm is developed for the presented MSSP as a solution approach in Section 6. The rolling horizon procedure including the simulation’s structure and corresponding formulation is described in Section 7. In Section 8, the reallife case study is explained in details. In Section 9, we present numerical implementations and sensitivity analyses. Derived managerial insights from this study are presented in Section 10. Finally, Section 11 contains conclusions, and future research directions. 2. Literature review Most of previous studies in the relevant area examined an individual part or process of biofuel production systems such as feedstock processing, transportation, and technology selection. However, in recent years, several research papers have dealt with planning of a biofuel SC as an integrated system from biomass resources to biofuel end users. In this section, a selective overview on biofuel SC planning is presented. Deterministic and single period models to minimize overall SC cost are presented by Ekşioğlu et al. (2010), Zamboni et al. (2009), and Elia et al. (2011). However, in order to capture time variable parameters and system dynamics in planning of biofuel industries, Lim & Ouyang (2016), Huang et al. (2010), Ekşioğlu et al. (2009), Xie et al (2014), and An et al. (2011) developed multi-period problems. Further, Hombach & Walther (2015) and Hombach et al. (2016) developed a multi-period model for design of SCs for second generation biofuels in which cultivation of biomass is considered in addition to the biofuel production system. These studies investigate the impact of European biofuel regulations on the design and performance of second generation biofuel SCs. Uncertainty in biofuel SC network design has been incorporated by several studies, and the major considered inherent 535

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

uncertainties include, but not limited to: (1) feedstock supply uncertainty (e.g., Huang et al., 2014; Chen & Fan 2012; Poudel et al., 2018), (2) demand uncertainty (e.g., Chen & Fan, 2012; Osmani & Zhang, 2013), (3) uncertainty of SC costs such as transportation and procurement cost (e.g., Osmani & Zhang, 2013; Dal-Mas et al., 2011), and (4) price uncertainty (e.g., Osmani & Zhang, 2013; DalMas et al., 2011). In addition to these inherent uncertainties in biofuel SCs, the disruption risk, rooted in man-made or natural disruptions, can affect the functionality of SC’s elements (Tang, 2006). A few studies such as Marufuzzaman & Ekşioğlu (2016) and Poudel et al. (2016) have addressed the disruption risk in planning of biofuel SC networks. In these studies, the probability failures are assumed for SC facilities or/and transportation links as pre-specified parameters. In this study, we have considered two types of uncertainty including the feedstock supply and disruption of SC facilities. Based on Govindan et al. (2017), decision-making environments under uncertainty can be categorized into three main groups as follows: 1. Decision-making environment with random parameters in which the probability distribution of stochastic parameters is known. Stochastic programming approaches, including two-stage stochastic programming, Multi-stage stochastic programming, and chance-constrained programming approach, belong to this group. 2. Decision-making environment with random parameters in which decision makers have not any information about the probability distribution of random parameters. Robust optimization models belong to this group. Several studies consider continuous uncertain parameters within pre-specified intervals, named as interval-uncertainty modeling, in this area. Other studies model discrete uncertain parameters using the scenario approach, and optimize the worst-case performance of the problem. 3. Fuzzy decision-making environment. Flexible and possibilistic programming are two well-known approaches to model ambiguity and vagueness under a fuzzy decision-making environment. Regarding robust optimization approach, Tong et al. (2014) proposed a robust optimization model for biofuel SC design in which the fuel demand and biomass supply varied within pre-defined intervals. Tay et al. (2013) proposed a scenario-based robust optimization problem that considers uncertainty in feedstock supply and demand. Hombach et al. (2018) presented a robust and multiobjective optimization framework with consideration of decision makers’ risk attitude, and applied their framework to the case study of German biodiesel SC. Using two-stage stochastic programming approach, several studies modeled the two-stage nature of decisions in designing biofuel SC networks. The studies that exploited this approach include: Dal-Mas et al. (2011), Chen & Fan (2012), Osmani & Zhang (2013), Kazemzadeh & Hu (2013), Li & Hu (2014), Huang et al. (2014), Gonela et al. (2015), Üster & Memişoğlu (2017), Quddus et al. (2017b), Poudel et al. (2017), and Quddus et al. (2017a). Among these studies, to capture both feedstock seasonality and uncertainty, several studies developed multi-period two-stage stochastic models (see e.g., Huang et al., 2014; Quddus et al., 2017b; Poudel et al., 2017; Quddus et al., 2017a). In the two-stage stochastic programs, it is assumed that there is a single moment for stochastic parameters to become realized (Govindan et al., 2017). However, the uncertainties have been realized progressively in more than one moment in real-life multi-period planning of biofuel SCs. For such a planning, extending MSSPs is more appropriate and it is one of the motivations for this study. Further, by using available historical data, biomass supply and disruption events are estimated as discrete scenarios, and hence a cost-efficient MSSP is proposed in this study. In scenario-based stochastic programs, a small set of scenarios might not correctly capture the existing uncertainty in the problem. To deal with this issue, the sample average approximation (SAA) approach has been used by Üster & Memişoğlu (2017), Quddus et al. (2017a), Quddus et al. (2017b), and Poudel et al. (2017). This study will exploit a scenario tree construction approach, developed by Heitsch & Romisch (2005) and Heitsch & Römisch (2009) to obtain an appropriate set of scenarios. Optimization problems for biofuel SC planning attain economic goals in terms of either profit maximization or cost minimization (Govindan et al., 2017). In the profit maximization, the SC’s profit is computed as revenues minus costs, and in the related literature, the cost minimization models are more popular than profit maximization ones. Further, SC’s costs usually include some components such as procurement cost, inventory cost, transportation cost, location cost and so on, which are different in various studies, and dependent on the SC planning decisions and structure. Besides economic goals, many biofuel industries consider the environmental performance and the social impact as other objectives to design a sustainable SC network. In the relevant literature, Elia et al. (2011), Giarola, et al. (2013), Babazadeh et al. (2017), You et al. (2012), Santibañez-Aguilar et al. (2014), Gonela et al. (2015), You & Wang (2011), Quddus et al. (2017b), Hombach & Walther (2015), Hombach et al. (2016), and Hombach et al. (2018) considered environmental goals/factors. In order to mitigate GHG life cycle emissions of biofuels, four main stages including cultivation of biomass, biofuel production, transportation, and usage can be considered. The biofuel production and transportation GHG emissions are taken into account by most of these studies such as our case. Santibañez-Aguilar et al. (2014), Gonela et al. (2015), Hombach & Walther (2015), Hombach et al. (2016), and Hombach et al. (2018) added the environmental impacts from the land use change to their life cycle analysis of biofuels. Finally, in the related literature, biomass cultivation and biofuel usage GHG emissions are considered by Hombach et al. (2018). In accordance with Jones (2010), CO2 is the primary and majority greenhouse gas, and SOx, NOx, and volatile organic compounds (VOCs) are other main emissions from biofuel industries that have major impacts on the air pollution or toxic (Seddighi & Ahmadi-Javid, 2015; Rabl et al., 2014). A few studies considered SOx, NOx, and VOC emissions in planning of biofuel SCs. In this study, by considering all these emission types, we calculate and limit the total environmental cost of the SC system such as Seddighi & Ahmadi-Javid (2015). A few studies addressed the social responsibility issues, as another dimension of sustainability, in biofuel SC planning. SantibañezAguilar et al. (2014) considered the job creation from biomass cultivation and Hombach et al. (2018) taken into account the competition with food/fodder industry and loss of biodiversity as social criteria. In this study, we derive the social responsibility criteria related to the establishment of biofuel production plants using social life cycle assessment (S-LCA) approach. The guidelines 536

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 1. A biofuel SC network.

for S-LCA approach are defined by UNEP/SETAC (2009) through the integration of social issues into the life cycle assessment, and the obtained social criteria are applied to calculate the social impact of the SC system. In terms of optimization aspect, most of research studies developed a multi-objective optimization problem for sustainable biofuel SC planning, and Gonela et al. (2015) limit the GHG life cycle emissions from the SC network via some constraints. Finally, the interested readers are referred to some survey studies that are recently published in this area: Awudu & Zhang (2012) categorized the related literature in terms of incorporating uncertainty and sustainability concepts. Sharma et al. (2013) and Shabani et al. (2013) presented a survey on biomass-biofuel SCs. These review papers analyzed several aspects in the existing studies such as SC decisions and structures, modeling approaches, and objective functions. Zandi Atashbar et al. (2017) investigated the literature on biomass SCs in terms of optimization aspects and modeling approaches. Based on the presented literature review, little work has been done to address: (1) the disruption risk, (2) the multi-period planning under uncertainty, and (3) the sustainability issues in biofuel SC planning and management. By simultaneously addressing all these aspects, a realistic MSSP is proposed in this study. Furthermore, the stochastic model is applied to investigate the economic feasibility and system effectiveness of converting cellulosic biomass including agricultural residues and forest biomass to ethanol in Iran. Finally, we contribute a framework for simulating and adjusting strategic and tactical planning decisions through a rolling horizon procedure. 3. Problem description In this section, the biofuel SC network structure, main assumptions, and planning decisions are described. The network includes feedstock suppliers, feedstock storage facilities, fuel production plants (bio-refineries), and demand centers. City gates are assumed as demand centers, and fuel dispensing to individual fuel stations of the cities is not considered. Fig. 1 shows the biofuel SC structure from cellulosic biomass resources to demand centers. A planning horizon divided into multiple periods is considered. Different types of feedstock should be forwarded from suppliers to storage facilities or biofuel production plants. It is worth noting that the seasonality of the feedstock supply causes temporal dimension of the problem. Feedstock storage facilities, as the second layer of the network, by maintaining inventories may impose an extra cost to the SC, but reduce shortage risk of future supply. As shown in Fig. 1, such as Xie et al (2014) and Marufuzzaman & Ekşioğlu (2016), a storage facility in the network serves as a shipment consolidation point, the first hub, and a production plant is a deconsolidation point as the second hub. Therefore, a multimodal transportation network is assumed between storages and plants to optimize the transportation cost and increase the network flexibility. In accordance with Marufuzzaman & Ekşioğlu (2016), various modes of transportation typically used in real-life biofuel SCs for high volume and long-haul transportation of bulk products. In this study, the multimodal transportation network includes rail and truck transportation. The transportation cost by rail is calculated as a function of traveled distance and a fixed transportation cost (Marufuzzaman & Ekşioğlu, 2016; Dal-Mas et al., 2011). In our case study, we have assumed that to use the rail transportation with a predetermined capacity at each facility and over a time period; a fixed cost is required. To obtain bioethanol from cellulosic biomass, an advanced biofuel conversion technology should be used in plants. Main 537

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

manufacturing processes at plants include pretreatment, fermentation, distillation, and solid recovery. Finally, produced bioethanol should be transported from plants to demand centers. It is not mandatory to fulfill the demand of customers and hence the penalty cost for per unit of unsatisfied demand that is equal to the price of a substitute product is considered. 3.1. Dealing with uncertainty Because of unpredictable weather conditions and changing regulations and policies, the amount of feedstock supply at various geographical areas has inherent uncertainty. The annual supply amount of feedstock fields in previous years can be available based on the historical information. However, annual yields of some feedstock types such as non-food energy crops, agricultural residues, and forest residues should be distributed to each time period in accordance with their harvesting season. We assume that each storage facility may be disrupted independently because of natural disasters. Therefore, the available capacity of each storage facility in each period is considered as a stochastic parameter (see e.g., Fattahi et al., 2017b; Peng et al., 2011; Mak & Shen, 2012). In principle, due to possible disruptions, the available capacity of storage facilities can be reduced in each period based on the magnitude and type of disruption events and the recovery time of disrupted facilities. Moreover, it is assumed that disrupted storages would be recovered after disruption events and their recovery costs can be obtained. A set of scenarios as a scenario tree is considered to model stochastic parameters including feedstock supply, storages’ capacity, and recovery cost. This study achieves the economic goal by minimizing the expected of SC costs. 3.2. Planning decisions In the problem, seven types of decisions have to be made: (i) location decisions, (ii) capacity decisions, (iii) activation/discontinued use of storage facilities, (iv) transportation mode selection, (v) storage level of feedstock and biofuel, (vi) shipment decisions, and (vii) production decisions. Location decisions should be determined in two layers of the network including storage facilities and biofuel production plants. Furthermore, the biofuel industry should determine the capacity of each established plant. However, storage facilities which are located by the biofuel company have predetermined capacity for handling (the capacity of receiving biomass) and holding biomass. Location and capacity decisions have to be made at the beginning of the planning horizon before uncertainty realization as strategic decisions that would be constant over the planning horizon. In each period, the company should decide about the activation of established storage facilities or discontinued use of them. In addition, various modes of transportation can be selected for feedstock shipment from storages to plants. A fixed transportation cost is required at storage facilities to use the rail transportation with a predetermined capacity over a time period. These planning decisions, decisions (iii) and (iv), are dynamic (time-dependent) and have to be made in each period before uncertainty realization. Finally, decisions (v), (vi), and (vii) would be made only after realization of uncertain parameters in each period. In order to hedge against supply uncertainty, feedstock inventory can be held in storage facilities. Feedstock and biofuel shipments through the network’s entities, production, and inventory level of biofuel at plants are other tactical planning decisions. In the network, the produced biofuel will be sent to demand centers or stocked at plants. Fig. 2 illustrates the decisions to be determined in a generic time period and over the planning horizon. 3.3. Environmental performance consideration The environmental performance of the SC is evaluated by the expected environmental cost due to different types of emission including CO2, NOx, SOx, and VOC. Unit environmental cost is considered for different types of emission such as Seddighi & AhmadiJavid (2015), and in each period, the amount of the SC expected environmental cost caused by processing, handling, and transportation of feedstock/biofuel is limited. The limit for the environmental cost can be obtained based on emissions’ trading that is government-mandated.

>ŽĐĂƟŽŶĂŶĚ ĂƉĂĐŝƚLJĚĞĐŝƐŝŽŶƐ

ĐƟǀĂƟŽŶͬĚŝƐĐŽŶƟŶƵĞĚ ƵƐĞŽĨƐƚŽƌĂŐĞƐ dƌĂŶƐƉŽƌƚĂƟŽŶŵŽĚĞ ƐĞůĞĐƟŽŶ

/ŶǀĞŶƚŽƌLJ͕ ƐŚŝƉŵĞŶƚ͕ ĂŶĚƉƌŽĚƵĐƟŽŶ

Time period Planning horizon

Beginning of the planning horizon Uncertainty realization Fig. 2. Planning decisions in each period and over the planning horizon. 538

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 3. The developed inventory indicators, impressionable stakeholders, and selected impact subcategories.

3.4. Social impacts consideration To assess social impacts of the SC network, the social scores of biofuel production plants in different potential locations and with various capacity levels are estimated by using the S-LCA methodology that is suggested by UNEP/SETAC (2009). Generally, the S-LCA methodology presents a framework for evaluation of social and socio-economic aspects of products or services through their complete life cycle in four main steps: 1- goal and scope definition, 2- determination of inventory indicators for the considered goal and scope, 3- assessment of social impacts based on the inventory indicators (social indicators in this study), and 4- interpretation of results. These steps are explained in details in Appendix A. In this study, we have assumed that the final product (bioethanol) and available technologies of different biofuel production plants are the same, and our objective is to estimate the social scores of potential geographical location and capacity level of plants. To achieve this objective, in accordance with the considered scope that includes the entire stages of the biofuel production, we have interviewed with a panel of experts, including managers and engineers in biofuel/energy companies and other stakeholders such as local communities, about all five stakeholder categories and 31 impact subcategories related to the S-LCA approach (see Appendix A). Fig. 3 shows the developed three social indicators, the considered impressionable stakeholders and selected impact subcategories. Consider c ∈ C , j ∈ J , and e ∈ E as an element of social indicators, potential locations for plants, and plants’ capacity levels set, respectively. In order to quantify the social scores of different potential locations and capacity levels of plants, parameter ρc, j, e is defined as the score of social indicator c for establishing biofuel production plant at location j with capacity level e. In step 3 of the SLCA methodology, fuzzy analytic hierarchy process (AHP) approach is used to obtain the quantitative values of this parameter based on linguistic preferences of the experts’ panel. The fuzzy AHP approach is also explained in details in Appendices B and C. Finally, the social score of the network would be guaranteed to be more than a predetermined target by imposing a constraint into the optimization problem. 4. Model formulation Based on the SC network’s assumptions, explained in the previous section, an MSSP is developed. In Table 1, the needed notations including sets, parameters, and decision variables are described. The complete mathematical formulation of the proposed model is as follows. 4.1. Objective function

Minimize :

O LCostP =

O DisCost (ω) + OUCostS (ω) + OCCostT (ω) + O ICostP (ω) + O ICostS (ω) ⎞ ⎫ ⎧ LCostP O + EΩ ⎜⎛ ⎟ , TCost (ω) + O ProcurementCost (ω) + O ProductionCost (ω) + O PenaltyCost (ω) ⎬ ⎨ ⎝+ O ⎠⎭ ⎩

∑ ∑ ηj, e Xj, e + ∑ γi Yi , j∈J e∈E

(1)

(1a)

i∈I

539

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Table 1 Description of the sets, parameters, and decision variables. Sets: M S (m) I J E K L (i) T C G Ω

Set Set Set Set Set Set Set Set Set Set Set

Costs ηj, e

Annualized fixed cost for opening biofuel production plant j with capacity level e,

of of of of of of of of of of of

feedstock (biomass) types, m ∈ M, feedstock fields (suppliers) of type m , s ∈ S (m), potential locations for feedstock storages, i ∈ I, potential locations for biofuel production plants, j ∈ J, capacity levels (small, medium, and large) for biofuel production plants, e ∈ E, demand centers (markets), k ∈ K, possible transportation modes at feedstock storage i , l ∈ L (i), time periods, t , t ′ ∈ T , social criteria, c ∈ C, emission types to the air, g ∈ G, scenarios, ω ∈ Ω .

fuit, l

Fixed cost for establishing feedstock storage at location i , Fixed cost for using transportation mode l at feedstock storage i in time period t,

γi gsit

Fixed recovery gain related to discontinued use of a feedstock storage at location i in time period t,

usit

Fixed cost for activating feedstock storage at location i in time period t,

dsit (ω) bcm pcj, m

Recovery cost related to disruption events at feedstock storage i in time period t under scenario ω,

ism ipj

Unit inventory cost ($/ton) of feedstock of type m during a time period, Unit inventory cost ($/gallon) of biofuel at production plant j during a time period,

uckt

Unit penalty cost ($/gallon) of unsatisfied biofuel at demand center k in period t,

τs, i τs, j

Unit flow cost of biomass ($/ton) in transportation path (s, i) using trucks, Unit flow cost of biomass ($/ton) in transportation path (s, j ) using trucks, Unit flow cost of biomass ($/ton) in transportation path (i, j ) using transportation mode l,

τi, j, l τ j, k αg

Unit procurement cost ($/ton) for feedstock of type m, Unit production cost ($/gallon) of biofuel at production plant j from feedstock type m,

Unit flow cost of biofuel ($/gallon) in transportation path (j, k ) using trucks, Unit environmental cost ($/kg emission) for emission of type g.

Other parameters t θm , s (ω)

Feedstock yield of type m at feedstock field s during time period t under scenario ω,

ψit (ω)

Available (non-disrupted) fraction of capacity for storage i in time period t under scenario ω,

cai chi βit, m

Capacity of storage facility i in holding feedstock at each time period, Handling capacity of storage facility i in receiving biomass from feedstock fields at each time period, Deterioration rate of feedstock of type m during time period t in feedstock storage i ,

cpj, e

The production capacity of fuel production plant j with capacity level e,

csj, e

The fuel storage capacity in production plant j with capacity level e,

dkt κm λm νm ctl, i cl l, i, j

Biofuel demand at center k in time period t,

φg, s, i (φg, s, j )

The amount of emission (kilogram/ton) of type g from transportation of biomass from feedstock field s to storage facility i (plant j) using trucks,

φl, g, i, j

The amount of emission (kilogram/ton) of type g from transportation of biomass through transportation link (i, j ) using transportation mode l,

φg, j, k

The amount of emission (kilogram/gallon) of type g from transportation of biofuel through transportation link (j, k ) using trucks,

δg, m

The amount of emission (kilogram/gallon) of type g from unit production of biofuel from feedstock of type m,

ρc, j, e

The score of social criterion c for establishing fuel production plant j with capacity level e.

Coefficient for using capacity of feedstock storages for holding each unit feedstock of type m, Conservation rate (gallon/ton) from feedstock of type m to biofuel, Moisture content of feedstock of type m, Capacity of transportation mode l at storage facility i in each time period, Capacity (ton) for transportation mode l between node i and j in each time period,

The importance weight of social criterion c. δc Decision variables Xj, e A binary variable that equals 1 if biofuel production plant is established at location j with capacity level e,

Zit (ω)

A binary variable that equals 1 if storage facility is established at location i, A binary variable that equals 1 if storage facility is active at location i in time period t under scenario ω,

Wit, l (ω)

A binary variable that equals 1 if transportation mode l is used at storage facility i in time period t under scenario ω,

Fst, i, m (ω)

A non-negative real-valued variable representing the amount of feedstock of type m transported through path (s, i) in time period t under scenario ω, A non-negative real-valued variable representing the amount of feedstock of type m transported through path (s, j ) in time period t under scenario ω, A non-negative real-valued variable representing the amount of feedstock of type m transported through path (i, j ) using transportation mode l in time period t under scenario ω,

Yi

Fst, j, m (ω) Fit, j, m, l (ω)

(continued on next page)

540

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Table 1 (continued) F tj, k (ω)

A non-negative real-valued variable representing the amount of biofuel transported through path (j, k ) in time period t under scenario ω,

Hit, m (ω)

A non-negative real-valued variable representing the amount of feedstock of type m held in storage facility i in time period t under scenario ω,

P tj, m (ω) V tj (ω)

A non-negative real-valued variable representing the amount of biofuel produced in production plant j from feedstock of type m in time period t under scenario ω, A non-negative real-valued variable representing the amount of biofuel held in production plant j in time period t under scenario ω,

Ukt (ω)

A non-negative real-valued variable representing the amount of unsatisfied demand at market k in time period t under scenario ω,

O DisCost (ω) =

∑ ∑ dsit (ω) Yi ,

(1b)

i∈I t∈T

OUCostS (ω) =

∑ ⎛⎜∑ usit Zit (ω)(1−Zit−1 (ω))− ∑ gsit Zit−1 (ω)(1−Zit (ω)) ⎞⎟, ⎝ i∈I

t∈T

OCCostT (ω) =

∑ ⎛⎜∑ ∑ fuit, l Wit, l (ω) ⎞⎟, ⎝ i∈I

t∈T

O ICostP (ω) =

⎛

⎞

⎝ j∈J

⎠

∑ ⎛⎜∑ ∑ ⎝ i∈I

t∈T

m∈M

(1e)

⎞ ism Hit, m (ω) ⎟, ⎠

(1f)

τs, i Fst, i, m (ω)

OTCost (ω) =

τs, j Fst, j, m (ω) ⎞

+ ∑ ∑ ∑ ⎛ ∑ ∑ ∑ m ∈ M s ∈ S (m) j ∈ J ⎜ m ∈ M s ∈ S (m) i ∈ I ⎟, t t + τ F ∑ ∑ ∑ ∑ i , j , l i, j, m, l (ω) + ∑ ∑ τ j, k F j, k (ω) ⎟ ⎜ i ∈ I j ∈ J m ∈ M l ∈ L (i) j∈J k∈K ⎝ ⎠

∑ t∈T

O ProcurementCost (ω) =

∑ ∑ t∈T m∈M

O ProductionCost (ω) =

⎛ ⎜bcm ⎝

∑∑ ∑

∑ s ∈ S (m)

⎛ Ft (ω) + ⎜∑ s, i, m ∈ i I ⎝

(1g)

⎞⎞

∑ Fst, j, m (ω) ⎟ ⎟, j∈J

⎠⎠

(1h)

pcj, m P tj, m (ω), (1i)

t∈T j∈J m∈M

∑ ∑ uckt Ukt (ω).

O PenaltyCost (ω) =

(1c)

(1d)

⎠

l∈L

∑ ⎜ ∑ ipj V tj (ω) ⎟,

t∈T

O ICostS (ω) =

⎠

i∈I

(1j)

t∈T k∈K

The objective function (1) minimizes the expected total SC cost. The capital cost associated with the establishment of biofuel production plants and feedstock storages is computed in Eq. (1a). Eq. (1b) calculates the cost of recovering feedstock storages after disruption events. The activation cost and gain of discontinued use of feedstock storages are computed by Eq. (1c). The fixed cost related to using transportation modes at storage facilities is computed by Eq. (1d). The tactical costs include fuel and feedstock inventory holding cost (Eqs. (1e) and (1f)), transportation cost (Eq. (1g)), feedstock procurement cost (Eq. (1h)), fuel production cost (Eq. (1i)), and penalty cost for unsatisfied demand (Eq. (1j)). 4.2. Constraints Constraints on feedstock sites:

∑ Fst, i, m (ω) + ∑ Fst, j, m (ω) ⩽ θmt ,s (ω) i∈I

∀ m ∈ M , ∀ s ∈ S (m), ∀ t ∈ T , ∀ ω ∈ Ω, (2)

j∈J

Constraints (2) ensure that total supplied feedstock at each field cannot exceed its availability. Constraints on feedstock storages:

Zit (ω) ⩽ Yi , Hit,+m1 (ω) =

∀ i ∈ I, ∀ t ∈ T, ∀ ω ∈ Ω ,

∑

Fst, i, m (ω)− ∑

s ∈ S (m)

Hit,+m1 (ω) = βit, m Hit, m (ω) +

∑

(3)

Fit, j, m, l (ω)

∀ i ∈ I , ∀ m ∈ M , t = 1, ∀ ω ∈ Ω, (4)

j ∈ J l ∈ L (i)

∑ s ∈ S (m)

Fst, i, m (ω)− ∑

∑

Fit, j, m, l (ω) ∀ i ∈ I , ∀ m ∈ M , ∀ t ∈ T {1}, ∀ ω ∈ Ω, (5)

j ∈ J l ∈ L (i)

541

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

∑

κm Hit, m (ω) ⩽ cai ψit (ω) Zit (ω),

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω, (6)

m∈M

∑ ∑

Fst, i, m (ω)

⩽

chi ψit

(ω) Zit

(ω).

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω, (7)

m ∈ M s ∈ S (m)

Under each scenario, constraints (3) guarantee that feedstock storage i can be active in period t, if the storage facility is established. Eqs. (4) and (5) compute the inventory level of feedstock of type m at the beginning of period t + 1 in facility i under scenario ω (Hit,+m1 (ω) ), which equals to the net input flow of feedstock in period t plus the inventory level at the beginning of period t that is reduced because of deterioration during the period. Constraints (6) and (7) ensure that each active storage facility cannot hold and handle feedstock more than its available capacity, successively. Further, due to the existing uncertainty in storages’ capacity, the right hand side of these constraints are multiplied by the amount of the available (non-disrupted) capacity fraction of the corresponding storage, ψit (ω) . Constraints on transportation modes:

Wit, l (ω) ⩽ Zit (ω),

∑ ∑

∀ i ∈ I , ∀ l ∈ L (i), ∀ t ∈ T , ∀ ω ∈ Ω ,

Fit, j, m, l (ω) ⩽ cti, l Wit, l (ω),

(8)

∀ i ∈ I , ∀ l ∈ L (i), ∀ t ∈ T , ∀ ω ∈ Ω, (9)

j∈J m∈M

∑

Fit, j, m, l (ω) ⩽ cli, j, l Wit, l (ω),

∀ i ∈ I , ∀ j ∈ J , ∀ l ∈ L (i), ∀ t ∈ T , ∀ ω ∈ Ω. (10)

m∈M

Constraints (8) indicate that available transportation modes, can be selected for active storage facility i in period t and under scenario ω . At each feedstock storage, constraints (9) assure that the total amount of feedstock forwarded using a transportation mode is limited to its capacity. Constraints (10) impose an upper bound limit for transportation links from a storage facility. The fixed cost is not considered for using the trucks, and there is not any limitation on their capacity. Constraints on biofuel production plants:

∑

P tj, m (ω) ⩽

m∈M

V tj (ω) ⩽

∑ cpj, e Xj, e ,

∀ j ∈ J , ∀ t ∈ T , ∀ ω ∈ Ω, (11)

e∈E

∑ csj, e Xj, e ,

∀ j ∈ J , ∀ t ∈ T , ∀ ω ∈ Ω, (12)

e∈E

P tj, m (ω) = λm (1−νm ) V tj + 1 (ω) =

∑

⎛ Ft (ω) + ⎜ ∑ s , j, m ∈ s S ( m ) ⎝

P tj, m (ω)−

∑ F tj, k (ω)

m∈M

V tj + 1 (ω) = V tj (ω) +

∑ ∑ i ∈ I l ∈ L (i)

⎞ Fit, j, m, l (ω) , ⎟ ⎠

∀ j ∈ J , ∀ t ∈ T , ∀ m ∈ M , ∀ ω ∈ Ω,

∀ j ∈ J , t = 1, ∀ ω ∈ Ω, (14)

k∈K

∑

P tj, m (ω)−

m∈M

(13)

∑ F tj, k (ω)

∀ j ∈ J , ∀ t ∈ T {1}, ∀ ω ∈ Ω, (15)

k∈K

Constraints (11) and (12) impose limitations on the biofuel production and storage capacity, respectively. Eqs. (13) ensure that the amount of biofuel production from feedstock type m equals to the available feedstock in plant j during period t and under scenario ω , which converted to the biofuel. Eqs. (14) and (15) compute the inventory level of biofuel at the beginning of period t + 1 in plant j under scenario ω (V tj + 1 (ω) ). Constraints on satisfying the demands:

∑ F tj, k (ω) + Ukt (ω) = dkt

∀ k ∈ K , ∀ t ∈ T , ∀ ω ∈ Ω, (16)

j∈J

Eqs. (16) ensure that during each time period and under each scenario, demand of each market should be satisfied either partially or completely. Constraints on the SC environmental impact: t t t ⎡ ⎛ ∑ ∑ ∑ φg, s, i Fs, i, m (ω) + ∑ ∑ ∑ φg, s, j Fs, j, m (ω) + ∑ ∑ ∑ ∑ φl, g, i, j Fi, j, m, l (ω) ⎞ ⎤ i ∈ I j ∈ J m ∈ M l ∈ L (i) m ∈ M s ∈ S (m) j ∈ J m ∈ M s ∈ S (m) i ∈ I ⎢ ⎟ ⎥ ⩽ UBEnvCost , EΩ ⎢ ∑ αg ⎜ ⎥ + ∑ ∑ φg, j, k F tj, k (ω) + ∑ ∑ δg, m P tj, m (ω) ⎟⎥ ⎜ g ∈ G ⎢ j∈J m∈M ⎝ j∈J k∈K ⎠⎦ ⎣

∀ t ∈ T.

(17) Constraints (17) guarantee that in each period, the expected of total environmental cost caused by biofuel production and transportation of feedstock and biofuel is restricted by a predetermined upper bound value.

542

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Constraint on the SC social responsibility score:

⎛

⎞

∑ ∂c ⎜ ∑ ∑ ρc,j,e Xj, e ⎟ ⩾ LB SocScore

c∈C

⎝ j∈J

e∈E

(18)

⎠

The social dimension of the sustainability is taken into account by Constraint (18). This constraint ensures that the social score related to the establishment of biofuel production plants should be greater than or equal a constant value. Non-anticipativity constraints: In MSSPs, a policy should be non-anticipative, which means the decisions determined at each stage cannot be dependent on the future realization of stochastic parameters. In the stochastic problem, decisions Z and W should be determined before knowing the uncertainty realization of stochastic parameters at each period. Here, the realization of stochastic parameters in period t and under scenario ω is presented by ξtω , ξtω = {θt (ω), ψt (ω)}. The uncertainty realization from the beginning of the planning horizon until period t under scenario ω is also presented by (ξ1ω, ξ2ω, ...,ξtω) . Eqs. (19)–(22) assure that these planning decisions are non-anticipative (Fattahi et al., 2017b; Fattahi et al. 2017a).

Zit (ω) = Zit (ω′) Wit, l (ω)

=

∀ i ∈ I , t = 1 , ∀ (ω, ω′) ∈ Ω,

Wit, l (ω′),

(19)

∀ i ∈ I , ∀ l ∈ L (i), t = 1 , ∀ (ω, ω′) ∈ Ω,

(20) ′

′

′

∀ i ∈ I , ∀ t ∈ T {1}, ∀ (ω, ω′) ∈ Ω: (ξ1ω, ξ2ω, ...,ξtω− 1) = (ξ1ω , ξ2ω , ...,ξtω− 1),

Zit (ω) = Zit (ω′), Wit, l (ω) = Wit, l (ω′),

∀ i ∈ I , ∀ l ∈ L (i), ∀ t ∈ T {1} , ∀ (ω, ω′) ∈ Ω: (ξ1ω, ξ2ω, ...,ξtω− 1) =

′ (ξ1ω ,

(21) ′ ξ2ω ,

′ ...,ξtω− 1).

(22)

Integer and non-negativity constraints:

F, H, P, V, U ⩾ 0

(23)

X , Y, Z , W ∈ {0, 1}

(24)

For convenience, the indices of variables are omitted in constraints (23) and (24). 4.3. Linearization of the objective function’s formulation Eq. (1c) is nonlinear because of the Zit (ω)(1−Zit − 1 (ω)) and Zit − 1 (ω)(1−Zit (ω)) mathematical expressions in which two binary decision variables are multiplied. In order to linearize Eq. (1c), a technique, proposed by Ghaderi et al. (2012), is used. In this t technique, binary variables Z¯it (ω) and Zi ̂ (ω) are defined, and Eqs. (25) and (26) are embedded into the optimization model as follows: t Zit − 1 (ω) + Z¯it (ω) = Zit (ω) + Zi ̂ (ω), t Zi ̂ (ω) + Z¯it (ω) ⩽ 1,

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω,

(25)

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω,

(26)

By using Eqs. (25) and (26), Z¯it (ω) equals 1 if non-active storage facility at location i is activated in time period t under scenario ω, t and Zi ̂ (ω) equals 1 if active storage facility at location i is discontinued use in time period t under scenario ω. Therefore, Eq. (1c) can be written as follows:

OOCostS (ω) =

∑ ⎛⎜∑ usit Z¯it (ω)− ∑ gsit Zi ̂ (ω) ⎞⎟. t

t∈T

⎝ i∈I

i∈I

(1k)

⎠

Further, Eqs. (27)–(30) can also be embedded into the problem that are resulted from relations (19), (21), (25), and (26).

Z¯it (ω) = Z¯it (ω′)

∀ i ∈ I , t = 1 , ∀ (ω, ω′) ∈ Ω,

(27)

t Zi ̂

∀ i ∈ I , t = 1 , ∀ (ω, ω′) ∈ Ω,

(28)

(ω) =

t Zi ̂

(ω′)

′ ′ ′ Z¯it (ω) = Z¯it (ω′) ∀ i ∈ I , ∀ t ∈ T {1} , ∀ (ω, ω′) ∈ Ω: (ξ1ω, ξ2ω, ...,ξtω− 1) = (ξ1ω , ξ2ω , ...,ξtω− 1),

t

t

Zi ̂ (ω) = Zi ̂ (ω′) ∀ i ∈ I , ∀ t ∈ T {1} , ∀ (ω, ω′) ∈ Ω: (ξ1ω, ξ2ω, ...,ξtω− 1) = (ξ1ω , ξ2ω , ...,ξtω− 1). ′

′

′

(29) (30)

5. Scenario tree construction for multivariate stochastic parameters One of the main challenges in MSSPs is to construct an efficient scenario tree which represents stochastic parameters, discretely. The scenario construction approach that is exploited in this paper has two main steps. In the first step, named as scenario generation step, discrete scenarios should be generated using the Monte-Carlo (MC) simulation as a scenario fan. In this step, some stochastic

543

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

processes are applied to model stochastic parameters over the planning horizon. In the second step, the generated scenario fan should be reduced and converted to a scenario tree. Stochastic feedstock supply: It is assumed that the supply of each type of feedstock in each geographical area follows a stochastic process. In order to handle feedstock supply fluctuations due to the seasonality, a modified autoregressive (AR) model of first order is applied. The model is proposed by Sodhi (2005) and then applied by several studies in different contexts (see e.g., Fattahi et al., 2017b; Fattahi, 2017; Fattahi & Govindan, 2017). In this study, seasons are considered as time periods and under each scenario ω and period t, supply amount of feedstock of type m at site s obtains by following relation: t t−1 θmt , s (ω) = θ¯m, s + Am, s (θmt −, s1 (ω)−θ¯m, s ) + εmt , s (ω),

(31)

where θ¯ is a period-dependent parameter representing the expected supply of different types of feedstock at each field and each period. Using this parameter, we can capture the corresponding temporal dimension comes from the seasonality of feedstock supply. Further, in relation (31), A is the autoregressive parameter and ε is the error term parameter following the normal distribution with mean zero and variance σε2 . Stochastic non-disrupted fraction of storage facilities’ capacity in face of disruptions: In this study, it is assumed that the disruption occurrence at each storage facility and in each period does not affect the probability of disruption in other geographical areas and subsequent periods. As pointed out by Snyder et al. (2016), for such a problem setting, disruptions can be considered as a special case of yield uncertainty in which the yield is a Bernoulli random variable. Therefore, we assume that probability of disruption at each storage facility and in each period follows a Bernoulli distribution with parameter 0 < pri, t < 1, which represents the disruption probability. Here, parameter Λ ti (ω) is a binary indicator which equals one if a disruption is occurred at facility i in period t and under scenario ω. Moreover, if a disruption event occurs at established storage facility i in period t and under scenario ω , we assume that the percentage of handling and holding capacity of the storage which would be disrupted over the time period, oit (ω) , is dependent on the recovery time of the storage and the disruption magnitude, and hence follows a normal distribution N (μi , σi ) . In this study, by assuming that a disrupted storage in each time period would be recovered over the same period, we use relation (32), to simulate parameter ψit (ω) .

o t (ω) ⎞ ψit (ω) = 1−Λ ti (ω) ⎛ i ⎝ 100 ⎠ ⎜

⎟

(32)

Finally, for each scenario ω , the amount of the capacity which should be recovered for storage i in each period can be computed after simulating parameter ψ . Therefore, in the biofuel SC planning, obtaining the value of parameter dsit (ω) can be simple by considering the recovery cost for each unit of capacity.

5.1. Scenario generation The MC simulation is applied to generate a set of discrete scenarios as a scenario fan. The discretization procedure for continuous probability density function of stochastic parameters via a set of discrete scenarios is called scenario generation. In relation (31), the error terms are time independent and are generated discretely based on their known distributions. Further, for the stochastic capacity of storage facilities due to the disruption events, the parameters Λ and o in relation (32) are randomly generated. In this step, the MC simulation generates a scenario fan and in the following sub-section, we will explain how to reduce the number of generated scenarios and convert the generated scenario fan into a scenario tree.

5.2. Scenario tree construction To avoid from computationally intractable stochastic programs, it is essential to efficiently reduce the number of generated scenarios. For this purpose, Dupačová et al. (2003) proposed two heuristics (forward and backward reduction) algorithms that are briefly explained in Appendix D. Generated scenarios for stochastic parameters using MC simulation are in the form of a fan of scenarios. In order to covert a scenario fan to a scenario tree and reduce the number of generated scenarios, a scenario tree construction method, proposed by Heitsch & Romisch (2005) and Heitsch & Römisch (2009), is exploited. The scenario tree construction methods are also divided into forward and backward construction procedures which are based on applying forward and backward reduction algorithms, respectively. The brief explanation regarding these methods is presented in Appendix E.

6. Benders decomposition In order to solve the proposed MSSP for a large number of scenarios in a reasonable time, an exact solution algorithm based on BD is developed for the problem. Benders (1962) presented BD algorithm to solve large-scaled linear programming (LP) and MILP problems with a specific structure, and this algorithm is widely applied for solving SC network design problems (See e.g., Santoso et al. 2005; Lin & Wang 2011; Fattahi et al. 2017a). A modified version of BD presented by Conejo et al. (2006) is applied in this study. Consider an optimization problem (OP) with a finite optimal solution as following: 544

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan T T ⎧ Min: z = c x + b y ⎪ s. t . Ax ⩾ d, OP: ⎨ Bx + Dy ⩾ h, ⎪ x ∈ X, y ⩾ 0, ⎩

where c ∈ Rn1, b ∈ Rn2 , d ∈ Rm1, h ∈ Rm2 , A, B, and D are m1 × n1, m2 × n1, and m2 × n2 , matrices, respectively, and X is a non-empty set of feasible combinations for decision vector x. In the OP, decision vector variables are partitioned into two sets x and y. According to BD algorithm, by considering decision vector x as a set of complicated decisions, this problem can be decomposed into a master problem (MP) and sub-problem (SP). At the first iteration of BD algorithm, the MP { 1} can be presented as follows: T ⎧ Min: c x + φ ⎪ s. t . φ LB ⩽ φ , MP { 1} : ⎨ Ax ⩾ d, ⎪ x ∈ X, φ ∈ R, ⎩

If we assume (x (1) , φ(1) ) as the optimal solution of MP { 1} , its objective function value, cT x (1) + φ(1) , is a lower bound for the OP. Further, the corresponding SP of BD algorithm at iteration κ is: T ⎧ Min: b y, ⎪ s. t . x = x (κ ) , SP { κ} : ⎨ Dy ⩾ h−Bx , ⎪ x ⩾ 0, y ⩾ 0, ⎩

In the optimization model SP { κ} , x (κ ) is the optimal solution for decision vector x from solving the MP at iteration κ . Consider y ( κ ) be the optimal solution of SP { κ} , then UB { κ} = cT x (κ ) + bT y ( κ ) is an upper-bound for the OP at iteration κ of BD algorithm. If Benders’ sub-problem be feasible for every solution of the feasible region of MP { 1} and λ(κ ) be the vector of optimal dual variables related to constraint x = x (κ ) from solving the sub-problem at iteration κ , an optimality cut of BD would be added to the MP at the next iteration of the algorithm as follows:

bT y (κ ) + (λ(κ ) )T (x −x (κ ) ) ⩽ φ

(33) UB ( κ ) − LB (κ )

⩽ ε , where ε is a constant value near to zero, and Benders’ solution for the OP is The algorithm converges whenever UB ( κ ) based on the best upper-bound of the algorithm until its termination. If the algorithm doesn’t terminate at iteration κ , we need to solve MP { κ + 1} at iteration κ + 1 as follows: T ⎧ Min: c x + φ ⎪ s. t . φ LB ⩽ φ , ⎪ MP { κ + 1} : Ax ⩾ d, ⎨ T (τ ) ⎪ b y + (λ(τ ) )T (x −x (τ ) ) ⩽ φ , τ = 1, 2, ..., κ ⎪ x ∈ X, φ ∈ R, ⎩

Afterward, BD algorithm is continued by solving SP { κ + 1} , and this procedure will be continued until the algorithm’s convergence. Our MSSP is decomposed into an MP and SP, and the MP includes binary decision variables, X , Y , Z , and W . The MP { κ + 1} and SP { κ} of BD algorithm are presented in Appendix F. In the biofuel SC network, it is not mandatory to fulfill the customers’ demand and in fact, it is possible to have stock-out for customers. Therefore, solving the SP after fixing binary variables X , Z , and W based on any solution from the feasible region of MP { 1} results in a feasible solution. Hence, to achieve BD algorithm’s convergence, only an optimality cut should be added to the MP at each iteration of the algorithm. At iteration κ of BD algorithm, let λ {X} (j,κe) , λ {Z} (i,κt), ω, and λ {W} i(,κl), t , ω be the optimal dual variables corresponding to constraints of (κ ) (κ ) , OFMP and φ(κ ) be the optimal objective of the SP { κ} , MP { κ} , and the optimal fixing decision vectors X , Z , and W in SP { κ} , and OFSP { κ } value of φ in MP , respectively. As a consequence, the upper and lower bound of the original problem, at iteration κ of BD algorithm (κ ) (κ ) (κ ) + OFMP −φ(κ ) and OFMP , respectively, and the optimality cut of BD algorithm is as follows (See Constraints can be computed as OFSP (F.30) in Appendix F): (κ ) OFSP + ∑ ∑ (λ {X} j(,κe) )(Xj, e −X (jκ, e) ) + ∑ ∑ ∑ (λ {Z} i(,κt), ω)(Zit (ω)−Zti,,ω(k ) )+ j∈J e∈E

∑ ∑ ∑ ∑

i∈I t∈T ω∈Ω

(λ {W} (i,κl), t , ω)(Wit, l (ω)−Wti,,l(,kω) )

⩽ φ,

i∈I t∈T l∈L ω∈Ω

(34)

7. Rolling horizon simulation A policy obtained from MSSPs is dependent on a scenario tree that is used to represent stochastic parameters. To make the policy implementable and test the solution empirically, the rolling horizon approach can be exploited. The rolling horizon approach allows decision makers to make use of data that is revealed as time progresses. The purpose of the rolling horizon simulation is to examine the outcomes of implementing the optimal solution over a planning 545

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 4. The pseudo code of the rolling horizon simulation.

horizon (Chand et al., 2002). Here, one generated path for |T| time periods, illustrated by ω , is considered as realized stochastic parameters. For implementation of a policy for period t, we solve the MSSP with a planning horizon including periods t, t + 1, …, |T| in which the optimal policy at period t is fixed and stochastic parameters for period t are known based on the realized path ω . It should be noted that the given optimal decisions that should be made before uncertainty realization at period t are considered as the optimal policy for period t, and obtained from solving the MSSP including |T|−t + 1 periods. At the beginning of the planning horizon, we solve the multi-stage stochastic programming model with a planning horizon of |T| periods. Therefore, by implementing the optimal policy from solving the model for the first period and realized path ω , the first period objective function obtains as follows: DisCost UCostS CCostT ICostP ICostS ⎧ LCostP ⎛ Ot = 1 (ω) + Ot = 1 (ω) + Ot = 1 (ω) + Ot = 1 (ω) + Ot = 1 (ω) ⎞ ⎫ , O +⎜ PenaltyCost TCost ProcurementCost ProductionCost ⎨ (ω) + Ot = 1 (ω) + Ot = 1 (ω) ⎟⎠ ⎬ ⎝+ Ot = 1 (ω) + Ot = 1 ⎭ ⎩

Then, for the next period, the stochastic model at the given current state of the problem, called rolling horizon model, is solved with a planning horizon of |T|−1 periods. The state of the problem includes predetermined decisions from implementing the optimal policy in the previous periods that influence on the current optimization model. The obtained optimal policy is again implemented and the objective function for period t can be obtained as follows: DisCost (ω) + OtUCostS (ω) + OtCCostT (ω) + OtICostP (ω) + OtICostS (ω) ⎞ ⎫ ⎧ ⎛ Ot ⎜ ⎨ + OtTCost (ω) + OtProcurementCost (ω) + OtProductionCost (ω) + OtPenaltyCost (ω) ⎟ ⎬ ⎠⎭ ⎩⎝

This process is repeated until we know the final state of the problem at the end of period |T|. The pseudo code for the rolling horizon simulation is explained in Fig. 4. In this study, we simulate the implementation of our model within a rolling horizon framework and compare the objective function with the model’s objective. 7.1. Rolling horizon model The rolling horizon model at period t, t > 1, is the same as the stochastic programming model with some modifications. In Table 2, notations for the rolling horizon model are presented and then, the modified constraints are formulated. In order to derive the rolling horizon formulation from the stochastic programming model, sets T and Ω are substituted by Th and Ωh , respectively. Further, constraints (3), (11), and (12) have been modified as follows:

Zit (ω) ⩽ Yi∗,

∀ i ∈ I , ∀ ω ∈ Ωh ,

(3-h)

Table 2 Notations for the rolling horizon model. Sets and parameters Th Set of time periods that the rolling horizon model will be solved for them, Set of scenarios that the rolling horizon model will be solved for them, Ωh t0 The virtual time related to the beginning of the planning horizon in the rolling horizon model, Binary indicator representing the fixed decisions related to the location and size of biofuel plants that are obtained from solving the stochastic model, X∗j, e

Y ∗i

Binary indicator representing fixed decisions related to the location of storage facilities that are obtained from solving the stochastic model,

Zti 0

Binary indicator that equals 1 if a storage facility i had been active in the previous period of the current time,

Hti,0m

Amount of feedstock of type m held in storage facility i at the beginning of the first period in the rolling horizon model based on implementation of the optimal policy for a particular scenario path in the previous time period, Amount of biofuel held in production plant j at the beginning of the first period in the rolling horizon model based on implementation of the optimal policy for a particular scenario path in the previous time period.

V tj0

546

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 5. The structure of the rolling horizon simulation for a planning horizon including three periods.

∑

P tj, m (ω) ⩽

m∈M

V tj (ω) ⩽

∑ cpj, e X∗j, e,

∀ j ∈ J , ∀ t ∈ T , ∀ ω ∈ Ωh ,

(11-h)

e∈E

∑ csj, e X∗j, e,

∀ j ∈ J , ∀ t ∈ T , ∀ ω ∈ Ωh ,

(12-h)

e∈E

Finally, for the first period of the rolling horizon model, constraints (4), (14), and (25) have been written as follows:

Hit,+m1 (ω) = βit, m Hti,0m +

∑

Fst, i, m (ω)− ∑

s ∈ S (m)

V tj + 1 (ω)

=

V tj0

+

∑

P tj, m (ω)−

m∈M t Zti 0 + Z¯it (ω) = Zit (ω) + Zi ̂ (ω),

∑

Fit, j, m, l (ω)

∀ i ∈ I , ∀ m ∈ M , t = 1, ∀ ω ∈ Ωh , (4-h)

j ∈ J l ∈ L (i)

∑

F tj, k (ω)

∀ j ∈ J , t = 1, ∀ ω ∈ Ωh ,

k∈K

∀ i ∈ I , t = 1, ∀ ω ∈ Ωh ,

(14-h) (25-h)

Fig. 5 shows the structure of the rolling horizon approach for decision making using the MSSP over a planning horizon with three periods. 8. Case study In this section, a real-life case study is developed using data from Iran. Here, for the case study, the planning horizon includes four time periods and each one takes one season (three months). Feedstock supply: As mentioned by Searcy & Hess (2010), by using an advanced uniform-feedstock preprocessing, it is possible to use single conversion technology for ethanol production from multiple types of feedstock. Two types of biomass resources in Iran, including agricultural residues and forest biomass, are considered. All types of agricultural residues such as corn Stover, wheat straw, and rice straw are taken into account, and forest biomass consists of logging residues and short rotation woody crops including poplar and eucalyptus. The supply amount assessment of biomass resources has been done by Iran Ministry of Energy, Renewable Energy and Energy Efficiency Organization (REEEO), and we have used the provided data by this organization. Generally, in according to the provided information by REEEO, and Maghanaki et al. (2013) and Najafi et al. (2009), the total supply amount of agricultural residues and forest biomass in Iran are approximated as 25.3 and 4.8 million tons (MT) per year, respectively. Further, forest biomass can be harvested all the year except the winter season, and the distribution of total yearly amount of agricultural residues is approximated to be 50% in the fall, 30% in the summer, 10% in the spring, and 10% in the winter season. Finally, the geographical distribution of them is illustrated in Fig. 6. Biofuel demand: Based on the available data from Iran Ministry of Energy, the total fuel consumption, except liquefied petroleum gas, is 22.7 BGY (Billion Gallons in a Year) in 2014. In accordance with available feedstock and conservation rate 72.6 gallons per dry ton biomass (Marufuzzaman & Ekşioğlu, 2016), the total biofuel production in Iran, by consideration of the moisture content of the feedstock, can be set to 1.2 BGY that is about 5.2% of the total reported fuel consumption in 2014. Iran is subdivided into thirty one provinces that are shown in Fig. 7, and are considered as demand centers. The provinces’ demand is set based on their population amount. As a consequence, the biofuel demand at each Iranian province as a demand center, in contrary to the feedstock supply, is relatively stable, and is modeled as a deterministic parameter.

547

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 6. Geographical distributions of the feedstock availability.

1

3

Caspian Sea

2

12 25

14

24

18

27

5

20

8

30

11

15

19

22 6

26

28 10 4

9

3 31

13 23 17

21 7

Potential locations for plants Potential locations for storage facilities

16

29

Fig. 7. Potential locations for plants and storage facilities. Table 3 The parameters related to different types of feedstock.

Agricultural residues Forest biomass

No. of suppliers

Moisture content (%)

Average deterioration rate in each period (%)

Average procurement cost ($/ton)

Average storage cost ($/ton)

31 30

15 50

10 12

35 30

8 2

548

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Table 4 The parameters related to production plants.

Small size Medium size Large size

Production capacity (in one year)

Biofuel storage capacity

Cost of production capacity

Cost of storage capacity

Total cost

50 MGY 120 MGY 200 MGY

50 thousand barrels* 100 thousand barrels 200 thousand barrels

$21.85 M $41.95 M $54.53 M

$765000 $1.26 M $2.01 M

$22.61 M $43.21 M $56.54 M

* (1 barrel = 42 Gallons).

Feedstock parameters: Conversion rate, deterioration factor, average procurement cost, moisture content, and storage cost for agricultural residues and forest biomass are reported in Table 3 (Huang et al., 2014; Marufuzzaman & Ekşioğlu, 2016; Parker et al., 2008). Investment costs: Three capacity levels (small, medium, and large), including production and biofuel storage capacity, are considered for the establishment of plants. In accordance with Parker et al. (2008) and Huang et al. (2014), the annualized fixed capital cost is approximated for the plant establishment with different capacity levels (reported in Table 4). It should be mentioned that these costs were estimated by these studies based on a project life of 20 years and an interest rate of 10%. The geographical area of our case study is partitioned based on 31 provinces of Iran. Based on opinions of experts in Iran Ministry of Energy, we have considered that the provinces related to the potential locations for the plant establishment should satisfy two criteria: 1- they have required infrastructure for establishing a biofuel production plant, 2- the biomass supply share of them be more than 5% of the total biomass supply of Iran. As a consequence, 12 potential locations based on Iranian provinces are selected for the plant establishment according to the provided information by Iran REEEO about the mentioned criteria. Fig. 7 illustrates these potential locations. Potential locations for storage facilities are considered in all districts in our case study’s geographical area, Iranian Provinces, and hence 31 potential locations are considered for establishing storage facilities. In Fig. 7, Iranian provinces are numbered and potential locations for storage facilities and plants are shown. Although the actual fixed cost for establishing facilities would vary by location, an average annualized fixed cost for a storage facility with 400,000 tons holding capacity and 1.2 M tons handling capacity in each season, based on a project life of 20 years and an interest rate of 10%, is set to $9.35 M according to the opinions of experts in Iran Ministry of Energy. Further, the activation cost and gain of discontinued use of feedstock storages are set to $174,000/season and $129,000/season, respectively. Transportation costs: In this real-life case study, trucks are available for transportation between SC entities. However, it is possible to use railway transportation between some feedstock storage facilities and plants based on Iran’s railway map. Fig. 8 shows a map of Iran's railway in accordance with Iranian provinces' number in Fig. 7. In this study, a fixed cost for using railway transportation with a predetermined capacity during one time period is assumed for each storage facility. The data related to the capacity of railway routes between storage facilities and plants, the routes’ length for rail transport, the fixed cost of using rail transportation in each storage facility, and the cost of transporting one-unit feedstock (ton) over a distance of one kilometer using railway transport are in according to the information from Iran’s Railway Transportation Company. We have assumed that the capacity of trucks for feedstock shipment is 25 tons, and for biofuel transportation is 8000 gallons. The cost of transporting one-unit feedstock (ton) and biofuel (gallon) between each two SC entities is calculated by multiplying the distance between the entities and the transportation rate per unit of distance (Kilometer). The length of roads between SC entities is obtained from Iran's Road Maintenance and Transportation Organization.

Fig. 8. Iran railway’s map.

549

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Table 5 The parameters related to emission amounts from different processes and their cost. Processes

Unit

CO2

NOx

SOx

VOC

Transportation by using trucks (GHG Protocol, 2015) Transportation by using rail transports (GHG Protocol, 2015) Biofuel production (Jones, 2010) Unit environmental cost (Rabl et al., 2014)

(Kg/ton-kilometer) (Kg/ton-kilometer) (Kg/gallon) ($/kg)

0.2002 0.0285 6.61 0.023

0.031 0.005 0.0117 5.7

0.00011 – 0.0019 6.6

0.0051 0.0008 0.0189 4.4

Disruption Probabilities of storage facilities: Earthquake, flood, and fire are main disaster events in Iran’s cities with considerable likelihoods (GAR, 2009). The probability of disruption occurrence in potential locations for feedstock storage facilities in each time period is approximated based on the Global Assessment Report on disaster risk reduction (GAR, 2009) and the data provided by Natural Disaster Research Institute in Iran. In addition, in accordance with Fattahi et al., (2017b), the capacity lost percentage of storage facilities due to any disruption event is assumed that follows a normal distribution with mean 60 and standard deviation 8. Emissions data: in accordance with Jones (2010), four main emission types including CO2, NOx, SOx and VOC are taken into account. The emission data related to the transportation and biofuel production and the unit environmental cost for these emission types are obtained from the existing research studies and scientific reports that are shown in Table 5. Social score of facilities: The social score of facilities are obtained by using S-LCA method. As mentioned in Section 3.4, three inventory indicators including social acceptance, creating job opportunities, and annual turnover are considered for scoring potential biofuel plants with different capacities. Fuzzy AHP approach, introduced by Chang (1992) and Chang (1996), is used to obtain the values of parameter ρc, j, e . The details of S-LCA method, fuzzy AHP approach, and the final social scores of potential plants for different inventory indicators are presented in Appendices A–C. 9. Computational study In this section, the results from the computational study are summarized and interpreted. The goal of this section is to evaluate the performance of the proposed model, BD algorithm, rolling horizon approach, and scenario tree construction method, and derive important managerial insights about the biofuel SC planning. The developed MSSP is solved using GAMS 24.1 by CPLEX solver. BD algorithm is coded with MATLAB R2014a and CPLEX12.0. It should be mentioned that a personal computer with Intel Core i7-640 M CPU (2.8 GHz), with 4.00 GB of RAM, is used for all implementations in this section. 9.1. Assessing the performance of the MILP model and BD In this section, the performance of the presented model and BD algorithm are investigated by using several generated test problems. The input parameters of the test problems are randomly generated based on the presented data for the case study. In Table 6 Results from solving the test problems with CPLEX and BD algorithm. Test problems

Test problem 1 (4, 10, 42) Test Problem 2 (5, 12, 41) Test Problem 3 (8, 18, 32) Test problem 4 (10, 20, 31) Test problem 5 (12, 22, 41) Test problem 6 (15, 24, 39) Test problem 7 (16, 25, 38) Test problem 8 (18, 28, 42) Test problem 9 (20, 30, 37) Test problem 10 (25, 35, 41)

CPLEX Solver

BD algorithm

Optimal solution

CPU time (S)

Upper bound

3.39E+08

316

3.41E+08

4.21E+08

602

4.47E+08

Optimality gap between CPLEX and BD Lower bound

CPU Time (S)

BD optimality gap

3.33E+08

452

2.40%

0.58%

4.23e+08

4.12E+08

923

2.66%

0.48%

806

4.51E+08

4.39E+08

2677

2.73%

0.89%

5.29E+08

957

5.30E+08

5.11E+08

3821

3.58%

0.19%

5.71E+08

2431

5.73E+08

5.58E+08

2672

2.68%

0.36%

6.50E+08

3604

6.56E+08

6.40E+08

3987

2.50%

0.92%

5.98E+08

4524

6.05E+08

5.89E+08

4950

2.71%

1.17%

5.93E+08

6112

5.99E+08

5.83E+08

5314

2.74%

1.01%

8.39E+08

5376

8.45E+08

8.21E+08

6202

2.93%

0.72%

1.08E+09

7542

1.10E+09

1.07E+08

6169

2.99%

1.85%

550

(

(UB − LB ) LB

× 100% )

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Scenario fan

Stage 1

Stage 2

Stage 3

Stage 4

Fig. 9. The scenario tree construction procedure.

Table 6, the considered characteristics for the test problems are reported as (N, M, S) in which N is the number of suppliers, potential storage facilities, and customers, M is the potential number of plants, and S is the number of scenarios. Further, the experimental results from solving test problems using CPLEX solver and BD are reported. In Table 6, the optimality gap between CPLEX and BD algorithm is calculated based on the relative difference between the best BD upper bound and the obtained optimal solution from CPLEX solver. It should be mentioned that for the scenario tree construction in the test problems, a scenario fan including 150 scenarios is generated and then, the scenario fan is reduced and converted to a scenario tree by considering εrel = 0.55. Further, in Table 6, for the implementation of BD algorithm, reaching 100 iterations or a Benders gap less than 3% is considered as the termination criterion. The results in Table 6 show the presented mathematical model is solvable by CPLEX solver in a reasonable CPU time. In terms of CPU time, we can observe that in a few test problems, including test problem 8 and 10, BD has a better performance in comparison with CPLEX to reach a Benders gap less than 3%. Moreover, it is highlighted that the number of scenarios has a main influence on the computational tractability of test problems, and hence it is important to generate an efficient set of scenarios for the stochastic problem. In solving the test problems, the average and maximum optimality gap between the CPLEX and BD are 0.82% and 1.85%, respectively that shows the acceptable performance of BD algorithm.

9.2. Assessing the performance of the stochastic solution To construct a scenario tree for the case study, a scenario fan including 100 scenarios is generated for the stochastic parameters, and then, the scenario tree construction approach is implemented in which εrel is set to 0.6. Finally, a scenario tree with 28 scenarios is obtained. The process of scenario tree construction is shown in Fig. 9.

Fig. 10. Simulation of the stochastic and deterministic solution. 551

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Table 7 Results from simulating the test problems and case study.

Case study Test problem Test problem Test problem Test problem

1 2 3 4

Optimal objective value (Z ∗ )

Expected of simulation response (ES)

|Z ∗ − ES| ES

1.29E+09 3.39E+08 4.21E+08 4.47E+08 5.29E+08

1.33E+09 3.44E+08 4.25E+08 4.60E+08 5.40E+08

3.01% 1.42% 0.94% 2.91% 2.03%

× 100%

The stochastic model is solved for the case study, and the optimal expected objective value is equal to 1.29E+09. It is worth mentioning that in the model, the unit penalty cost for unsatisfied demands is set at 4 $/gallon, and constraints (17) and (18), constraints related to the environmental cost and the social score of the network, are not considered. Then, the model is solved in which the fall season is considered as the first period. In this sub-section, the stochastic model’s solution is evaluated using the rolling horizon procedure. For this purpose, 100 realizations of stochastic parameters are considered from the first period until the last one, and the average objective function from the simulation is obtained. In Fig. 10, the cumulative strategic investment cost and the average of tactical cost is presented in each period. Moreover, to compare the stochastic and deterministic solution, the solution from solving the deterministic model in which the stochastic parameters are set to their expected value is also simulated by the same realizations of the stochastic parameters. In Fig. 10, the expected response from simulating the deterministic solution is also presented. As shown by Fig. 10, the expected of simulation responses for the stochastic and deterministic solution are equal to 1.33E+09 and 1.40E+09, respectively. Therefore, in our case study, the relative saving percentage by using the stochastic programming model in comparison with the deterministic one is 5.11%. Although the superiority of the stochastic model in compared with the deterministic model is meaningful, the value of the stochastic solution can be more for other case studies with higher deviations of stochastic parameters. In the objective function of the stochastic program, the investment cost related to the establishment of biofuel production plants and storage facilities is deterministic and the tactical cost related to the operation of the SC is scenario dependent. The expected of the tactical cost for the solution of the model and simulation response are 6.81E+08 and 7.14E+08, respectively. In Table 7, in addition to the case study, the results from simulating test problems 1, 2, 3, and 4 are also reported. In the rolling horizon procedure, the stochastic model should be solved several times for each scenario path, and hence we have selected these test problems because of their small run time using CPLEX solver. The presented rolling horizon framework can be used to make the optimal policy from the stochastic program implementable, and this approach is applied to simulate the SC system’s cost under uncertainty and over a multi-period horizon. As shown in Table 7, the insignificant difference between the expected response from the rolling horizon simulation and optimal solution confirms that the proposed MSSP is applicable as an approximation of the reality.

9.3. Assessing the performance of the scenario construction approach A scenario tree construction approach has to be evaluated in terms of its quality. In this regard, two main stability criteria are introduced in the context of the stochastic programming including: in-sample and out-of-sample stability (Kaut & Wallace, 2007). In this study, the MC simulation is used to generate discrete scenarios related to stochastic parameters of the problem and a forward scenario construction approach is exploited for constructing a scenario tree. In the scenario construction approach, by increasing the value of εrel , the number of obtained scenarios decreases relatively. Although a small number of scenarios results in a stochastic problem with better computational tractability, it may not correctly capture the underlying uncertainty in the problem. Further, when we generate scenarios several times with the same input Table 8 In-sample stability analysis of the scenario tree construction for the case study. Number of scenarios

12 14 11 18 16 22 28 31 26

εrel

0.7 0.7 0.7 0.65 0.65 0.65 0.6 0.6 0.6

CPLEX results

in-sample stability error

Objective value

CPU time (S)

1.24E+09 1.21E+09 1.27E+09 1.29E+09 1.28E+09 1.26E+09 1.29E+09 1.30E+09 1.27E+09

729 783 698 838 842 933 1282 1421 1198

552

4.84%

2.35%

2.33%

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Value of upper and lower bound

6

x 10

9

5 4 3

Upperbound

2 1 0

Lowerbound 0

10

20

30

40

50

60

70

80

90

BD iterations

Fig. 11. BD convergence for the case study.

parameters, different sets of discrete scenarios will be generated and hence different scenario trees will be obtained. In-sample stability guarantees that whichever of these scenario trees is used in the MSSP, the optimal value of the objective function is approximately the same. In this sub-section, different scenario fans including 100 scenarios are generated by MC simulation, and then reduced and converted into a scenario tree by the presented scenario tree construction approach and assuming various values for εrel. Results from solving the MSSP related to the case study are reported in Table 8. In Table 8, the in-sample stability error for each value of εrel is calculated as (Max of objective values - Min of objective values) (Average of objective values) × 100%. As reported in Table 8, in the case study, the difference between the objectives of the stochastic problem with different scenario trees is insignificant when εrel is equal to 0.6, and the in-sample stability error is 2.33%. Further, it is highlighted that the stochastic model is solvable by using the commercial solver CPLEX for the case study. On the other hand, out-of-sample stability guarantees that true objective value of the problem that can be obtained from the simulation of the optimal policy is also approximately the same as the objective value of the stochastic model. In this problem, the expected of rolling horizon simulation’s response that is equal to 1.33E+09 is considered as the true objective value of the problem. Therefore, the out-of-sample stability error, the relative difference between the simulation’s response and the optimal objective value, is 3.1%. Therefore, in the case study, we can conclude that the scenario tree construction approach has reasonable in-sample and outof-sample stability for the MSSP, and the considered scenario tree, including 28 scenarios, is suitable.

9.4. Results from solving the stochastic model for case study In this sub-section, the results from solving the stochastic model for the case study by CPLEX solver and BD algorithm are reported. Fig. 11 shows the convergence of BD algorithm for the case study, and in Table 9, the obtained results are reported. It should be mentioned that in Table 9, for the implementation of BD algorithm, a Benders gap less than 1% is considered as the termination criterion. The optimal solution related to the location of plants and storage facilities is illustrated in Fig. 12. As shown by Fig. 12, 9 biofuel production plants with highest capacity level are selected out of 12 potential locations and 11 feedstock storage facilities are established from 31 candidate sites. We can see that the most of storages and plants are located near to the feedstock fields, and the plants and storages are established in the neighborhood of each other to reduce the transportation cost. The optimal expected objective value is equal to 1.29E+09, Strategic cost of the system related to the establishment of plants and storage facilities equals to 6.12E+08, about 47% of the total system cost. Here, the storage facilities’ cost is equal to the sum of their investment cost, activation cost, and recovery cost due to disruption events, and is 1.17E+08. The tactical cost caused by the operation of the biofuel SC system over a one-year planning horizon takes more than 50% of the total system cost, and are scenariodependent. The breakdown of the total system tactical cost is illustrated in Fig. 13. Among these costs, production cost of the biofuel is recognized as the main cost, about 50% of the total tactical cost. Procurement, transportation, and inventory cost take 20%, 25%, and 6% of the total tactical cost, respectively. These results highlight the importance of logistics planning in the system and the major

Table 9 Results from solving the case study by CPLEX solver and BD algorithm. CPLEX solver

BD Algorithm

Optimal objective function

CPU time

Upper bound value

Lower bound value

Benders optimality gap (%)

CPU time

# of iterations

1.29E+09

842

1.2985E+09

1.2867E+09

0.91%

2395

88

553

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Caspian Sea

Optimal location of plants Optimal location of storage facilities

Fig. 12. Optimal location of plants and storage facilities.

Stock-out cost 0%

Inventory costs 6% Transportation cost 25%

Production cost 49%

Procurement cost 20% Fig. 13. Breakdown of the system’s tactical cost in the optimal objective value.

role of stocking feedstock and biofuel through the SC system. By holding inventory for biofuel and feedstock, the system will be enabled to adjust the balance between the supply and demand over time periods (seasons). In Fig. 14, the expected values of SC’s decisions are illustrated. This figure illustrates how feedstock and biofuel should be stocked in storage facilities and plants, respectively and how feedstock and fuel transport through the network. In Fig. 14, in the network, the importance of direct biomass flows from suppliers to production plants is highlighted. Moreover, it is illustrated that stocking feedstock in storages and biofuel in production plants make balance between supply and demand over the planning horizon due to the seasonality of biomass resources. As illustrated by Fig. 14, the storage facilities are not used in the summer season, and the average of active storages in fall, winter, spring, and summer equals to 11, 11, 9.9, and 0, respectively. The expected of available and used different types of feedstock are illustrated by Fig. 15. Fig. 15 presents the amount of available feedstock types in different seasons for the SC. In the case study, agricultural residue has higher availability in comparison with forest biomass. However, the choice of consuming one type of feedstock is dependent on its availability and characteristics including the moisture content, deterioration rate, procurement cost, and geographical distribution. Further, it is shown by Fig. 15 that the SC in the case study would be able to response to the demand growth relatively because of unused available resources.

554

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan Biomass resources

Storage facilities

Biofuel production plants

Demand centers

A: Agricultural residues

Fall F: Forest biomass

A: Agricultural residues

Winter F: Forest biomas

A: Agricultural residues

Spring F: Forest biomass

A: Agricultural residues

Summer F: Forest biomass

Biomass transportation by rail transports

Storage amount of agricultural residues at the beginning of the period

Storage amount of forest biomass at the beginning of the period

Storage amount of biofuel at the beginning of the period

Fig. 14. Expected value of flows and inventory amounts through the SC network.

Fig. 15. Expected value of available and used biomass resources over the planning horizon.

9.5. Sensitivity analyses In this section, we discuss about the impact of some main parameters’ value on the results of the stochastic model. 9.5.1. Analyzing the impact of the unit penalty cost Based on the optimal costs of the SC network, by setting the unit penalty cost for unsatisfied biofuel at demand centers to 4 $/gallon, the stock-out amount would be obtained at its possible minimum level. In this subsection, to investigate the importance of this parameter value, we consider different unit penalty costs ranging from $0.4 to $4, then, Fig. 16 illustrates the average nonsatisfied percentage of demand for different values of the unit penalty cost. The unit penalty cost can be interpreted as the unit cost of a substitute fuel for biofuel. Therefore, as shown by Fig. 16, the 555

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 16. Expected value of the stock-out amount in different values for unit penalty cost.

Fig. 17. The sensitivity of optimal objective value to the unit penalty cost.

Fig. 18. Impact of demand level on the optimal objective value and non-satisfied demand.

economic potential of biofuel production is dependent on the unit cost of available substitute fuels. Further, in Fig. 17, the sensitivity of the optimal objective value to the unit penalty cost is illustrated for the case study and several test problems. As shown by Fig. 17, the optimal objective value is sensitive to the unit penalty cost for non-satisfied demand meaningfully. Further, there is a limit for the penalty cost in biofuel SCs that by decreasing the penalty cost less than this limit; it is not economically attractive to satisfy the demand of all customer centers.

9.5.2. Analyzing the impact of the demand level In this sub-section, we investigate the impact of the demand level on the optimal objective value and the expected of stock-out amount. By increasing the demand level from 0% to 40%, Fig. 18 shows the percentage increase in the optimal objective value and the expected non-satisfied percentage of the demand. As shown by Fig. 18, it is possible to design the network in order to cover 20% increase in the demand level without stock-out. However, for more than 20% increase, the stock-out amount increases because of limited biomass resources. The sensitivity of the economic objective function to the demand level can be different in various case studies based on the available biomass resources and the penalty cost for non-satisfied demand. In Fig. 19, the sensitivity of the optimal objective function to the changes of the demand level is illustrated for three test problems. 556

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 19. Impact of demand level on the optimal objective value in test problems.

Fig. 20. Impact of disruptions’ characteristics on the investment cost and the expected of tactical cost.

9.5.3. Analyzing the impact of disruptions on the SC costs In this study, it is assumed that if a disruption occurs at a storage facility, the percentage of handling and holding capacity of the storage which would be disrupted over one time period follows a normal distribution. The values for the moments of the normal distribution can be approximated by the recovery time of the storage and the disruption magnitude. In fact, for longer recovery time and higher magnitude of disruptions, the mean value of the normal distribution increases. Fig. 20 illustrates the impact of the normal distribution’s mean on the SC costs. Fig. 20 illustrates that by increasing the mean of the normal distribution related to the disrupted capacity of storage facilities under disruption events, the optimal objective of the stochastic problem including the investment cost and the expected of tactical

Fig. 21. the impact of disruption probability and magnitude on the optimal objective value. 557

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Table 10 Payoff table for the case study. Optimization problem

Z1

Z2

Z3

Min Z1 Min Z2 Max Z3

1.29E+09 2.60E+09 2.60E+09

9.98E+09 6.86E+09 1.17E+10

859.9 905.15 1255.6

cost relatively increases. However, as illustrated by Fig. 20, in some situations, by increasing the mean value, the investment cost remains unchanged and for these situations, the increase of the expected of tactical cost is higher. Therefore, for SC planning under disruption events with higher impact, increasing the investment on the SC structure is dependent to the problem’s state and parameters. Further, three different conditions are considered for several test problems as follows: Condition (1): if a disruption event occurs at a storage facility, whole capacity of the facility would be disrupted, Condition (2): In the scenario tree construction, the probability of disruption events at various storage locations is multiplied by 2.5, Condition (3): both of conditions (1) and (2) are considered, In Fig. 21, the percentage increase of the optimal objective value for each condition is illustrated. The impact of disruption events on the optimal objective value is dependent on the disruption modeling in the stochastic problem. Fig. 21 shows how the objective value is sensitive to the increase of the magnitude and probability of disruption events. 9.6. The ε-constraint method In this study, in addition to the economic objective, social scores of the plant establishment and environmental impacts of the SC are considered. To deal with these issues as other problem’s objectives, ε-constraint method is applied. By using this method, we optimize the economic objective function by considering the other objectives as constraints. Therefore, the model is written as follows: Economic ⎧ Min Z1 = Z ⎪ S. t: ⎪ Z2 ⩽ ε Env , ⎨ Soc ⎪ Z3 ⩾ ε ⎪Constraints (2−30) ⎩

where UBEnvCost and LB SocScore are presented by Z2 and Z3 as decision variables, respectively. In order to properly apply the ε-constraint method, the ranges of Z2 and Z3 are needed. The most common approach is to calculate these ranges from the payoff table. In this study, firstly the optimal value of Z1 is obtained individually without consideration of constraints related to other objective functions, and is presented as Z1∗. Because of the importance of the economic objective, we have prevented that the value of Z1 becomes more than Z1U in obtaining the ranges of Z2 and Z3, and here, we have set Z1U to 2 × Z1∗. As a consequence, in the payoff table, the optimal ⎧ Min Z2 ⎧ Max Z3 , and Z3∗ = S. t: Z1 ⩽ Z1U , respectively. Further, the values values of Z2 and Z3 are calculated by solving Z2∗ = S. t: Z1 ⩽ Z1U ⎨ ⎨ − Constraints (2 30) Constraints (2 30) − ⎩ ⎩ of other objectives in according to the optimal solution of each optimization problem are obtained. For our case study, Table 10 shows the obtained payoff table.

Fig. 22. Impact of minimum acceptable rate for the social score on the objective function. 558

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Fig. 23. Impact of the allowable bound for the expected environmental cost on the objective function.

Optimal objective function

Firstly, we analyze how the minimum acceptable rate for the social score in constraint (18) influences on the optimal objective value. It should be noted that in the optimization model, the scores of social criteria for biofuel plants that are obtained from fuzzy AHP are multiplied by 1000, and the importance weight of each social criterion is considered to be one. The lower and upper values for the minimum acceptable rate of the social score are obtained based on the payoff table, and equal to 859.9 and 1255.6, respectively. Fig. 22 illustrates the sensitivity of the optimal economic objective value to the minimum acceptable rate of the social score in which seven discrete values from 859.9 to 1255.6 are assumed for LB SocScore . The obtained results presented by Fig. 22, show that to increase the social score of biofuel production plants by 5%, the optimal objective value increases about 1.2%. Moreover, the maximum acceptable rate for the social score in constraint (18) so that the stochastic model leads to a feasible solution is 1255.6. Therefore, by the presented setting for the case study, it is possible to increase the social score of the network about 45%. In the optimal solution of the stochastic model without consideration of the environmental cost bound, the average environmental cost of the biofuel SC network in each season has been obtained around $ 9.98E+06 M. Although this cost is high, the limit for this cost is based on the emissions’ trading that is government-mandated. Here, we analyze how the upper-bound limit for the expected environmental cost in each season influences on the expected objective value of the stochastic problem. To analyze the sensitivity of the optimal economic objective to the maximum allowable rate for the environmental cost, based on the payoff table, the lower and upper value for UBEnvCost equal to 6.86E+09 and 9.98E+09, respectively. Fig. 23 shows the sensitivity of the objective function value to the maximum allowable rate for the expected environmental cost in each season in which seven discrete values from 6.86E+09 to 9.98E+09 are taken into account for UBEnvCost . As shown by Fig. 23, if the allowable bound for the expected environmental cost in each season decreases by 5%, the optimal objective value of the stochastic problem increases about 3%. Finally, we have analyzed the combined impact of social responsibility and environmental cost limits on the economic objective to have a sustainable SC. Here, five values are considered for LB SocScore , and for each of them, Fig. 24 shows how the economic objective changes at different levels for UBEnvCost . As it is expected, the higher values for LB SocScore and lower values for UBEnvCost lead to harder constraints in the optimization problem and increase the optimal value of the economic objective.

1.55E+09

Social LB=850 Social LB=1050 Social LB=1250

Social LB=950 Social LB=1150

1.50E+09 1.45E+09 1.40E+09 1.35E+09 1.30E+09 1.25E+09 6.50E+09

7.50E+09

8.50E+09

9.50E+09

maximum allowable rate for the expected environmental cost in each season Fig. 24. Combined impact of LB SocScore and UBEnvCost on the economic objective.

559

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

10. Managerial insights The proposed MSSP is applicable in evaluating the economic potential and infrastructures needs to establish a sustainable biofuel SC system under possible disruptions and feedstock uncertainty and seasonality. The model can help managerial board of biofuel companies to integrate the related tactical and strategic planning decisions over a multi-period setting under uncertainty, and thus by employing the stochastic model which can be easily adapted in accordance with their biofuel SC structure, they would improve their SC’s efficiency. The economic potential of planning a biofuel SC network is strongly dependent on the SC’s demand level and the price of substitute products. By decreasing the substitute products’ price that is considered as the penalty cost for unsatisfied demands, the biofuel companies may wish to not satisfy a part of SC’s customers. In such a situation, servicing this part of customers yields additional costs that are higher than the corresponding price of substitute products. Three important issues are highlighted by this study to mitigate the potential risk related to supply seasonality of feedstock: 1 – consideration of a biofuel SC as an integrated and entire system with vast geographically dispersed biomass resources, 2 – use of different types of feedstock in biofuel production plants with advanced technologies, and 3 – storage function of feedstock and biofuel, although taking an insignificant share of the total SC cost. To deal with the uncertainty related to disruption events, the MSSP is flexible to adjust tactical decisions and decisions related to activation/discontinued use of storage facilities dynamically under disruption events. In this study, the impact of disruption events on the capacity of storages is parameterized based on their magnitude and recovery time. Generally, disruption consideration results in increasing cost of SC planning; however, in according to our computational results on the case study, it is not always optimal to increase the investment to cope with disruption events. The presented framework for the sustainable planning of biofuel SCs will help decision makers to make balance between the economic objective, the environmental cost, and the social impact of the SC. For the presented case study, to separately reduce the network’s environmental cost and improve the social impact by 5%, the optimal economic objective value increases about 3% and 1.2%, respectively. 11. Conclusion This study is a new attempt to propose an MSSP for integrated strategic and tactical planning in biofuel SC networks with disruption risk, feedstock seasonality and uncertainty, and sustainability consideration. To consider the environmental performance and social impact, greenhouse gas emissions, economic development, life-threatening issues, and social expectations are taken into account in the presented stochastic model by related costs and constraints to create a framework for sustainable planning in this context. The developed model is meaningfully sensitive to policies related to the SC sustainability. The multi-stage stochastic model has been applied to investigate the economic potential of converting biomass resources, including agricultural residues and forest biomass, to ethanol in Iran. The computational results show that among different types of tactical cost related to the SC operation, production and transportation cost have the highest portions. Moreover, the high impact of the feedstock seasonality and uncertainty on the strategic and tactical planning decisions is highlighted. In the stochastic model, in addition to the inherent uncertainty evoked by uncertain supply of biomass resources, disruption events results in stochastic reduction of facilities’ capacity. To hedge well against the disruption risk, it is important to consider the disruption risk in SC planning models. In this study, to construct an efficient scenario tree for presenting the stochastic parameters in the MSSP, a MC simulation method is used for obtaining a fan of scenarios, and then the number of these scenarios is reduced by applying a forward scenario selection technique successively to construct a suitable scenario tree. Computational results show that this approach has in-sample and out-ofsample stability. For the first time in the context of SC planning, a rolling horizon approach is presented to simulate the optimal policy of the model and make planning decisions of the MSSP implementable. In this approach, a rolling horizon model is presented to obtain decisions in each time period over the planning horizon. The results from the rolling horizon approach demonstrate that the MSSP is suitable for biofuel SC planning under uncertainty. It is also shown that for the real life case study, the percentage cost saving from simulating the stochastic model’s policy in comparison with deterministic one is 5.11%. There are many opportunities to extend the presented model: 1. the integration of other planning decisions and uncertainties such as biofuel demand and cost parameters to the stochastic model can lead to a more comprehensive model for biofuel SC planning, 2. the development of mitigation and contingency strategies in face of SC risks can result in an efficient risk-averse model, 3. depending on the used biomass, considering other criteria for the environmental performance and/or social impacts related to the cultivation stage of biofuel SCs would be interesting. Generally, developing MSSPs and presenting efficient solution approaches for them can be suitable for biofuel SC planning under uncertainty and over a multi-period planning horizon, which can be paid more attentions by the scholars in other real-life case studies. Finally, the robust optimization is a suitable tool for biofuel SC planning in which historical data may be scarce and cannot be described by probability distributions accurately, and will be another promising research direction. Appendix A. Implementation of the S-LCA methodology As presented by UNEP/SETAC (2009), the implementation of the S-LCA methodology includes four main steps that are explained in this section to obtain social scores of different location and capacity levels for establishing biofuel production plants. 560

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Table A1 Five stakeholder categories and 31 impact subcategories (UNEP/SETAC, 2009). Stakeholder categories

Impact subcategories

Worker

Child labor Fair salary Working hour Forced labor Equal opportunities Freedom of association and collective bargaining Health and safety Feedback mechanism Consumer privacy Transparency Health and Safety End of life responsibility Access to material resources Access to immaterial resources Delocalization and Migration Cultural Heritage Safe & healthy living conditions Respect of indigenous rights Community engagement Local employment Secure living conditions Public commitments to sustainability issues Contribution to economic development Prevention & mitigation of armed conflicts Technology development Corruption Fair competition Promoting social responsibility Supplier relationships Respect of intellectual property rights

Consumer

Local Community

Society

Value chain actors (not including consumers)

The first step is to define the objective and scope of the S-LCA methodology. In this study, our objective is to obtain the social scores of potential location and capacity level of biofuel production plants. Based on our goal, the scope is the entire stages of the biofuel production system from biomass suppliers to biofuel markets. In the second step, to reach our objective, several inventory indicators (social indicators) should be determined. Inventory indicators provide the most direct evidence of the condition or result they are measuring, and can be qualitative or quantitative. As suggested by UNEP/SETAC (2009), all five stakeholder categories and 31 impact subcategories (See Table A1) are presented to the panel of experts. In our study, the panel includes managers and engineers in the area of biofuel/energy companies and other stakeholders such as local communities. They are interviewed and the goal and scope are described. Based on their opinions three main impressionable stakeholders and six impact subcategories that have more influence on them are taken into account as a generic list for selecting location and capacity level of plants (See Fig. 3). Finally, in order to measure the social scores of different potential locations and capacity levels for plants, three inventory indicators, including job creation, annual turnover, and social acceptance are developed. In the third step, to calculate the social scores of different potential locations and capacity levels for the plant establishment, we have applied the fuzzy AHP approach. In this approach, the problem hierarchy is constructed based on the inventory indicators that are developed in the second step. By using fuzzy AHP, pairwise comparison matrices for different inventory indicators are obtained based on linguistic preferences of an experts’ panel, and then we have obtained quantitative scores. The details of the fuzzy AHP implementation are explained in Appendix B. Finally, in the fourth step, we can compare the social scores of different potential locations and capacity levels for establishing biofuel production plants. In this study, we have used the obtained scores in the third step as parameter ρc, j, e , the score of social indicator c for establishing biofuel production plant j with capacity level e, to guarantee that the social score of established plants in the network be more than a predetermined target. This target is imposed by a constraint into the optimization problem. Appendix B. Implementation of the fuzzy AHP The fuzzy AHP approach is implemented through four main steps as follows: Step 1: the problem hierarchy is constructed based on the three developed inventory indicators (annual turnover, job creation, and social acceptance) through the S-LCA methodology. Fig. B1 shows the considered problem hierarchy. Step 2: The pairwise comparison of location and capacity level of plants is done for each social indicator by using a balanced multidimensional experts group. The comparisons are done by asking the question which of two potential location j with capacity level e and potential location j’ with capacity level e’ has more importance for establishing a biofuel production plant in terms of a social 561

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Social indicators

location 1, capacity level 1

Annual turnover

Job creation

location 1, capacity level 2

location j, capacity level e

Social acceptance

Location |J |, capacity level |E|

Fig. B1. Hierarchy of our decision making problem.

Table B1 Triangular fuzzy numbers corresponding to the linguistic preferences (source: Tolga et al. (2005)). Linguistic preference

Demonstrated importance (DI)

Very strong importance (VSI)

Strong Importance (SI)

Moderate importance (MI)

Equal importance (EI)

Triangular fuzzy scale Triangular fuzzy reciprocal scale

(2, 5 2, 3) (1 3, 2 5, 1 2)

(3 2, 2, 5 2) (2 5, 1 2, 2 3)

(1, 3 2, 2) (1 2, 2 3, 1)

(1 2, 1, 3 2) (2 3, 1, 2)

(1, 1, 1) (1, 1, 1)

Table B2 The final quantitative social scores of the plant establishment in different potential locations with various capacity levels (×1000). Potential location

Capacity level

Indices

Provinces

1

Tabriz

2

Isfahan

3

Tehran

4

Khorasan razavi

5

Khuzestan

6

Shiraz

7

Ghazvin

8

Kerman

9

Kermanshah

10

Mazandaran

11

Arak

12

Yazd

1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

Social indicators

Sum

Social acceptance

Job creation

Annual turnover

34.58 34.58 34.58 40.23 40.23 40.23 0.00 0.00 0.00 13.91 13.91 13.91 50.40 50.40 50.40 8.86 8.86 8.86 24.10 24.10 24.10 10.15 10.15 10.15 24.73 24.73 24.73 11.09 11.09 11.09 50.47 50.47 50.47 64.81 64.81 64.81

2.55 9.11 18.10 23.33 32.79 45.77 7.25 9.96 13.67 17.23 23.67 32.50 33.15 44.79 60.74 18.95 26.04 35.74 8.70 11.95 16.41 31.15 42.79 58.74 31.15 42.79 58.74 28.03 38.51 52.87 19.57 23.15 28.05 13.70 16.95 21.41

21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97 21.15 22.21 39.97

562

58.28 65.90 92.65 84.71 95.24 125.97 28.40 32.17 53.65 52.29 59.79 86.38 104.70 117.40 151.12 48.96 57.10 84.58 53.95 58.26 80.48 62.45 75.14 108.86 77.02 89.72 123.44 60.27 71.81 103.93 91.20 95.83 118.50 99.66 103.97 126.19

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

indicator. The experts only choose the corresponding linguistic preference to evaluate this question. Then, the linguistic preferences are transformed into the triangular fuzzy numbers that are developed by Tolga et al. (2005). The positive triangular fuzzy numbers related to the linguistic preferences are shown in Table B1. The pairwise comparison matrix for a social indicator results in an N × N matrix, and N equals to the number of potential locations multiplied by the number of capacity levels (|J | × |E|). In our case study, the number of potential locations and capacity levels are 12 and 3, respectively. Therefore, the main pairwise comparison matrix is a 36 × 36 matrix. Step 3: in this step, we check for the consistency of pairwise comparison matrices. In this study, the degree of consistency in comparison matrices are calculated based on Gogus & Boucher (1998) by using MATLAB R2014a software. In this approach, two separate matrices AM and AG corresponding to a fuzzy pairwise comparison matrix are formed based on the mean values of preferences and the geometric means of lower and upper bounds, respectively. The consistency ratios for these matrices are calculated, and if both of them for a comparison matrix are less than 0.1, as suggested by Gogus & Boucher (1998), the matrix is consistent. For our case study, the consistency ratios for each comparison matrix are less than 0.06. Therefore, in our implementation, it is not required to revise the experts’ judgments. For more information about AM, AG, and the calculation of their consistency ratios in pairwise comparison matrices, one can refer to Gogus & Boucher (1998) and Saaty & Vargas (2001). Step 4: in this step, using extent analysis method, proposed by Chang (1992) and Chang (1996), we calculate the score of social indicators for potential locations and capacity levels of biofuel production plants. The final results of our implementation for the case study are illustrated in Table. B2. Moreover, the extent analysis approach that is suggested by Chang (1992) is explained in Appendix C. Appendix C. The extent analysis method on fuzzy AHP Consider X = {x1, x2 , ...,x n} be a set of objects, and U = {u1, u2, ...,um} be a set of goals. Based on the extent analysis of Chang (1992) and Chang (1996), each object is considered and for each goal, the extent analysis is done. As a consequence, m extent analysis value can be obtained for each object. Assume Mg1i , Mg2i , ...,Mgmi , i = 1, 2, ...,n be the signs of m extent analysis for object i. All Mgji (i = 1, 2, ...,n, j = 1, 2, ...,m) are triangular fuzzy number as Mgji = (lij , mij , uij ) . Generally, the value of fuzzy synthetic extent corresponding to ith object can be calculated as follows:

m

Si =

∑ j=1

n

Mgji

m

⎛

⎡ ⎤ ⊗ ⎢∑ ∑ Mgji ⎥ ⎣ i=1 j=1 ⎦

=

m

m

m

⎞

⎜ ∑ lij , ∑ mij , ∑ uij ⎟ j=1 j=1 ⎝ j=1 ⎠

−1 n

m

n

m

n

m

= (li , mi , ui ).

⎛ ⎞ ⎜ ∑ ∑ lij , ∑ ∑ mij , ∑ ∑ uij ⎟ i=1 j=1 i=1 j=1 ⎝i = 1 j = 1 ⎠

In the next step, the degree of possibility of Si = (li , mi , ui ) ⩾ Sk = (lk , mk , uk ) for triangular fuzzy numbers is computed using the following equation:

if mi ⩾ mk , ⎧1 ⎪0 if lk ⩾ ui , V (Si ⩾ Sk ) = ⎨ lk − ui others ⎪ (mi − ui) − (mk − lk ) ⎩ Further, the degree of possibility for a convex fuzzy number s to be greater than k other convex fuzzy number si (i = 1, 2, ...,k ) can be obtained as:

V (s ⩾ s1, s2, ...,sk ) = V [s ⩾ s1, s ⩾ s2, ...,s ⩾ sk ] = min V (s ⩾ si ) , i = 1, 2, ...,k Finally, the weight vector for n objects can be obtained as W ′ = (d′ (A1 ), d′ (A2 ), ...,d′ (An ))t in which d′ (Ai ) = min V (si ⩾ sk ) , k = 1, 2, ...,n and i ≠ k . Via normalization, the normalized final weights can be obtained as W = (d (A1 ), d (A2 ), ...,d (An ))t in which the elements of weight vector are non-fuzzy numbers. In our case study, for each social indicator, the weight vector for the objects, including potential locations and capacity levels for the plant establishment, are separately obtained based on the explained extent analysis method on the fuzzy AHP. The final results are illustrated in Table B2. Appendix D. Scenario reduction heuristics We consider an n-dimensional stochastic process ξ = {ξt }|tT=|1 with distribution function χ and a finite support. Its support can be presented by NS discrete scenarios as supp(ξ ) = {ξ (1) , ξ (2) , ...,ξ (NS ) }, ξ (ω) = {ξt(ω) }|tT=|1, and ω = 1, 2, ...,NS. The scenarios’ probability is NS ∼ ∼ χ is the distribution function of another n-dimensional stochastic process ξ = {ξt }|tT=|1. Let denoted by π (ω) and ∑ π (ω) = 1. Assume ∼ ω=1

∼ ∼(1) ∼(2) ∼(NS′) ∼(ω′) , ω′ = 1, 2, ...,NS′, supp(ξ ) = {ξ , ξ , ...,ξ } and NS′ be its number of scenarios with corresponding probabilities π NS′

∼(ω′) = 1. The optimal solution of linear transportation problem (S. 1) is equal to Kantorovich distance (KD) between χ and ∼ χ. ∑ π

ω′= 1

563

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

KD(χ , ∼ χ ) = min

NS

⎧ ∑ ⎨ ω=1 ⎩

NS′

∑

∼(ω′) ηωω′ c|T| (ξ (ω) , ξ )

NS′

NS

ηωω′ ⩾ 0,

ω′= 1

⎫ ηωω′ = π (ω) , ω = 1, 2, ...,NS, ω′ = 1, 2, ...,NS′ , ⎬ ω′= 1 ⎭ (D1)

∑ ηωω′ = π∼(ω′) , ∑ ω=1

t ∼(ω′) ∼(ω′) where ct (ξ (ω) , ξ ) = ∑ ∥ξ τ(ω) −ξ τ ∥, t = 1, 2, ...,|T| and ∥. ∥ is a norm over Rn . As a consequence, the distance between scenarios τ=1

over the planning horizon equals to c|T | . ∼(ω′) ∼(ω′) χ is the reduced distribution function of χ , its discrete support contains scenarios ξ If we consider that ∼ = {ξ t }|tT=|1, ω′ ∈ {1, 2, ...,NS } DS, and DS is the set of deleted scenarios. In such a situation, for a predetermined set DS ⊂ {1, 2, ...,NS } , the χ ) , can be obtained as follows: χ , KD(χ , ∼ minimum Kantorovich distance between χ and ∼

KD(χ , ∼ χ)=

∑ ω ∈ DS

∼(ω′) π (ω) min c|T| (ξ (ω) , ξ ), ω′∉ DS

∼(ω′) , ω′ ∈ {1, 2, ...,NS } DS , is determined by π ∼(ω′) = π (ω′) + χ, π where the probability of scenarios related to ∼

∑

π (ω) where

ω ∈ DS ω′

∼(ω′) DS ω′ = {ω ∈ DS: ω′ = ω′ (ω)}, and ω′ (ω) ∈ arg min c|T| (ξ (ω) , ξ ), ∀ ω ∈ DS, i.e., a selection of the index set of closest scenarios to ω′∉ DS

ξ (ω) ,∀ ω ∈ DS. In order to optimally select set DS with fixed cardinality #DS , the following reduction problem should be solved: Min

∼(ω′) ⎧ ⎫ c|T| (ξ (ω) , ξ ): DS ⊂ {1, 2, ...,NS }, #DS = NS−N ′ , ∑ π (ω) ωmin ′∉ DS ⎨ ω ∈ DS ⎬ ⎩ ⎭

(D2)

in which N ′ is the number of preserved scenarios after reduction, N ′ = NS−#DS . Dupačová et al. (2003) presented two heuristics for solving reduction problem (S.2) because of its NP-hardness. These heuristics called backward reduction and forward selection that are based on deleting one scenario (#DS = 1) or keeping one scenario (N ′ = 1), successively. In the backward scenario reduction algorithm, until deleting NS−N ′scenarios, optimal deletion of one scenario should be recursively repeated, but in the forward scenario selection algorithm, until reaching N ′ scenarios, optimal selection of one scenario should be recursively repeated. Appendix E. Scenario tree construction approach If we consider scenarios of multivariate stochastic parameters related to distribution function χ in a form of a scenario fan, each scenario ω = 1, 2, ...,NS can be represented as ξ (ω) = (ξ0(ω) , ξ1(ω) , ..., ξ|T(ω| ) ) with occurrence probability π (ω) where |T| is the number of time periods. In the scenarios fan, all scenarios are the same at the first node, i.e. ξ0(1) = ξ0(2) = ...=ξ0(NS ) = ξ0∗, and the total number of nodes is NS × |T| + 1. In the scenario tree construction procedure, the objective is to obtain a scenario tree with probability disχ is less tribution χε based on the scenario fan in which the number of nodes is less and Kantorovich distance (KD) between χ and ∼ than ε (a predetermined value). In this study, to convert the scenario fan into a scenario tree based on the forward scenario tree construction method, proposed by Heitsch & Romisch (2005), Heitsch & Römisch (2009), and Growe-Kuska et al. (2003), from period t = 1 until t = |T|, successive clustering of scenarios is performed. To satisfy KD(χ , χε ) < ε , at each period t, εt is considered for the |T | reduction algorithm under the condition ∑t = 1 εt ⩽ ε . In other words, maximal reduction strategy is performed in period t such that

∑ π (ω) min ct (ξ (ω) , ξ (ω′) ) ⩽ εt .

ω ∈ DS

ω′∉ DS

Moreover,

in

period

t,

the

distance

between

two

scenarios

is

obtained

via

t

ct (ξ (ω) , ξ (ω′) ) = ∑ ∥ξ τ(ω) −ξ τ(ω′) ∥. τ=1

According to Heitsch & Romisch (2005), relation ε = εrel. εmax is typically applied to compute ε . In this relation, 0 < εrel < 1, is a constant value representing a scale for the reduction amount compared with the scenario fan and εmax is the best distance between probability distribution of the initial scenario fan and one of its scenarios with probability one. At period t, for constructing the scenario tree, εt is computed via the following relation (Fattahi et al. 2017b):

εt =

ε ⎡ 1 + q¯ ⎛1− t ⎞⎤ , |T| + 1 ⎢ ⎝ |T| + 1 ⎠ ⎥ ⎣2 ⎦

t = 1, 2, ..,|T|,

(E1)

where q¯ ∈ [0, 1] is a constant parameter, which is considered to be 1. Appendix F. Master and sub-problem related to BD algorithm Benders sub-problem at iteration κ of BD algorithm, SP {κ} , is as follows:

Minimize : EΩ ⎛⎜ ⎝

O ICostP (ω) + O ICostS (ω) + OTCost (ω) + O ProcurementCost (ω) ⎞ ⎟, + O ProductionCost (ω) + O PenaltyCost (ω) ⎠

564

(F.1)

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

∑ Fst, i, m (ω) + ∑ Fst, j, m (ω) ⩽ θmt ,s (ω) i∈I

∀ m ∈ M , ∀ s ∈ S (m), ∀ t ∈ T , ∀ ω ∈ Ω, (F.2)

j∈J

Hit,+m1 (ω) =

Fst, i, m (ω)− ∑

∑ s ∈ S (m)

Hit,+m1 (ω)

βit, m Hit, m (ω)

=

Fit, j, m, l (ω)

∑

(F.3)

Fst, i, m (ω)−

∑

+

s ∈ S (m)

∑

κm Hit, m (ω)

⩽

∀ i ∈ I , ∀ m ∈ M , ∀ t = 1, ∀ ω ∈ Ω,

j ∈ J l ∈ L (i)

cai ψit

(ω) Zit

∑ ∑

Fit, j, m, l (ω)

∀ i ∈ I , ∀ m ∈ M , ∀ t ∈ T {1}, ∀ ω ∈ Ω, (F.4)

j ∈ J l ∈ L (i)

(ω),

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω, (F.5)

m∈M

Fst, i, m (ω) ⩽ chi ψit (ω) Zit (ω).

∑ ∑

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω, (F.6)

m ∈ M s ∈ S (m)

Fit, j, m, l (ω) ⩽ cti, l Wit, l (ω),

∑ ∑

∀ i ∈ I , ∀ l ∈ L (i), ∀ t ∈ T , ∀ ω ∈ Ω, (F.7)

j∈J m∈M

∑

Fit, j, m, l (ω) ⩽ cli, j, l Wit, l (ω),

∀ i ∈ I , ∀ j ∈ J , ∀ l ∈ L (i), ∀ t ∈ T , ∀ ω ∈ Ω. (F.8)

m∈M

∑

P tj, m (ω)

∑ cpj, e Xj, e ,

⩽

m∈M

∀ j ∈ J , ∀ t ∈ T , ∀ ω ∈ Ω, (F.9)

e∈E

V tj (ω) ⩽

∑ csj, e Xj, e ,

∀ j ∈ J , ∀ t ∈ T , ∀ ω ∈ Ω, (F.10)

e∈E

P tj, m (ω) = λm (1−νm ) V tj + 1 (ω)

∑

=

⎛ Ft (ω) + ⎜ ∑ s , j, m ⎝ s ∈ S (m)

P tj, m (ω)−

F tj, k (ω)

∑

m∈M

⎞ Fit, j, m, l (ω) , ⎟ ⎠

∑ ∑ i ∈ I l ∈ L (i)

∀ j ∈ J , ∀ t ∈ T , ∀ m ∈ M , ∀ ω ∈ Ω, (F.11)

∀ j ∈ J , ∀ t = 1, ∀ ω ∈ Ω, (F.12)

k∈K

V tj + 1 (ω) = V tj (ω) +

∑

P tj, m (ω)−

m∈M

∑ F tj, k (ω) + Ukt (ω) = dkt

∑ F tj, k (ω)

∀ j ∈ J , ∀ t ∈ T {1}, ∀ ω ∈ Ω, (F.13)

k∈K

∀ k ∈ K , ∀ t ∈ T , ∀ ω ∈ Ω, (F.14)

j∈J

∑ ∑ ∑ φg, s, i Fst, i, m (ω) + ∑ ∑ ∑ φg, s, j Fst, j, m (ω) ⎡ ⎛ ⎞⎤ m ∈ M s ∈ S (m) i ∈ I m ∈ M s ∈ S (m) j ∈ J ⎢ ⎟ ⎥ ⩽ UBEnvCost , EΩ ⎢ ∑ αg ⎜ + ∑ ∑ ∑ ∑ φl, g, i, j Fit, j, m, l (ω) + ∑ ∑ φg, j, k F tj, k (ω) + ∑ δg P tj (ω) ⎟ ⎥ ⎜ g ∈ G ⎥ ⎢ j∈J k∈K j∈J ⎝ i ∈ I j ∈ J m ∈ M l ∈ L (i) ⎠⎦ ⎣ Xj, e = X (jκ, e) ,

(λ {X} j(,κe) ),

Zit (ω) = Zti,,ω(κ )

(λ {Z} (i,κt), ω)

Wit, l (ω) = Wti,,l(,kω)

∀ j ∈ J , ∀ e ∈ E,

(F.17)

∀ i ∈ I , ∀ l ∈ L (i), ∀ t ∈ T , ∀ ω ∈ Ω.

If we is as follows:

(F.18) (F.19)

X , Z , W , F , H , P, V , U ⩾ 0 (κ ) denote OFSP

(F.15) (F.16)

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω,

(λ {W} (i,κl), t , ω)

∀ t ∈ T.

as the optimal objective value of

SP {κ} ,

then Benders master problem at iteration κ + 1of BD algorithm, MP {κ + 1} ,

Minimize : {O LCostP + EΩ (O DisCost (ω) + OUCostS (ω) + OCCostT (ω))} + φ ,

(F.20)

Zit (ω) ⩽ Yi ,

(F.21)

Wit, l (ω)

⩽

Zit

∀ i ∈ I, ∀ t ∈ T, ∀ ω ∈ Ω , (ω),

∀ i ∈ I , ∀ l ∈ L (i), ∀ t ∈ T , ∀ ω ∈ Ω ,

t Zit − 1 (ω) + Z¯it (ω) = Zit (ω) + Zi ̂ (ω), t Zi ̂ (ω) + Z¯it (ω) ⩽ 1,

⎛

(F.22)

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω,

∀ i ∈ I , ∀ t ∈ T , ∀ ω ∈ Ω,

(F.24)

⎞

∑ ∂c ⎜ ∑ ∑ ρc,j,e Xj, e ⎟ ⩾ LB SocScore

c∈C

⎝ j∈J

(F.23)

e∈E

(F.25)

⎠

565

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

Zit (ω) = Zit (ω′)

∀ i ∈ I , ∀ t = 1 , ∀ (ω, ω′) ∈ Ω,

Wit, l (ω) = Wit, l (ω′), Zit (ω) = Zit (ω′), Wit, l (ω) = Wit, l (ω′), (γ ) OFSP +

∑∑ j∈J e∈E

(F.26)

∀ i ∈ I , ∀ l ∈ L (i), ∀ t = 1 , ∀ (ω, ω′) ∈ Ω,

(F.27) ′

′

′

∀ i ∈ I , ∀ t ∈ T {1}, ∀ (ω, ω′) ∈ Ω: (ξ1ω, ξ2ω, ...,ξtω− 1) = (ξ1ω , ξ2ω , ...,ξtω− 1),

(F.28) ′

′

′

∀ i ∈ I , ∀ l ∈ L (i), ∀ t ∈ T {1} , ∀ (ω, ω′) ∈ Ω: (ξ1ω, ξ2ω, ...,ξtω− 1) = (ξ1ω , ξ2ω , ...,ξtω− 1),

(λ {X} (j,γe) )(Xj, e −X (jγ, e) ) +

∑∑ ∑

(λ {Z} (i,γt), ω)(Zit (ω)−Zti,,ω(γ ) ) +

i∈I t∈T ω∈Ω

∑∑∑ ∑

(λ {W} i(,γl), t , ω)(Wit, l (ω)−Wti,,l(,γω) ) ⩽ φ ,

(F.29) γ = 1, 2, ..., κ

i∈I t∈T l∈L ω∈Ω

(F.30)

X , Y , Z , Z¯ , Z ,̂ W ∈ {0, 1}, φ ⩾ 0.

(F.31)

References An, H., Wilhelm, W.E., Searcy, S.W., 2011. A mathematical model to design a lignocellulosic biofuel supply chain system with a case study based on a region in Central Texas. Bioresour. Technol. 102 (17), 7860–7870. Awudu, I., Zhang, J., 2012. Uncertainties and sustainability concepts in biofuel supply chain management: a review. Renew. Sustain. Energy Rev. 16 (2), 1359–1368. Babazadeh, R., Razmi, J., Pishvaee, M.S., Rabbani, M., 2017. A sustainable second-generation biodiesel supply chain network design problem under risk. Omega 66, 258–277. Benders, J.F., 1962. Partitioning procedures for solving mixed-variables programming problems. Numer. Math. 4 (1), 238–252. Searcy, E., Hess, J.R., 2010. Uniform-format feedstock supply system design for lignocellulosic biomass: A commodity-scale design to produce an infrastructurecompatible biocrude from lignocellulosic biomass. Idaho National Laboratory, Bioenergy Program Report, Idaho Falls. Chand, S., Hsu, V.N., Sethi, S., 2002. Forecast, solution, and rolling horizons in operations management problems: a classified bibliography. Manufact. Service Operat. Manage. 4 (1), 25–43. Chang, D.-Y., 1992. Extent analysis and synthetic decision. Opt. Tech. Appl. 1 (1), 352–355. Chang, D.-Y., 1996. Applications of the extent analysis method on fuzzy AHP. Eur. J. Oper. Res. 95 (3), 649–655. Conejo, A.J., Castillo, E., Minguez, R., Garcia-Bertrand, R., 2006. Decomposition Techniques in Mathematical Programming: Engineering and Science Applications. Springer Science & Business Media. Chen, C.-W., Fan, Y., 2012. Bioethanol supply chain system planning under supply and demand uncertainties. Transport. Res. Part E: Logist. Transport. Rev. 48 (1), 150–164. Dal-Mas, M., Giarola, S., Zamboni, A., Bezzo, F., 2011. Strategic design and investment capacity planning of the ethanol supply chain under price uncertainty. Biomass Bioenergy 35 (5), 2059–2071. Dupačová, J., Gröwe-Kuska, N., Römisch, W., 2003. Scenario reduction in stochastic programming. Math. Program. 95 (3), 493–511. Ekşioğlu, S., Li, S., Zhang, S., Sokhansanj, S., Petrolia, D., 2010. Analyzing impact of intermodal facilities on design and management of biofuel supply chain. Transport. Res. Rec.: J. Transp. Res. Board 2191, 144–151. Ekşioğlu, S.D., Acharya, A., Leightley, L.E., Arora, S., 2009. Analyzing the design and management of biomass-to-biorefinery supply chain. Comput. Ind. Eng. 57 (4), 1342–1352. Elia, J.A., Baliban, R.C., Xiao, X., Floudas, C.A., 2011. Optimal energy supply network determination and life cycle analysis for hybrid coal, biomass, and natural gas to liquid (CBGTL) plants using carbon-based hydrogen production. Comput. Chem. Eng. 35 (8), 1399–1430. Eskandarpour, M., Dejax, P., Miemczyk, J., Péton, O., 2015. Sustainable supply chain network design: an optimization-oriented review. Omega 54, 11–32. Farrell, A.E., Plevin, R.J., Turner, B.T., Jones, A.D., O’hare, M., Kammen, D.M., 2006. Ethanol can contribute to energy and environmental goals. Science 311 (5760), 506–508. Fattahi, M., 2017. Resilient procurement planning for supply chains: a case study for sourcing a critical mineral material. Resour. Policy. https://doi.org/10.1016/j. resourpol.2017.10.010. Fattahi, M., Govindan, K., 2017. Integrated forward/reverse logistics network design under uncertainty with pricing for collection of used products. Ann. Oper. Res. 253 (1), 193–225. Fattahi, M., Govindan, K., Keyvanshokooh, E., 2017a. A multi-stage stochastic program for supply chain network redesign problem with price-dependent uncertain demands. Comput. Oper. Res. https://doi.org/10.1016/j.cor.2017.12.016. Fattahi, M., Govindan, K., Keyvanshokooh, E., 2017b. Responsive and resilient supply chain network design under operational and disruption risks with delivery leadtime sensitive customers. Transport. Res. Part E: Logist. Transport. Rev. 101, 176–200. GAR, 2009. Global Assessment Report on Disaster Risk Reduction: I.R. Iran’s Natural Hazards Profile. Report, available on line at < http://www.grid.unep.ch > . Ghaderi, A., Boland, N., JabalAmeli, M.S., 2012. Exact and Heuristic Approaches to the Budget-Constrained Dynamic Uncapacitated Facility Location-Network Design Problem. The University of Newcastle, Centre for Optimal Planning and Operations. GHG Protocol, 2015. Emission Factors from Cross-Sector Tools. Available on line at. http://www.ghgprotocol.Org/calculationtools/all-tools. Giarola, S., Bezzo, F., Shah, N., 2013. A risk management approach to the economic and environmental strategic design of ethanol supply chains. Biomass Bioenergy 58, 31–51. Gogus, O., Boucher, T.O., 1998. Strong transitivity, rationality and weak monotonicity in fuzzy pairwise comparisons. Fuzzy Sets Syst. 94, 133–144. Gonela, V., Zhang, J., Osmani, A., 2015. Stochastic optimization of sustainable industrial symbiosis based hybrid generation bioethanol supply chains. Comput. Ind. Eng. 87, 40–65. Govindan, K., Fattahi, M., Keyvanshokooh, E., 2017. Supply chain network design under uncertainty: a comprehensive review and future research directions. Eur. J. Oper. Res. 263, 108–141. Growe-Kuska, N., Heitsch, H., Romisch, W., 2003. Scenario reduction and scenario tree construction for power management problems. Paper presented at the Power Tech Conference Proceedings, 2003 IEEE Bologna. Heitsch, H., Romisch, W., 2005. Generation of multivariate scenario trees to model stochasticity in power management. Paper presented at the Power Tech, 2005 IEEE Russia. Heitsch, H., Römisch, W., 2009. Scenario tree reduction for multistage stochastic programs. CMS 6 (2), 117–133. Hombach, L.E., Büsing, C., Walther, G., 2018. Robust and sustainable supply chains under market uncertainties and different risk attitudes – a case study of the German biodiesel market. Eur. J. Oper. Res. 269 (1), 302–312. Hombach, L.E., Walther, G., 2015. Pareto-efficient legal regulation of the (bio) fuel market using a bi-objective optimization model. Eur. J. Oper. Res. 245 (1), 286–295. Hombach, L.E., Cambero, C., Sowlati, T., Walther, G., 2016. Optimal design of supply chains for second generation biofuels incorporating European biofuel

566

Transportation Research Part E 118 (2018) 534–567

M. Fattahi, K. Govindan

regulations. J. Cleaner Prod. 133, 565–575. Huang, Y., Chen, C.-W., Fan, Y., 2010. Multistage optimization of the supply chains of biofuels. Transport. Res. Part E: Logist. Transport. Rev. 46 (6), 820–830. Huang, Y., Fan, Y., Chen, C.-W., 2014. An integrated biofuel supply chain to cope with feedstock seasonality and uncertainty. Transport. Sci. 48 (4), 540–554. Jones, D.L., 2010. Potential air emission impacts of cellulosic ethanol production at seven demonstration refineries in the United States. J. Air Waste Manage. Assoc. 60 (9), 1118–1143. Kaut, M., Wallace, S.W., 2007. Evaluation of scenario-generation methods for stochastic programming. Pacific J. Optim. 3 (2), 257–271. Kannan, D., 2018. Role of multiple stakeholders and the critical success factor theory for the sustainable supplier selection process. Int. J. Prod. Econ. 195, 391–418. Kazemzadeh, N., Hu, G., 2013. Optimization models for biorefinery supply chain network design under uncertainty. J. Renew. Sustain. Energy 5 (5), 053125. Li, Q., Hu, G., 2014. Supply chain design under uncertainty for advanced biofuel production based on bio-oil gasification. Energy 74, 576–584. Lim, M.K., Ouyang, Y., 2016. Biofuel supply chain network design and operations. Environmentally Responsible Supply Chains. Springer. pp. 143–162. Lin, C.C., Wang, T.H., 2011. Build-to-order supply chain network design under supply and demand uncertainties. Transport. Res. Part B: Methodol. 45 (8), 1162–1176. Maghanaki, M.M., Ghobadian, B., Najafi, G., Galogah, R.J., 2013. Potential of biogas production in Iran. Renew. Sustain. Energy Rev. 28, 702–714. Mak, H.-Y., Shen, Z.-J., 2012. Risk diversification and risk pooling in supply chain design. IIE Trans. 44 (8), 603–621. Marufuzzaman, M., Ekşioğlu, S.D., 2016. Designing a reliable and dynamic multimodal transportation network for biofuel supply chains. Transport. Sci. 51 (2), 494–517. Naik, S.N., Goud, V.V., Rout, P.K., Dalai, A.K., 2010. Production of first and second generation biofuels: a comprehensive review. Renew. Sustain. Energy Rev. 14 (2), 578–597. Najafi, G., Ghobadian, B., Tavakoli, T., Yusaf, T., 2009. Potential of bioethanol production from agricultural wastes in Iran. Renew. Sustain. Energy Rev. 13 (6–7), 1418–1427. Osmani, A., Zhang, J., 2013. Stochastic optimization of a multi-feedstock lignocellulosic-based bioethanol supply chain under multiple uncertainties. Energy 59, 157–172. Parker, N., Tittmann, P., Hart, Q., Lay, M., Cunningham, J., Jenkins, B., … Gray, E., 2008. Strategic development of bioenergy in the western states development of supply scenarios linked to policy recommendations Task 3: Spatial analysis and supply curve development. University of California, Davis. Final report, Western Governors’ Association, June 2008: PC. Peng, P., Snyder, L.V., Lim, A., Liu, Z., 2011. Reliable logistics networks design with facility disruptions. Transport. Res. Part B: Methodol. 45 (8), 1190–1211. Poudel, S.R., Marufuzzaman, M., Bian, L., 2016. Designing a reliable bio-fuel supply chain network considering link failure probabilities. Comput. Ind. Eng. 91, 85–99. Poudel, S., Marufuzzaman, M., Quddus, M., Chowdhury, S., Bian, L., Smith, B., 2018. Designing a reliable and congested multi-modal facility location problem for biofuel supply chain network. Energies 11 (7), 1682. https://doi.org/10.3390/en11071682. Poudel, S.R., Quddus, M.A., Marufuzzaman, M., Bian, L., 2017. Managing congestion in a multi-modal transportation network under biomass supply uncertainty. Ann. Oper. Res. 1–43. Quddus, M.A., Chowdhury, S., Marufuzzaman, M., Yu, F., Bian, L., 2017a. A two-stage chance-constrained stochastic programming model for a bio-fuel supply chain network. Int. J. Prod. Econ. 195, 27–44. Quddus, M.A., Hossain, N.U.I., Marufuzzaman, M., Jaradat, R.M., Roni, M.S., 2017b. Sustainable network design for multi-purpose pellet processing depots under biomass supply uncertainty. Comput. Ind. Eng. 110, 462–483. Rabl, A., Spadaro, J.V., Holland, M., 2014. How Much is Clean Air Worth? Calculating the Benefits of Pollution Control. Cambridge University Press. Saaty, T.L., Vargas, L.G., 2001. Models, Methods, Concepts and Applications of the Analytic Hierarchy Process. Springer, New York. Santibañez-Aguilar, J.E., González-Campos, J.B., Ponce-Ortega, J.M., Serna-González, M., El-Halwagi, M.M., 2014. Optimal planning and site selection for distributed multiproduct biorefineries involving economic, environmental and social objectives. J. Cleaner Prod. 65, 270–294. Santoso, T., Ahmed, S., Goetschalckx, M., Shapiro, A., 2005. A stochastic programming approach for supply chain network design under uncertainty. Eur. J. Oper. Res. 167 (1), 96–115. Seddighi, A.H., Ahmadi-Javid, A., 2015. A sustainable risk-averse approach to power generation planning with disruption risk and social responsibility considerations. J. Cleaner Prod. 105, 116–133. Shabani, N., Akhtari, S., Sowlati, T., 2013. Value chain optimization of forest biomass for bioenergy production: a review. Renew. Sustain. Energy Rev. 23, 299–311. Shalaby, E.A., 2013. Biofuel: Sources, Extraction and Determination. In: Liquid, Gaseous and Solid Biofuels-Conversion Techniques. InTech. Sharma, B., Ingalls, R.G., Jones, C.L., Khanchi, A., 2013. Biomass supply chain design and analysis: basis, overview, modeling, challenges, and future. Renew. Sustain. Energy Rev. 24, 608–627. Snyder, L.V., Atan, Z., Peng, P., Rong, Y., Schmitt, A.J., Sinsoysal, B., 2016. OR/MS models for supply chain disruptions: a review. IIE Trans. 48 (2), 89–109. Sodhi, M.S., 2005. Managing demand risk in tactical supply chain planning for a global consumer electronics company. Prod. Oper. Manage. 14 (1), 69–79. Tang, C.S., 2006. Perspectives in supply chain risk management. Int. J. Prod. Econ. 103 (2), 451–488. Tay, D.H., Ng, D.K., Tan, R.R., 2013. Robust optimization approach for synthesis of integrated biorefineries with supply and demand uncertainties. Environ. Prog. Sustain. Energy 32 (2), 384–389. Tolga, E., Demircan, M.L., Kahraman, C., 2005. Operating system selection using fuzzy replacement analysis and analytic hierarchy process. Int. J. Prod. Econ. 97 (1), 89–117. Tong, K., You, F., Rong, G., 2014. Robust design and operations of hydrocarbon biofuel supply chain integrating with existing petroleum refineries considering unit cost objective. Comput. Chem. Eng. 68, 128–139. UNEP/SETAC, 2009. Guidelines for social life cycle assessment of products. United Nations Environment Programme and the Society of Environmental Toxicology and Chemistry, Belgium. Üster, H., Memişoğlu, G., 2017. Biomass logistics network design under price-based supply and yield uncertainty. Transport. Sci. U.S. Environmental Protection Agency, 2007. Renewable Fuels: Regulations and Standards. Xie, F., Huang, Y., Eksioglu, S., 2014. Integrating multimodal transport into cellulosic biofuel supply chain design under feedstock seasonality with a case study based on California. Bioresour. Technol. 152, 15–23. You, F., Tao, L., Graziano, D.J., Snyder, S.W., 2012. Optimal design of sustainable cellulosic biofuel supply chains: multiobjective optimization coupled with life cycle assessment and input–output analysis. AIChE J. 58 (4), 1157–1180. You, F., Wang, B., 2011. Life cycle optimization of biomass-to-liquid supply chains with distributed–centralized processing networks. Ind. Eng. Chem. Res. 50 (17), 10102–10127. Zamboni, A., Shah, N., Bezzo, F., 2009. Spatially explicit static model for the strategic design of future bioethanol production systems. 1. Cost minimization. Energy Fuels 23 (10), 5121–5133. Zandi Atashbar, N., Labadie, N., Prins, C., 2017. Modelling and optimisation of biomass supply chains: a review. Int. J. Prod. Res. 1–25. https://doi.org/10.1080/ 00207543.2017.1343506.

567