Designing a Reliable and Congested Multi-Modal Facility Location Problem for Biofuel Supply Chain Network

This study presents a mathematical model that designs a reliable multi-modal transportation network for a biofuel supply chain system while site-dependent facility failure and congestion are taken into consideration. The proposed model locates the multi-modal facilities and biorefineries and determines the optimal production, storage, and routing plans in such a way that the overall system cost is minimized. We propose a hybrid Constraint generation-based Rolling horizon algorithm to solve this challenging NP-hard problem. The performance of this algorithm is tested in a example case study with numerical analysis showing that the hybrid algorithm can find near-optimal solutions to large-scale problem instances in a reasonable amount of time. Results indicate that the effect of congestion reduces the usage of multi-modal facilities in the biofuel supply chain network while bio-refineries and multi-modal facilities tend to move away from coastal areas when disruption probabilities are taken into consideration.


Introduction
The biofuel industry in the United States is growing at a fast pace.Bioethanol production increased from 1.15 billion gallons in 2000 to 15.08 billion gallons in 2017 [1].Simultaneously, biodiesel production is increasing sharply, producing almost 107 million gallons only in September, 2015.The new US Renewable Fuels Standard (RFS) sets an annual production target of 36 billion gallons of biofuel, consisting of mainly ethanol and biodiesel, by 2022 [2].Designing a robust biofuel supply chain network is of special importance, given this projected increase in production.The efficiency and reliability of a biofuel supply chain network will be determined by the performance of an integrated biofuel production and transportation system that not only operates well under normal conditions but also hedges against risk of disruption and/or congestion [3,4].
Designing an efficient and robust biofuel supply chain is challenging due to the physical characteristics of biomass (e.g., seasonality, uncertainty, and widely dispersed supply) [5].For instance, corn stover is only available from early September to November, while woody biomass is available all year round with the exception of three months in the winter.Overall, early September through late November is the peak biomass production season in the US, and more biomass is expected to ship through multi-modal facilities during those time periods of the year.This peak seasonality not only congests the facilities, but also impacts overall biomass supply chain efficiency [6].
Moreover, there is a significant increase in grain transportation for U.S. railways during the last few years.Grain transportation increased by 26 percent in 2014 compared to 2013 and 25 percent compared to the average of the three years (2011-2013) [7].This increase in grain transportation has a ripple effect on the overall supply chain structure, resulting in increased road traffic and air pollution.Another detrimental effect of congestion is that the quality of biomass degrades with time and may become unusable if the time needed to reach destination points increases [6].Therefore, congestion must be taken into consideration when designing a biofuel supply chain.
Congestion management has received attention from the research community over the last few years.Grove and O'Kelly [8] first investigated the relationship between hub-and-spoke networks and congestion by simulating the daily operations of a single assignment hub-and-spoke network.Elhedhli and Hu [9] introduced a non-linear congestion term in the objective function of an uncapacitated hub location problem which they later extended to a capacitated hub-and-spoke network design problem [10].Most recently, Camargo et al. [11] modeled the congestion as a convex cost function and developed a mixed integer nonlinear programming model for a multiple allocation hub-and-spoke network design problem.The authors successfully solved as many as 81 nodes by using a generalized Benders decomposition algorithm.A handful of studies [12,13] focused on the impact of congestion under demand uncertainty; however, this area remains vastly uninvestigated by the biofuel supply chain research community where feedstock seasonality and uncertainty can cause severe network congestion.Bai et al. [14] first introduced the traffic congestion factor in the traditional facility location model to decide the optimal location of refineries and the flow of biomass and ethanol in the transportation network.Later, Hajibabai and Ouyang [15] proposed an integrated mathematical model that focused on minimizing the total cost of facility construction, roadway capacity expansion, and biomass/biofuel transportation delay due to congestion.Few studies elaborated on the significance of congestion pricing on various aspects of decision making.Interested readers can refer to the studies by De Palma and Lindsay [16] and Friesz et al. [17] for more details.
Similar to facility congestion, transportation infrastructure, particularly those bearing multi-modal traffic, may be vulnerable to various disruption risks, such as natural disasters (e.g., 2005 Hurricane Katrina, 2008 China and 2009 Haiti earthquake [18,19]) and human-caused disasters (e.g., 2003 US Northeast blackout, 2010 Gulf of Mexico oil spill [20,21]).Furthermore, some areas are recognized as disaster prone areas.For instance, Figure 1 shows that in total 39 different storms affected North Carolina between 2000 and 2008 [22].Daskin [23] was the first to consider facility unavailability in a maximal covering location problem.The author then extended this study to formulate a maximum expected covering location problem where the demand is not readily met as the facilities capable of fulfilling the demand is busy in serving other customers [24].This work was further extended by Drezner [25] to develop a reliable p-median facility location problem.Snyder and Daskin [26] proposed an integer programming model for the stochastic fixed charge p-median location problem assuming that the facility disruptions occur independently and with equal probability.Later, Cui et al. [27] extended this work by creating both discrete and continuous models where disruption probabilities are dependent on specific sites.Location design models explored by An et al. [28] for several types of emergency service facilities under probabilistic facility disruption risks to reduce the expected cost of losing service.Parvaresh et al. [29] developed a bi-level game model by formulating a multiple allocation p-hub median problem under intentional disruptions.The primary objective of the study was to determine the hubs, failure of which would diminish service efficiency to the greatest extent.Darayi et al. [30] integrated a multi-commodity network flow formulation with an economic interdependency model to quantify the multi-industry impacts of a disruption in a transportation network.The goal of the study was to measure the vulnerability of a multimodal freight transportation network and assess the importance of network components using the vulnerability analysis.
In the biofuel supply chain literature, the number of studies that consider the effect of biorefinery disruption are scarce.Li et al. [31], Wang and Ouyang [32], and Huang and Pang [33] studied failure risks at biorefineries but ignored the fact that disruption can simultaneously occur at multi-modal facilities.Please note that the biofuel supply chain system from biomass production to bio-refineries can be viewed as a hub-and-spoke transportation network.Please note that most of the prior studies (e.g., [34][35][36][37][38][39][40]) assume that the transportation hubs are always functioning and will never fail, which cannot adequately describe the real-world scenarios.A few studies focus on managing and rescheduling rail and port operations during different disrupted scenarios (e.g., [41][42][43]).A system dynamics approach was proposed by Peng et al. [44] to analyze the behavior of disrupted disaster relief supply chain by simulating the uncertainties associated with predicting post-seismic road network and delayed information.Kim and O'Kelly [45] proposed a single and multiple allocation reliable p-hub location model where each arc and each hub is assigned a reliability factor.Marufuzzaman et al. [46] and Marufuzzaman and Eksioglu [47] proposed an optimization model to examine the impact of disruption risks at multi-modal facilities and biorefineries in a biofuel supply chain network.Poudel et al. [48] proposed a pre-disaster planning model while considering the failure probability of links between multi-modal facilities.Bai et al. [49] investigated the impact of failure risks at biorefineries and at intermediate transportation hubs in the biofuel supply chain network.Through experimentation it was observed that the number of biorefineries was directly proportional to the probability of disruption.Authors also found that the cost of a biorefinery had significant impact on the optimal decisions.
To fill this gap in the literature, this paper studies the impact of facility disruption along with network congestion in the context of a biofuel supply chain network.We propose an optimization model that determines the optimal location of biorefineries and multi-modal facilities while simultaneously determining the production, storing, and routing plans in such a way that the overall system cost is minimized under normal and disrupted scenarios.To tackle the complexity of this problem, we develop a hybrid decomposition algorithm that combines Constraint generation algorithm with a Rolling horizon heuristic.The proposed algorithm is validated by developing a example case study.We use the states of Mississippi and Alabama as a testing ground for this study.The outcomes of this study provide several managerial insights such as the optimal deployment of multi-modal facilities and biorefineries, amount of biomass transported via multi-modal facilities/highways, and the impact of disruption and congestion on biofuel supply chain performance.
The paper is organized as follows: Section 2 formulates the mathematical model; Section 3 introduces the solution algorithms; Section 4 presents numerical results and draw managerial insights; and finally, Section 5 provides conclusions and future research directions.

Problem Description and Model Formulation
The main objective of this paper is to build a reliable biofuel supply chain network that takes multi-modal facility disruption along with facility congestion into account.The aim is to build biorefineries and use multi-modal facilities among candidate facility locations to produce, store, and transport biomass/biofuel in such a way that the overall system cost is minimized under both normal and disrupted scenarios.Additionally, our model attempts to minimize facility congestion caused by feedstock seasonality and facility unavailability.Figure 2 presents the structure of the supply chain network consisting of two biomass supplier, three multi-modal facilities, two biorefineries, and three markets for biofuel production.

Nonlinear Problem Formulation
Let us denote a biofuel supply chain network G(N , A), where N denotes the set of nodes and A denotes the set of arcs.Set N consists of the set of harvesting sites I, the set of candidate multi-modal facilities J , the set of biorefineries K, and the set of markets for biofuel G i.e., N = I ∪ J ∪ K ∪ G. Similarly, set A can be partitioned into four disjoint subsets i.e., A = A 1 ∪ A 2 ∪ A 3 ∪ A 4 where, A 1 represents the set of arcs connecting harvesting sites I with multi-modal facilities J , A 2 represents the set of arcs connecting multi-modal facilities J with biorefineries K, A 3 represents the set of arcs that directly connect harvesting sites I with biorefineries K, and A 4 represents the set of arcs connecting biorefineries K with markets G.
Transportation distance along arcs A 1 are relatively short; therefore, trucks are preferred to ship biomass along these arcs and its unit transportation cost is denoted by c ijt .On the other hand, transportation distance along with shipping volumes in arcs (j, k) ∈ A 2 are high.Thus, transportation mode such as rail can be used to ship biomass along these arcs.Cargo containers are usually used to transport biomass between the multi-modal facilities, and the loading and unloading of biomass in the cargo container incurs a fixed cost of ξ jkt .The unit cost along arcs (j, k) ∈ A 2 are denoted by c jkt .Therefore, the unit transportation cost along arcs {(i, j), (j, k)} can be represented by c ijkt =c ijt +c jkt .We further denote c ikt as the unit truck transportation cost of biomass along arcs (i, k) ∈ A 3 .These arcs are preferred when the harvesting sites i ∈ I are located close to the biorefineries k ∈ K and thus, direct shipments of biomass using trucks are considered cheaper over multi-modal facilities.Finally, trucks are used to ship biofuel along arcs A 4 by incurring a unit cost of c kgt .
Let us consider l ∈ L be the set of capacities for multi-modal facilities and biorefineries and T be the set of time periods.In this model, we denote s it as the amount of biomass available at site i ∈ I in time t ∈ T and b gt as the demand for biofuel at market g ∈ G in time period t ∈ T .The unmet biofuel demand at market g ∈ G in time period t ∈ T can be substituted by paying a unit penalty cost of π gt .If the unit cost of delivering biofuel exceeds this threshold cost (π gt ), it will be beneficial to use substitution rather than producing biofuel at a higher cost.Locating a biorefinery of capacity l ∈ L at each location k ∈ K incurs a fixed set up cost Ψ lk .Similarly, using a multi-modal facility of capacity l ∈ L at each location j ∈ J in time period t ∈ T costs η ljt .Some other key parameters used in the model formulation are: φ as the conversion rate (ton/gallon) from raw biomass to biofuel; biofuel inventory holding cost h kt , and biofuel production cost p lkt .
We assume that each facility j ∈ J K can disrupt independently in a given time period t ∈ T .The corresponding site-dependent failure probability is denoted by {q jt } j∈J K,t∈T .When either multimodel facility or biorefinery disrupts, we assume that the biomass will be routed through an emergency carrier by paying an additional penalty cost which is much higher than the regular transportation cost.To capture this case, we assume that the emergency carrier will cost β (β > 1) times more than the regularly scheduled cost c ijkt .Tables 1 and 2 summarize the notation that we use in this model.We now define the following location and allocation decision variables for our optimization problem.The second set of variables Z := {Z jkt } j∈J ,k∈K,t∈T identifies the number of containers used between multi-modal facilities j and k in period t.The remaining set of decision variables determine how to route, produce, and store the biomass in the biofuel supply chain network.Let X := {X ijkt } i∈I,j∈J ,k∈K,t∈T ∪ {X ikt } i∈I,k∈K,t∈T ∪ {X kgt } k∈K,g∈G,t∈T denote the flow of biomass through the multi-modal facilities and highways to biorefineries, and biofuel from biorefineries to the customer demand points.Variables P := {P lkt } l∈L,k∈K,t∈T represent the amount of biofuel produced at biorefinery k of size l in period t.Variables H := {H kt } k∈K,t∈T represent the amount of biomass stored at biorefinery k in period t, and U := {U gt } g∈G,t∈K represents the amount of unsatisfied demand in market g in period t.Multi-modal facilities are congested during the peak harvesting seasons of biomass (early September till late November).This congestion in multi-modal facility leads to a dramatic increase in the total supply chain cost.The impact of congestion on total cost becomes more severe when the total flow of feedstock X ijkt approaches the capacity c cap lj of a multi-modal facility j ∈ J which is modeled here as M/M/1 queuing system.Under steady-state conditions, the system-wide average waiting time for the entire network can be represented as [10].Here, for a facility j ∈ J , when the amount of feedstock X ijkt increases, the ratio of the equation will increase exponentially which will address the impact of congestion realistically in the model formulation.Now, considering c 0 as the congestion factor, the system-wide congestion cost becomes ∑ j∈J ∑ t∈T c 0 The objective function of [DR] minimizes the total expected system cost under both normal and disrupted conditions.More specifically, the first, second, and third terms represent respectively the total set-up cost of locating biorefineries, usage cost of multi-modal facilities, and the fixed cost associated with transporting cargo containers between the multi-modal facilities.The fourth and fifth terms represent the production and inventory holding cost at the biorefinieries.The sixth term is the regular transportation cost, which is weighted by (1 − q jt )(1 − q kt ), the probability that both facilities operate normally along each route (i, j, k).When either or both multi-modal facility j ∈ J and biorefinery k ∈ K disrupt, which occurs at a probability of (q jt + q kt − q jt q kt ), biomass X ijkt will flow through that route at a higher variable cost βc ijkt .This is reflected by the seventh term of the objective function.The eighth term is regular transportation cost along the arc A 3 and A 4 .The ninth and tenth terms represent the congestion cost and the penalty cost for the biomass supply shortage.

Symbol Description η ljt
Fixed cost of using a multi-modal of capacity l ∈ L at location j ∈ J in time period t ∈ T Fixed cost of a cargo container for transporting biomass along arc (j, Biomass storage/handling capacity of an multi-modal facility of size l ∈ L at location j ∈ J Failure probability of multi-modal facility j ∈ J K in time period t ∈ T Constraints (1) indicate that the amount of biomass shipped from supplier i ∈ I in period t ∈ T is limited by its availability.Constraints (2) are the flow conservation constraints at biorefineries.Constraints (3) indicate that the amount of biofuel delivered to the market is limited by the amount produced in period t ∈ T .Constraints (4) indicate that the demand for biofuel can be fulfilled either through the distribution network or through substitute products available in the market.Constraints (5) indicate that the total amount of biomass shipped through a multi-modal facility is limited by its capacity c cap lj .Constraints (6) count the number of containers needed for shipping biomass on each arc (j, k) ∈ A 2 .Constraints (7) set biofuel production capacity (p cap lk ) limitations at a biorefinery.Constraints (8) set biomass storage limitations (h cap lk ) at a biorefinery.Constraints ( 9) and (10) ensure that, at most one multi-modal facility/biorefinery of capacity l ∈ L is operating in a given location j ∈ J K. Finally, constraints (11) and ( 12) are the binary constraints, (13) are the integer constraints, and ( 14) are the standard non-negativity constraints.

Model Linearization
The objective function of model [DR] contains a nonlinear congestion cost function.We use the technique proposed by Elhedhli and Wu [10] to linearize this congestion term.To do this, let us now consider a new decision variable M := {M jt } j∈J ,t∈T as follows.
Equation ( 15) can be further reduced as follows: We now introduce another continuous variable R := {R ljt } l∈L,j∈J ,t∈T as follows: Given that ∑ l∈L Y ljt = 1, the above equation becomes In the condition when Y ljt = 0 constraints (17) forces R ljt = 0, which is enforced by using an additional constraints 0 Proof.While differentiating the function R ljt w.r.t.M jt , we get the first derivative, The first derivative is positive while the second derivative is negative and it proves that the function Lemma 1 shows that function R ljt (M jt ) is concave and can be approximated by a set of tangent cutting planes as shown below [50] which is equivalent to where {M h jt } j∈J ,t∈T ,h∈H are the set of points used to approximate Equation (18).The value of {Y ljt } l∈L,j∈J ,t∈T is finite; therefore, the value that variable {M jt } j∈J ,t∈T provides is also finite.This implies that the set H should be finite.Equation ( 19) is now derived from ( 17) and ( 18) as follows: The linearized objective function Subject to (1)-( 4), ( 6)-( 13), (19), and

Solution Methodology
This section discusses solution techniques for solving the linearized model [LDR].Our resulting model [LDR] is an extension of the fixed charge, uncapacitated facility location problem which is known to be an N P-hard problem [51].Therefore, commercial solvers, such as CPLEX, fail to solve the large scale instances of this problem.In this section, we propose solution techniques such as Constraint generation algorithm, Rolling horizon algorithm, and an integration of Constraint generation with Rolling horizon to solve model [LDR].The aim is to produce a high-quality solution in solving model [LDR] in a reasonable amount of time.

Constraint Generation Algorithm
In model [LDR], (19) generate large number of constraints.Therefore, it is really challenging to solve model [LDR] while considering all constraints at once.This motivates us to develop a Constraint Generation (CG) algorithm that can efficiently solve model [LDR] despite generating a large number of constraints through (19).The algorithm, extensively studied in [52], proceeds by sequentially solving a series of integer programs with a subset of the constraints obtained from (19) and added thereafter as needed.
The process stops when the algorithm finds a solution for the sub-problem which does not violate any constraints within some accepted tolerance in the full problem [DR].Otherwise, a new set of points and thus a new set of constraints/cuts are generated which are added to [LDR] in the next iteration.
Let LB q and UB q denote a lower and upper bound of the original problem at iteration q.We further let v[LDR] be the solution of the objective function value of [LDR] and (Y q , Z q , P q , H q , U q ) be its optimal solution.The following proposition provides the lower bound of the CG algorithm.Proposition 1.For any given subset of points {M h jt } H q ⊂H (23) provides a lower bound of the optimal objective function value of [DR].
Proof.[LDR(H q )] is the relaxed version of problem [DR].Thus, the optimal objective function value v[LDR(H q )] obtained from Equation (23) provides the lower bound to the optimal objective value of [LDR].Since problem [LDR] is an approximation for problem [DR], solution v[LDR(H q )] will also provide a valid lower bound for optimal objective function value of [DR].Proposition 2. For any given subset of points {M h jt } H q ⊂H (24) provides an upper bound for the optimal objective function value of [DR].
Proof.Any solution feasible to [LDR(H q )] also provides a feasible solution to [DR] since all the constraints of [DR] are also contained by [LDR(H q )].Thus, the objective function value of [DR] evaluated at (Y q , Z q , P q , H q , U q ) as shown in Equation ( 24), provides an upper bound for the optimal objective value of [DR].
The algorithm continues until the gap between the lower and upper bound falls below a tolerance level ; otherwise, a new set of points {M h new jt } are generated using the current solution (shown below) and the process continues.The overall algorithm is shown in Algorithm 1.
q ← 1; UB q−1 ← +∞; LB q−1 ← −∞ Choose an initial set of points {M h jt } h∈H q to approximate function ] and obtain its optimal solution (Y q , Z q , X q , P q , H q , M q , U q ) Update the lower bound: LB q ← [LDR(H q )] Update the upper bound UB ← min{UB q−1 , UB(Y q , Z q , X q , P q , H q , M q , U q )} using Equation (24) Generate a new set of points

Rolling Horizon Heuristics
Solving [LDR(H q )] is still considered challenging since the problem is a special case of a capacitated facility location problem which is known to be an N P-hard problem [51].Therefore, in this section we introduce a heuristic approach based on the rolling horizon scheme suggested by Balasubramanian and Grossmann(2004) [53] and Kostina et al. (2011) [54].This approach decomposes problem [LDR(H q )] into a series of smaller subproblems with a few consecutive time periods and solves each subproblem with initial set of points {M h jt } j∈J ,t∈T ,h∈H .The algorithm terminates once all the subproblems are investigated.The overall algorithm is shown in Algorithm 2.

Let t 0
s denote the starting time for subproblem s and N s be the total number of time periods on subproblem s.We can set a fixed or variable size of N s for each subproblem.Each approximate subproblem of the Rolling horizon algorithm is denoted by [LDR(H q )(s)].Now, for each iteration approximate subproblems are solved by setting variables as: The values of Y ljt and Z jkt for t < t 0 s are fixed to the values found when solving approximate subproblem s + 1, and the steps involved on the iteration are as follows.

Step 1:
Initialize starting time period to t 0 s , length of time interval to N; set s → 1.

Constraint Generation with Rolling Horizon Heuristic
We now propose a hybrid algorithm which integrates the Rolling horizon approach within the Constraint generation framework.Solving the subproblems [LDR(H q )] of the Constraint generation algorithm is still considered challenging.The motivates us to employ Rolling horizon algorithm to solve the subproblems [LDR(H q )] of the Constraint generation algorithm.The aim is to provide high quality feasible solution for subproblem [LDR(H q )] in a reasonable amount of time.The steps involved in solving model [DR] using the proposed Constraint Generation with Rolling Horizon Heuristic is shown in Algorithm 3.

Algorithm 3 Constraint Generation with Rolling Horizon Heuristic.
q ← 1; UB q−1 ← +∞; LB q−1 ← −∞ ,t r 0 = 0, r ← 1, N r , terminate ← false Choose an initial set of points {M h jt } h∈H q to approximate function M h jt 1+M h jt while (UB q−1 − LB q−1 )/UB q−1 > ξ do while (terminate = false) do Set: Fixing the values of {Y ljt } l∈L,j∈J ,t∈T and {Z jkt } j∈J ,k∈K,t∈T for t < t 0 s end if r ← r + 1 end while Obtain optimal solution of [LDR(q)] (Y q , Z q , X q , P q , H q , M q , U q ) Update the lower bound: LB q ← [LDR(H q )] Update the upper bound UB ← min{UB q−1 , UB(Y q , Z q , X q , P q , H q , M q , U q )} using Equation (24) Generate a new set of points

Case Study
This section presents a computational study on model [DR] to test our algorithms and to draw managerial insights.We use the state of Mississippi and Alabama as a testing ground for our study.The case study demonstrates how the model and the solution algorithms proposed in this study can be applied to a real-world supply chain network design problem.

Data Description
Biomass Supply and Demand: Biomass availability (i.e., corn stover and woody biomass) in Mississippi and Alabama are obtained from the Knowledge Discovery Framework (KDF) database of United States Department of Energy [55].Every year these two states produce about 6.85 million tons (MT) of biomass from 99 different counties.The counties that produce more than 10,000 tons of biomass each year are only considered in this study.Figure 3a shows the distribution of the densified biomass available in this region.
We set the total biofuel demand as 500 million gallons per year (MGY) for the planning horizon [7].We have identified 149 counties from the region as the demand area and the centroid of those counties is considered as the demand point.
Investment Costs: In this study we consider a total of 53 potential multi-modal facilities and 43 potential biorefinery locations.We used existing multi-modal and inter-modal facilities as potential multi-modal facilities.The potential biorefinery locations are selected based on the prior literatures applied to our test region ( [37,38]).Figure 3b shows the locations and distribution of potential multi-modal facilities and biorefineries in this region.The annualized fixed cost for a rail ramp of capacity 1.05 million ton per year (MTY) is estimated to be $54,949 per year [56].Based on this information, we estimate the investment costs for other rail ramp capacities.We consider five different rail ramp capacities (l = 0.55 MTY, 0.65 MTY, 0.75 MTY, 0.85 MTY, and 1.00 MTY).The annualized fixed cost for facilities are considered based on a lifetime span of 30 years and assuming a discount factor of 10%.The annualized fixed costs for a biorefinery with a 45 million gallon per year (MGY) capacity is set to be $1.59 million [46,57].This cost is estimated based on a lifetime of 20 years, and a discount factor of 15% is assumed.We consider five different biorefinery sizes (l = 20 MGY, 40 MGY, 60 MGY, 100 MGY, and 150 MGY).Please note that the actual fixed cost may vary by location; however, we use a common fixed cost as a reasonable approximation.Transportation Costs: We consider that truck and rail are the two modes of transportation available for transporting biomass from multiple sources to destination points.Trucks are used to transport biomass from a feedstock supplier i ∈ I to a multi-modal facility j ∈ J .Trucks can further be used to transport biomass directly from a feedstock supplier i ∈ I to a biorefiney k ∈ K and biofuel from biorefinery k ∈ K to market g ∈ G. Table 3 presents the cost estimations and some of the technical factors which are considered in this study.The unit truck transportation (c ij ) costs are calculated based on the following Equation [58]: where t d represents distance-dependent transportation cost ($/mile/truckload) such as expenses including fuel, insurance, maintenance, and permitting costs, t t represents time-dependent cost ($/hr/truckload) including labor and capital costs.Both costs t d and t t vary with traveling distance d ij .Please note that the loading and unloading cost γ is treated as fixed transportation cost and does not depend on traveling distance.All of these costs are estimated for a semi-truck with a load capacity δ cap 1 = 18 tons [59,60] and average traveling speed s 1 = 40 miles per hour is assumed.In case of disrupted scenarios, we multiply a penalty term (i.e., 1.5) with c ij in order to calculate the unit cost for truck transportation (c ik ).We also calculate the unit transportation cost (c kg ) of biofuel from biorefinery k ∈ K to market g ∈ G. Finally, rail transportation cost is obtained from a study by Gonzales et al. [61] where the authors set $2248 as the fixed shipment cost per rail car and $1.12 as the unit variable transportation cost per mile traveled.We have used Arc GIS Desktop 10 to present the existing transportation network in this region.The network includes railways as well as local roads, urban roads, and major highways in the states of Mississippi and Alabama.Estimating Disruption Probabilities: Failure probabilities of multi-modal facilities are obtained from a study by Marufuzzaman et al. [46].The authors consider three major types of natural disasters i.e., hurricane, floods, and droughts that frequently impact the Southeast region of United States.Figure 4 demonstrates the estimated disruption probabilities of the candidate multi-modal facilities where the size of the circle is proportional to its disruption probability.Results show that the failure probability of multi-modal facilities are not uniform throughout the year.For instance, late August until late September is the peak hurricane season in this region.Therefore, the multi-modal facilities located near the coastal regions of Mississippi and Alabama are more prone to disaster during this time period of the year.

Analyzing the Impacts of Congestion Cost on System's Performance
Figure 5 shows the impact of congestion cost (c 0 ) on a biofuel supply chain network performance.It is observed that the number of multi-modal facilities, containers transported between the facilities, and total amount routed through highways are significantly impacted by the changes in congestion cost.For instance, the number of multi-modal facilities used decrease by 12%, 22%, and 43% if congestion cost (c 0 ) increases from $0 to $10 thousand, $100 thousand, and $1 million, respectively.Please note that no multi-modal facilities are used if the congestion cost increases to $10 million.We observe the similar trend for the number of containers transported between the multi-modal facilities and biorefineries.Experimental results indicate that the number of containers transported between the multi-modal facilities and biorefineries decreases by 75% for an increase in congestion cost from $10 thousand to $1 million (shown in Figure 5b).Moreover, it is observed that the number of containers used during peak season (i.e., (t 5 − t 7 )) decreases sharply with an increase in congestion cost.This in turn increases the biomass transportation through highways.As it can be observed from Figure 5c that as the congestion cost increases, routing more biomass through highways are preferred compared to rail transportation.

Impacts of Different Risk Levels on System's Performance
This section presents the impact of different risk levels (i.e., different value of {q jt } j∈J K,t∈T ) on the biofuel supply chain network performance.To evaluate the impact of risk levels, we conducted experiments considering the original disruption probabilities from the study by Marufuzzaman et al. [46] for the risk averse case and disregarded the disruption probability for the risk-neutral case on different biomass available seasons (i.e., peak and low production season).In the experiment we set t = 1 for the month of May.Corn stover is typically harvested during September until November in a given year.Forest residues are harvested all year around, except during the winter (December to February) due to humid weather.Thus, September to November are considered as the peak biomass production season for this region since both feedstock types are available during this time period of the year.We consider March to mid-August as a low biomass production season since only woody biomass is available during that time period of the year.Figure 6 depicts the risk-neutral scenario of the biofuel supply chain network under various biomass available seasons: peak and low biomass production season.It is obvious from Figure 6 that when the decision maker adopts a risk-neutral attitude, bio-refineries away from the multi-modal facilities are selected.Additionally, it is important to note that as there is no extra cost in the risk-neutral scenario, the model selects bio-refineries away from the multi-modal facilities.Results indicate that at the peak biomass production season more multi-modal facilities are open than low biomass production season (shown in Figure 6a,b).More specifically, the model opens three additional multi-modal facilities in Lawrence, Lamar, and Pike county to cope with high feedstock availability in the peak biomass production season.
Figure 7 demonstrates the risk-averse scenario of the biofuel supply chain network under peak and low biomass production seasons.It is clear from Figure 7 that when the decision maker adopts a risk averse attitude, the model decides to select the bio-refineries in northwest Mississippi that are close to their multi-modal facilities.The low disaster prone area (Figure 4) and high availability of biomass (Figure 3a) are the main reasons for selecting bio-refineries and multi-modal facilities in the northwest.Moreover, in the low biomass production season the model decides to close one multi-modal facility in Newton county (Figure 7b).
It is clear from Figures 6 and 7 that when the decision maker adopts the risk averse scenario the amount of unsatisfied demand U gt increases.The model decides to fulfill some portion of the demand using substitute products rather than establishing biorefinaries near the higher disaster prone areas.We observe that at risk averse scenario, the amount of unsatisfied demand U gt increases by 17.5% and 15.1% during peak season (t = 7) and low season (t = 12), respectively, over the risk neutral scenario.In summary, these results indicate that adopting different risk levels has impact on the biofuel supply chain under peak and low biomass production seasons.

Impact of Different Disruption Scenarios on System Performance
To quantify the benefits of designing a reliable and congested biofuel supply chain system, we create different realistic tornado scenarios (shown in Figure 8).We assume that the tornado hits in the coastal region of Mississippi and Alabama.The intensity of color determines the impact of the tornado and higher intensity represents higher chance of facility failure in that zone.We assume that the multi-modal facilities and biorefineries located within the affected region will fail and cannot be used for production/distribution for the entire season.The length and area of the affected region will vary based on the strength and path of the tornado strikes.This experiment sets the length at 200 miles on a straight path by assuming that a medium tornado strikes the coast as illustrated in Figure 8.
We now solve our model

Impacts of Supply Change on System's Performance
Figure 9 shows the impact of supply changes on the biofuel supply chain network performance.It is observed that the number of multi-modal facilities, the containers transported between the multi-modal facilities and biorefineries, and total amount routed through highways are significantly impacted by the biomass supply changes.For instance, while considering the facility congestion cost (c 0 ) of $20,000 in the modeling process, a 20% increase in supply increases the total number of multi-modal facilities by 30.61% while a 20% decrease in supply results in a 28.95% decrease of the total number of multi-modal facilities used by the biomass supply chain netwrok.A similar trend can be observed for the number of containers transported between the multi-modal facilities.A 20% increase in supply results in a 20.14% increase in the total number of containers while a 20% decrease in supply results in a 22.7% decrease in total number of containers transported between the multi-modal facilities and biorefineries.The results further indicate that the biomass transportation via highway increases with an increase in biomass supply.Moreover, it is observed that during the peak biomass production season (t 5 − t 7 ), the total number of containers and the amount routed through highways increases sharply from previous time periods.10a).However, if failure probabilities are taken into account, biorefineries tend to be located further away from the coasts and become more centralized around the multi-modal hubs (shown in Figure 10b).Moreover, if congestion of multi-modal facilities is taken into consideration, it is observed that the number of multi-modal facilities decrease while increasing their capacity (shown in Figure 10c).When both failure probabilities and effect of congestion are taken into consideration, then most of the multi-modal facilities do not change their location but biorefineries shift their location away from coastal regions (shown in Figure 10d).
The above experiments have a clear impact on the unit delivery cost for biofuel.The unit cost of biofuel is calculated as $4.41/gallon while failure probability and congestion of multi-modal facilities are not taken into consideration.However, the unit cost increases to $4.49/gallon if congestion is considered and $4.58/gallon if failure probability of facilities are taken into consideration.When both the effect of congestion as well as probability of disruption are taken into account, unit cost jumps to $4.74/gallon.Results indicate that a supply chain structure robust against congestion and disruption incurs 7.3% additional unit cost of biofuel than if neither were considered while designing the supply chain network.Including the effect of disruption in the model increases unit cost more than that of congestion.Specifically, unit cost of biofuel increases 3.85% when disruption probability is considered, while this value is only 1.82% when congestion is considered.This is a clear indication that more investment is required to design a reliable biofuel network which is not vulnerable to sudden disruption than designing a network not susceptible to congestion of facilities.

Analyzing the Performance of the Solution Algorithm
We now present our computational performance in solving model [LDR] using the algorithms proposed in Section 3. The algorithms are terminated when at least one of the following criteria is met: (a) the gap between the upper and lower bound falls below a threshold limit , i.e., = |UB − LB|/UB × 100%, (b) the maximum time limit is reached i.e., t max = 3600 s, and (c) the maximum iteration limit is reached i.e., n = 500.We

Conclusions
This study investigates the impact of facility congestion and disruption in the context of a biofuel supply chain network.A mixed-integer non-linear programming model is developed and later linearized to determine the optimal location of facilities (i.e., multi-modal, biorefineries) for production and storage and shipment routes of biomass in such a way that the cost of congestion and disruption are minimized.A hybrid Constraint generation coupled with Rolling horizon algorithm is proposed to solve the optimization problem in a reasonable amount of time.Computational results indicate that the model locates biorefineries at areas with low disruption risks.Moreover, as the disruption probability becomes high, the model tends to locate fewer facilities and satisfy the remaining demand either by trucks or purchasing from outside market.Similarly, with the increase in congestion cost, the number of multi-modal facilities as well as the total number of containers used to transport biomass decreases and impacts the overall supply chain activities.
In summary, the major contributions of this study to the existing literature are manifold.First, our optimization approach considers both facilities disruption along with associated congestion under the same decision making framework.Second, we propose customized solution approaches to solve our optimization model.Finally, the modeling results are experimentally validated by developing a example case study.The findings can be used by decision makers to design and manage a reliable, non-congested biofuel supply chain network.
This research opens up several future research opportunities.It will be interesting to investigate how our model behaves under different system uncertainties (e.g., biomass supply and demand uncertainty, technology uncertainty).Further, in addition to considering facility uncertainty, the transportation link that connects multiple facilities may also be disrupted and congested during the time of a natural catastrophe (e.g., hurricane, tornado).These issues will be addressed in future studies.

Figure 2 .
Figure 2. Network representation of biofuel supply chain.

.
This non-linear cost function is now added in the objective function to account for the congestion cost.The following is an Mixed-Integer Nonlinear Programming (MINLP) formulation of the problem referred to as model [DR].
Biomass storage capacity of a biofuel plant of size l ∈ L at location k ∈ K v cap Cargo container capacity p cap lk Production capacity of a biofuel plant of size l ∈ L at location k ∈ K φ Conversion rate from biomass to biofuel c 0 Congestion factor q jt

(a) Supply sites 1 (Figure 3 .
Figure 3. Biomass distribution and potential location of multi-modal facilities and bio-refineries.

Figure 6 .
Figure 6.Network representation under risk neutral scenario.

Figure 8 .
Figure 8. Network representation of different disruption scenarios on system performance.

Figure 9 .
Figure 9. Impact of supply changes on system performance.(a) Multi-modal facilities; (b) Containers transported; (c) Biomass shipped via highways.

2 (Figure 10 .
Figure 10.Network representation of minimum cost, reliable, congested, and reliable and congested model.
generated various problem sizes by varying the number of supply sites | I |, potential multi-modal facilities | J |, potential biorefineries | K |, total customers | G |, and the length of the planning horizon | T |.

Table 1 .
Description of the sets.{Y ljt } l∈L,j∈J ,t∈T ∪ {Y lk } l∈L,k∈K determine the size and location to use/open multi-modal facilities and biorefineries, i.e.,
1} andZ jkt ∈ Z + for t 0 s ≤ t ≤ t 0 s + N s 0 ≤ Y ljt ≤ 1 and Z jkt ∈ R + for t > t 0 s + N s Solve the approximate sub-problem [LDR(H q )(s)] using CPLEX if(t 0 > |T |) then stop ← true elseFixing the values of {Y ljt } l∈L,j∈J ,t∈T and {Z jkt } j∈J ,k∈K,t∈T for t < t 0

Table 3 .
Description of the sets and parameters.
Table 4gives an overall idea about the size and characteristics of the overall problem solved.Table 5 reports the error gap (Gap), and the running time (CPU) of the algorithms proposed in this study.Note that the boldface values in the table indicate the best solutions for that instance among the proposed approaches.We calculate the error gaps of Rolling Horizon [RH] algorithm and hybrid produces the reported solution 9.73% faster than algorithm [CG].Moreover, it is observed that a balance between faster computation time and solution quality can be achieved by solving problem [DR] using [CG+RH] algorithm.Experimental results indicate that algorithm [CG+RH] solves 22 out of 24 problem instances with less than 1% error gap, and on average provides solutions 46.62% faster than [CG] and 40.86% faster than [RH].
Constraint Generation with Rolling Horizon algorithm [CG+RH] using the lower bound obtained from the Constraint Generation [CG] algorithm.This is because algorithms [RH] and [CG+RH] only generate upper bounds for model [DR].Results indicate that the [RH] algorithm provides a solution with less than 1.0% error gap for only 5 out of 24 problem instances.On the other hand, [CG] solves 19 out of 24 problem instances with less than 1.0% error gap.However, algorithm [RH]

Table 4 .
Problem size of the test instances.

Table 5 .
Comparison between different solution approaches.