A multi-stage stochastic programming approach for blood supply chain planning

Computers & Industrial Engineering 122 (2018) 1–14

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

A multi-stage stochastic programming approach for blood supply chain planning

T

⁎

B. Zahiria, S. Ali Torabib, , M. Mohammadic, M. Aghabegloob a

Department of Operations Management & Strategy, School of Management, State University of New York at Buffalo, Buffalo, NY, USA School of Industrial Engineering, College of Engineering, University of Tehran, Tehran, Iran c IMT Atlantique, Lab-STICC, UBL, F-29238 Brest, France b

A R T I C LE I N FO

A B S T R A C T

Keywords: Blood supply chain planning Vehicle routing problem Multi-stage stochastic programming Meta-heuristic algorithms

Perishability of blood products and uncertainty in donation and demand sizes complicate the blood supply chain planning. This paper presents a novel bi-objective mixed-integer model for integrated collection, production/ screening, distribution and routing planning of blood products, and seeks to simultaneously optimize the total cost and freshness of transported blood products to hospitals. To cope with inherent uncertainty of input data, a multi-stage stochastic programming approach with a combined scenario tree is presented. Due to the high complexity of the problem, a novel hybrid multi-objective self-adaptive differential evolution algorithm is developed, which benefits from the variable neighborhood search with fuzzy dominance sorting (hereafter it is briefly called MSDV). MSDV is validated through comparing its performance with two of the most common multi-objective evolutionary algorithms (i.e. MOICA and NSGA-II). Applicability of the proposed decision model is also tested through a real case study. Our results show that the solution efficiency of a network can be balanced with its effectiveness through customer satisfaction. Further, several sensitivity analyses are carried out to provide valuable managerial insights.

1. Introduction During the last 50 years, health spending has risen dramatically. In the U.S. for example, healthcare spending accounted for 5% of the gross domestic product (GDP) in 1960 which rose to 14.1% of GDP by the year 2001 (Margaret, Brandeau, & Pierskalla, 2004). In 2010, it accounted for 17.9% of GDP (The World Bank, 2015) which is expected to climb to near one fifth of the U.S. economy in 2021 according to the updated estimates from the federal government (Kaiser Health News, 2015). In Iran, the health services accounted for 5.6% GDP by the year 2010 and it is increasing likewise. The growing number of aging population, increasing costs of prescription drugs and emerging new health-related technologies and procedures are the main contributors to the escalating healthcare costs (Margaret et al., 2004). Owing to the fact that in many countries public resources are not always adequate for affording healthcare costs, governments and healthcare providers must find effective (in terms of meeting patients' requirements) and efficient (in terms of the lowest cost) ways of delivering health services. In this paper, we focus on a subset of healthcare industry (i.e. a blood supply chain), and optimization of the related network. Blood supply chain, which deals with time-sensitive products and in particular

perishable products, pose specific challenges. By definition, a perishable product has a limited lifetime, during which it can be used, and after which it should be discarded (Nagurney, Masoumi, & Yu, 2012). Derivation of blood, which is mainly referred to the whole blood (WB) results in three main products: red blood cells (RBC), platelets (PLT) and plasma (PLS), which have specific characteristics and usages (Katsaliaki, 2008). PLT’s function is to stop the bleeding by binding together when damage is recognized in blood vessels. This product is regarded as the most perishable and expensive blood component and is expired in 5–7 days. Owing to this fact, making a large inventory of PLT to meet the demand in far later succeeding periods is rather impossible. RBC is the most demanded component; whose shelf life is between 35 and 42 days. Lastly, PLS is the largest component of blood and has a long shelf life. Blood donation rate is an uncertain variable and quite difficult to predict. In the UK, about 60% of the population could be considered as the potential donors, while just 6% of total population donates blood. In the U.S., even though blood donation is voluntarily and purchased commercially, only 5% of population donates blood and only 8% of donors donate regularly, and approximately 62% of them never return

⁎

Corresponding author. E-mail addresses: behzadza@buffalo.edu (B. Zahiri), [email protected] (S.A. Torabi), [email protected] (M. Mohammadi), [email protected] (M. Aghabegloo). https://doi.org/10.1016/j.cie.2018.05.041 Received 15 February 2017; Received in revised form 27 September 2017; Accepted 25 May 2018 Available online 26 May 2018 0360-8352/ © 2018 Elsevier Ltd. All rights reserved.

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

decision support system in which a non-convex integer model is used for platelet production and blood mobile scheduling in blood centers. Jabbarzadeh, Fahimnia, and Seuring (2014) proposed a robust mathematical model for blood supply chain network design for post-disaster situations. The model tries to specify the number and location of both fixed and temporary blood facilities as well as coordinating supply with demand at the minimum cost for a set of disaster scenarios. Gunpinar and Centeno (2015) presented several district and stochastic models to minimize total cost, shortage and wastage level of RBCs and platelets in hospitals. Abdulwahab and Wahab (2014) applied approximate dynamic programming for modeling blood platelet inventory management. Four measures are considered in the stochastic model including shortages, outdating, inventory level, and reward gained. Zahiri, Torabi, Mousazadeh, and Mansouri (2015) presented a mixed integer linear programming model to address the strategic and tactical decisions in a blood collection and distribution network design problem considering both temporary and permanent collection sites. Furthermore, a robust possibilistic programming approach was applied to cope with inherent epistemic uncertainty in input parameters. Şahinyazan, Kara, and Taner (2015) addressed a new selective vehicle routing problem for bloodmobiles and shuttles to determine their tours and daily locations. The paper presents two mathematical models. The first model finds the best tours for shuttles and bloodmobiles while minimizing the logistics costs of collecting a pre-determined blood level. Also, the second model maximizes the total amount of blood that can be collected with a given number of vehicles during the planning horizon. Osorio, Brailsford, and Smith (2015) addressed a structured review of the literature on quantitative modeling for the blood supply chains. According to their study, majority of the literature in this field is focused on an individual echelon and does not consider relationships between the different echelons of a blood supply chain. In order to have a structured gap analysis and recognizing a number of future research directions in the context of blood supply chains, Table 1 provides a coding system for classifying the relevant literature according to eight major characteristics. Special features of the relevant

to a blood center. In some countries like Japan, blood is even being imported (Drackley, Newbold, Paez, & Heddle, 2012; Beliën & Forcé, 2012). As the aging population grows and new medical procedures need more blood components, the demand for blood products in certain countries is increasing considerably. Moreover, the increasing number of deferrals from donor pool exacerbates the difficulties of blood supply (Van Der Bom et al., 2009). In this paper, we develop a novel comprehensive mathematical model for integrated blood collection, production, distribution and routing planning (BCPDRP) in blood centers, considering perishability of the products. 2. Literature review Study on blood supply chain peaked in 1976–1985 and dropped until the turn of the century. Most of the recent studies focus on RBCs and WB (Katsaliaki, 2008). A simulation model was proposed for RBC procurement by Madden, Murphy, and Custer (2007). They considered both double RBC and WB collection and the impact of travel deferral on units available for transfusion. Heddle et al. (2009) studied the disposition data of RBC’s inventory of several hospitals using a benchmarking tool to identify those factors influencing the red blood outdating. Prastacos and Brodheim (1980) provided a decision support system for a regional blood bank in order to specify the minimum achievable outdate and shortage targets, required distribution policy to meet the targets, and supply level for achieving alternative targets. Mustafee, Taylor, Katsaliaki, and Brailsford (2009) proposed a discrete event simulation model to assess blood ordering policies in the UK national blood service. Stanger, Yates, Wilding, and Cotton (2012) studied blood inventory management in the UK hospitals to determine critical factors to improve inventory management of perishable items. Stanger et al. (2012) studied RBC inventory management in transfusion laboratories to develop the principles, which ensure the minimum spoilage rate. Some other recent papers which have studied red blood cells and/or whole blood supply chains could be traced in (Seifried et al., 2011; Eyles, Webert, Arnold, & McCurdy, 2011; Hussein & Teruya, 2012; Kamp, Heiden, Henseler, & Seitz, 2010). PLT is another blood product, which has been studied in this context. Katsaliaki (2008) conducted a specific analysis using the data related to the UK national blood service and developed a simulation model for the UK’s supply chain of RBC and PLT to determine those management policies which are more practical and cost effective. Rytilä and Spens (2006) proposed a computer simulation model to increase the blood supply chain’s efficiency of Finish healthcare system. Haijema, van der Wal, and van Dijk (2007) proposed a combined Markov dynamic programming and simulation approach for production management and inventory optimization of PLT and applied it to a Dutch blood bank. Haijema, van Dijk, van der Wal, and Smit Sibinga (2009) studied the production and inventory management of PLT in breaks such as Easter and Christmas. Van Dijk, Haijema, Van Der Wal, and Sibinga (2009) combined stochastic dynamic programming and simulation to optimize PLT production using a five-step procedure. Zhou, Leung, and Pierskalla (2011) studied PLT inventory management to find the optimal policy that minimizes the expected cost. Other papers addressing the PLT inventory issues are (Stanger et al., 2012; Chapman, 2007; Blake, 2009; Perera, Hyam, Taylor, & Chapman, 2009). Among the mathematical programming based models, Nagurney et al. (2012) developed an optimization model for procurement, testing, processing, and distribution of blood using the variational inequality approach. Hemmelmayr, Doerner, Hartl, and Savelsbergh (2010) developed an integer programming model to find delivery routes of blood products from blood center(s) to hospitals. The model accounts for natural uncertainty of blood usage in hospitals. A variable neighborhood search (VNS) solution method was also devised to find optimal or near-optimal solutions. Ghandforoush and Sen (2010) developed a

Table 1 Classification of blood supply chains’ literature.

2

Blood type Platelet Plasma Red blood cell Whole blood Planning horizon

PLT PLS RBC WB

Operational Tactical

OPR TAC

Strategic Modeling/Method Simulation Dynamic programming Stochastic dynamic programming Mixed integer nonlinear programming Mixed integer linear programming Integer nonlinear programming Discrete event simulation Statistical analysis Evaluation Benchmarking Robust Outputs Number of Production/ Collection centers Inventory level Location/allocation Ordering policy

STR

UR A SR DD DR BT

SM DP SDP

Usage rate Availability percentage Supply rules Donor deferral Delivery routes Amount of blood transfusion Different blood risks Factors for red blood cells’ outdating Objectives Min cost Min distance Min delivery time Min max unmet demand

MINLP

Min coefficient of variance

MV

MILP

Min risk

MR

INLP DES S E B RB

MO MRS

NP

Min outdated bloods Max remaining shelf life Type of Uncertainty Stochastic Robust Number of Time Periods Single period Multi period

SG M

IL L OM

Solution Method Exact Heuristic/Meta-heuristic

Ex He

DBR FRBCO

MC MD MDT MMU

ST R

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

Table 2 A summary of the reviewed articles in the context of blood supply chains. Reference articles

Blood type

Planning horizon

Modeling/ Method

Output

Objective function

Type of Uncertainty

Time Period

Solution Method

Haijema et al. (2009) Van Dijk et al. (2009) Ghandforoush and Sen (2010) Şahin, Süral, and Meral (2007) Hemmelmayr et al. (2010) Vera Hemmelmayr, Hartl, and Savelsbergh (2009) Sha and Huang (2012) Cetin and Sarul (2009) Drackley et al. (2012) Kamp et al. (2010) Eyles et al. (2011) Madden et al. (2007) Beliën and Forcé (2012) Chapman (2007) Perera et al. (2009) Stanger et al. (2012)

PLT PLT PLT WB WB WB

OPR TAC OPR STR OPR OPR

SM, SDP SM, SDP MINLP MILP MILP MILP

NP IL, NP NP, L L DR DR

MC

ST ST

SG SG SG SG SG M

He Ex Ex Ex He He

WB WB RBC, WB RBC RBC RBC, WB WB RBC RBC, PLT RBC, PLT, WB RBC, PLT WB PLT, RBC PLT WB PLT, PLS, RBC, WB WB PLT, PLS, RBC, WB

TAC STR OPR OPR STR STR STR STR STR OPR

MINLP BNLGP DES SM S SM S E S S

L, NP L IL, OM, OP BT DBR A DD

MC MC, MD, MV

M SG SG SG

He He Ex Ex

SG

Ex

STR STR OPR OPR STR, TAC STR, TAC

LP RB MILP SDP MINLP, RB MILP

NP, L NP, L, IL OM IL L L

MC, MR MC MC MO OM MC, MMU

ST ST R R

SG M M M M M

He Ex Ex He Ex Ex

OPR OPR, TAC

MILP MINLP

DR NP, IL, OM, DR

MC MC, MRS

ST

SG M

Ex He

Nagurney et al. (2012) Jabbarzadeh et al. (2014) Gunpinar and Centeno (2015) Abdulwahab and Wahab (2014) Zahiri et al. (2015) Zahiri and Pishvaee (2017) Şahinyazan et al. (2015) Our paper

papers in the literature are provided in Table 2. As shown in Table 2, majority of the reviewed articles are dedicated to a specific product (mostly RBC and WB), while the multi-product aspect of such chains has remained quite unexplored in the literature. Moreover, covering all major planning-related processes in a blood supply chain through an integrated mathematical model over a mixed operational and tactical planning horizon has not been addressed yet. Such an integrated model will provide the decision makers (DMs) with the global optimal solution as all interconnections between the major processes are considered concurrently within a single comprehensive decision model. In addition, coping with the inherent uncertainty of input parameters has not widely been considered in the literature. Accordingly, the main contributions of this paper that differentiate it from the other relevant papers are listed as follows:

ST

ST

Fig. 1. Graphical statement of the problem.

commonly taking place via bloodmobiles (BMs) or directly at the main (fixed) centers (MCs). Compared to the MCs, BMs are more extensive in both quantity and diversity. Consequently, their main responsibility is the local collection of blood. Since they have quite limited capacities, at the end of each period (mainly day), all the collected blood packs are transferred to the MCs. Noteworthy, establishing an MC is significantly more expensive than a BM. Therefore, utilizing BMs to cover sporadic demands seems justifiable in practice. All the received blood packs are sent to lab centers, where different components are derived afterwards. Therefore, the concerned network is a single-product network before the MCs, while thereafter it is transformed into a multi-product chain. As the network includes multiple MCs, in order to meet the dynamic demand more effectively, transshipment between MCs are allowed. More details about the advantages of transshipment in multi-site supply chains can be found in Torabi and Moghaddam (2012) and Mousazadeh, Torabi, and Zahiri (2015). Finally, all the derived components are sent to the demand points (i.e. hospitals/clinics). Notably, transfusion of blood products requires several time consuming processes. Since the corresponding time significantly varies among recipients, a congestion in this situation is very common (more information is provided in Section 3.1). In this study, based on the emergency level of recipients, each recipient is associated with a priority, and is transfused according to him/her priority. Since the

• Developing an integrated multi-product BCPDRP model considering perishability of products. • Developing a multi-priority queuing system to model the congestion of blood products at each hospital. • Formulating a bi-objective model to optimize the total costs and the •

MC MD MC MDT

remaining shelf life of products simultaneously, while making a possible trade-off analysis between these objectives. Developing a novel hybrid meta-heuristic algorithm to efficiently solve the large size instances.

The remainder of this paper is organized as follows. In Section 3, the problem statement and model formulation are provided. Section 4 is dedicated to the proposed solution approach. In Section 5, the developed metaheuristic algorithm is elaborated. Computational experiments are given in Section 6. Finally, Section 7 represents concluding remarks and some avenues for further research. 3. Problem statement Fig. 1 shows a schematic view of a general blood supply chain network and all interactions between its elements. Donation (D) is 3

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

transfusion takes a significant time for each recipient, waiting times of recipients should be taken into consideration in order to maximize their satisfaction level. On the other hand, as the waiting time for transfusion directly affects the survival chance of the recipients as well as the perishability of the derived products, controlling this time can significantly elevate the efficiency and effectiveness of the proposed network. Accordingly, the proposed decision model seeks to minimize the total cost of the blood supply chain network including the transportation, inventory holding, production and imposed cost by patients’ waiting times in queue as the cost efficiency measure of the model, while maximizing the freshness of transported blood products to the hospitals as the efficacy measure. In this section, we present a generalized mathematical model for collection, production, distribution and routing planning problem for blood components. In doing so, we incorporate the three following critical measures in our formulation:

Dnit

Dn′jt vcpv ctijp

crhl cm jj′ cppi

cwhpn Cost associated with waiting of product type p for emergency level n at hospital h ghpnt Unit penalty cost of unfulfilled demand of product type p with emergency level n at hospital h at time period t mjp Unit inventory cost of blood component type p in MC j up Maximum shelf life of blood component type p B Number of hospitals and MCs

• Reliability: the blood network should provide fresh products to recipients, in order to reduce the discarded products. • Efficiency: The total cost of the network should be minimized. • Efficacy: The network should be effective from the recipients’ point

Decision variables: x ijt Amount of WB transported from BM i to MC j at time period t Amount of WB which is required from direct donation in yjt MC j at time period t wjhpnt Amount of blood component p transported from MC j to node h for emergency level n at time period t Invjprt Stock level of blood component p at time period t that is produced in MC j at time period r rhlvt 1; if hospital h precedes hospital l immediately in the route of vehicle v at time period t, 0; otherwise ajh 1; if hospital h is assigned to MC j, 0; otherwise vm jj′ prt Amount of blood component p produced in MC j at time period r and transshipped to MC j′ at time period t ⩾ r um jj′ pt Amount of blood component p transshipped from MC j to MC j′ at time period t ejhprt Amount of blood component p produced in MC j at time period r and transported to hospital h at time period t ⩾ r Sub-tour elimination variable in the route of vehicle v Mhvt serving hospital h at time period t Amount of blood component p that is loaded on vehicle v in qjpvt MC j at time period t Amount of blood component p that is unloaded in hospital ′ qhlpvt h on the way to hospital l in the route of vehicle v at time period t t Arrival rate of product type p with emergency level n to λhpn hospital h at time period t t Waiting time in queue of product type p with emergency Wqhpn level n at hospital h at time period t t Total waiting time in queue of product type p with Whpn emergency level n at hospital h at time period t

of view as well. Noteworthy, in order to elevate the patients’ satisfaction level (in light of the aforementioned shortcomings), their waiting times in the queue should be incorporated through a suitable queuing system.

The following assumptions are considered to formulate the problem under investigation:

• Both MCs and BMs are capacitated • Based on the emergency level of recipients, they are classified into three priority groups namely; Low, Normal, and Urgent Each unit of unfulfilled demand in each period incurs penalty cost • • Dynamic nature of the supply and demand side data necessitate considering a dynamic planning horizon including multiple periods • Perishability of the blood products is taken into account by considering the production time of each unit • Transshipment between MCs are allowed in order to meet the •

fluctuating and dynamic demands through required redistribution of inventories in the whole system Products’ deliveries to the recipients as well as the transshipment flows between MCs are done via optimal routes in each period

To formulate the problem mathematically, the following indices, parameters, and decision variables are introduced. Indices and sets: i Index of BMs, i ∈ I j, j′ Index of MCs, j ∈ J Index of hospitals, h, l ∈ H h, l M Aggregate set of hospitals and MCs (h ∪ j ) p Index of blood components (i.e. blood product types), p ∈ P t Index of time periods, t ∈ T r Index of production periods, r ∈ T v Index of available vehicles, v ∈ V n Index of emergency level of recipients, n ∈ N

3.1. Priority queuing system Based on a report by California Pacific Medical Center (CPMC), transfusion of a blood unit usually takes two hours. Therefore, high number of recipients, congestion and long shortage of blood products at hospitals may increase the risk for both products’ quality and recipients. To cope with this challenge, the derived products are prioritized based on different recipients’ emergency levels. Accordingly, a priority queuing system is developed to model the arrival flow of products to hospitals. To do so, an Exponential distribution is considered for modeling the inter-arrival times of products arriving at recipients. Moreover, it is assumed that during the peak hours, the average arrival

Parameters: dhnpt Demand of blood component p with emergency level n at hospital h at time period t capi Storage capacity of BM i (in terms of WB unit) ′ Storage capacity of MC j for blood component p capjp

lup

Maximum number of receivable blood packs by BM i at time period t Maximum number of receivable blood packs by MC j through direct donation at time period t Capacity of vehicle v for blood component p Unit transportation cost of blood component p from BM i to MC j Unit transportation cost of products from hospital h to hospital l Unit transshipment cost of products from MC j to MC j′ Unit production cost of product type p at BM i

Production rate of blood component p per WB unit 4

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

problem in the considered blood network.

and service rates are both constant and the arrival of products to hospitals follows the Poisson distribution. This allows us to model the queue formed by derived products entering to each hospital as an M/M/ c queuing system. Accordingly, hospitals act as multiple servers with infinite capacities in this queuing system, thus no balking of arrivals is considered. In the proposed system, the priority of each recipient is addressed based on its numbering index. That is, recipients with priority index n = 1 have the highest priority, and lower priorities are numbered in an ascending order. Therefore, a priority queuing system is developed to model the arrival flow of blood products to the hospitals. Based on the above notations, the arrival rate of products to hospital h is calculated by: t λhpn =

∑

Min∑

t t Whpn = Wqhpn +

t μhpn

,

= 1−

∀ h, p, n, t , (3) Ψh Ψh − 1

∑ e=0

t

t

e

⎡ (λhpn / μhpn ) ⎤ ⎢ ⎥ e! ⎣ ⎦

∑

dhnpt crhl rhlvt

h∈H l∈H v∈V t∈T

cm jj′ um jj′ pt +

∑ ∑ ∑ ∑ h∈H

p∈P t∈T

t t cwhpn Whpn λhpn

p∈P n∈N t∈T

max ∑

∑ ∑ ∑

i∈I

p∈P r∈T t∈T

1 ⎛ r + up−t × e ⎞ ⎜ jhprt ⎟ u ∑ p ( ⎝ ⎠ n ∈ N dhnpt )

(9)

(∑n ∈ N dhnpt )

∀ h, p, n, t , (4)

t ∑n′= 1 λ hpn ′

∑

t Ψhμhpn

x ijt ⩽ capi ,

∀ i ∈ I, t ∈ T,

x ijt ⩽ Dnit ,

∀ i ∈ I, t ∈ T,

(10)

j∈J

,

∀ h, p, n, t , (5)

∑

= 1. Furthermore, by substituting Eqs. (4) and (5) in Eqs. where (3), the queue waiting time of arrived products in hospitals is presented by:

⎣

t ⎤⎛ ⎜ Ψhμhpn ⎥ 1−

t ∑n n′= 1 λ hpn′

t e t Ψh − 1 ⎡ (λhpn / μhpn ) ⎤

∑e = 0

⎢ ⎣

n−1

e!

,

∀ j ∈ J, t ∈ T,

(12)

∀ h, p, n, t ,

⎥ ⎦

t−1

∑

t

λ ∑ ⎞ ⎜⎛1− n′= 1 hpn′ ⎟⎞ t Ψhμhp ,n−1 ⎠⎝ ⎠

⎛ ∑ xijt + yjt ⎞⎟ + ⎜ ⎠ ⎝ j∈J ∀ j ∈ J , p ∈ P, t ∈ T ,

Invjprt − 1 + lup

r=0

⎟

t Ψhμhpn

⎦⎝

(11)

j∈J

yjt ⩽ Dn′jt ,

1 Ψh t ⎡Ψ ! (Ψ μ t −λ t ) ⎜⎛ μhpn ⎟⎞ t ⎢ h h hpn hpn λhpn ⎝ ⎠

+

p∈P r∈T t∈T

normalize the product flow based on the demand. Therefore, the objective function tries to maximize the defined ratio (i.e. the product flow over the demand) for each demand zone and time period.

t Bhp 0

t Wqhpn =

j∈J

The second objective (9) targets to maximize the remaining life of the delivered blood products to the hospitals. Remaining shelf life of a product is expressed by r + up−t , where r is the production date, up is the shelf life of the product p, and t accounts for the current time period (it also represents the date when the blood product is transported to a hospital since ejhprt exists in the formulation). Each type of blood products has its own shelf life. Therefore, the remaining shelf life in the objective function (9) is divided by the corresponding initial shelf life to normalize the impact of shelf life variation between different products. 1 Noteworthy, is incorporated in the formulation in order to

n

t Bhpn

∑

The first objective (8) minimizes the total cost of the network, which includes production cost; transportation cost of blood products from the BMs to the MCs; cost of failure to meet demand, inventory cost, routing and transshipment cost between the MCs, respectively.

(2)

t ⎡ ⎛ μhpn ⎞ t t t Ahpn = ⎢Ψh ! (Ψhμhpn −λhpn )⎜ t ⎟ λ ⎢ ⎝ hpn ⎠ ⎣

⎤ t + Ψhμhpn ⎥, ⎥ ⎦

+

∑ ∑ ∑

ctij. x ijt

j∈J t∈T

⎞ ⎛ ghpnt dhnpt − ∑ wjhpnt ⎟ ⎜ n∈N p∈P t∈T j∈J ⎠ ⎝ ∑ ∑ ∑ mjp Invjprt + ∑ ∑ ∑

∀ h, p, n, t ,

1 , t t t Ahpn × Bhpn × Bhp ,n−1

i∈I

(8)

In order to calculate the queue waiting time of arrived products in hospitals, Eqs. (3) are presented, in which Athpn and Bthpn are calculated as Eqs. (4) and (5), respectively (Zahiri, Tavakkoli-Moghaddam, Mohammadi, & Jula, 2014). t Wqhpn =

∑∑ ∑

∑ ∑ ∑ ∑

j , j′ ∈ J

Alternatively, the total service times, which is the sum of waiting times of blood products in the queue and their operating times in labs, is calculated by:

1

cppi (x ijt + yjt ) +

p∈P t∈T

+

+

(1)

j∈J

j∈J

h∈H

∀ h, p, n, t ,

wjhpnt

∑ ∑ ∑

i∈I

cap′jp ,

∑

um j′ jpt ⩽

j′

(6)

Constraints (10)–(13) are capacity constraints. Constraint set (10) ensures that the storage level of each BM at each period should be lower than its relative capacity. Constraint set (11) implies that the amount of blood products transferred from each bloodmobile to the MCs at each period are less than the potential units of donations. Similarly, constraint set (12) states that the amount of blood components at each MC should be less than the potential units of direct donations. Constraint set (13) indicates the maximum capacity of each MC.

t Finally, the total service time, Whpn , is computed by:

1

t Whpn =

t

⎡Ψ ! (Ψ μ t −λ t ) ⎛ μhpn ⎞ t ⎢ h h hpn hpn λhpn ⎝ ⎠ ⎣ ⎜

n

⎟

Ψh

Ψ −1 ∑e =h 0 ⎡ ⎢ ⎣ n−1

t

t / μ t )e (λhpn hpn

⎤ ⎥ ⎦

t

λ ∑n′= 1 λ hpn′ ∑ t ⎤⎛ ⎞ ⎛⎜1− n′= 1 hpn′ ⎞⎟ + Ψhμhpn t t ⎥ 1− Ψhμhpn Ψhμhp , n−1 ⎠⎝ ⎠ ⎦⎝ 1 + t , ∀ k, g, p, t , μhpn ⎜

e!

(13)

⎟

t

(7)

∑ r=0

3.2. The proposed mathematical model

⎛ ∑ x ijt + yjt ⎞⎟− ∑ ∑ wjhpnt− ∑ um jj′ pt ⎜ j′ ⎠ n∈N h∈H ⎝ j∈J + ∑ um j′ jpt , ∀ j ∈ J , p ∈ P , h ∈ H , r ∈ T , t = 1,

Invjprt = lup

j′

In this section, a bi-objective mathematical model based on the aforementioned assumptions is formulated for the integrated BCPDRP

(14) 5

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al. t−1

t

∑

Invjprt =

r=0

∑

⎛ ∑ xijt + yjt ⎞⎟− ∑ ∑ wjhpnt ⎜ ⎠ n∈N h∈H ⎝ j∈J + ∑ um j′ jpt , ∀ j ∈ J , p ∈ P , h ∈ H , r

⎛ ∑ xijt + yjt ⎞⎟− ∑ ejhprt− ∑ vm jj′ prt ⎜ j′ ⎠ h∈H ⎝ j∈J + ∑ vm j′ jprt , ∀ r , t ∈ T , p ∈ P , j ∈ J , r = t ,

Invjprt − 1 + lup

r=0

− ∑ um jj′ pt j′

Invjprt = lup

j′

j′

∈ T , 1 < t < up ,

(15)

(21)

∑

Invjprt = Invjprt − 1−

ejhprt − ∑ vm jj′ prt + j′

h∈H t−1

t

∑

Invjprt =

r = t + 1 − up

∑ r = t + 1 − up

j′

∑

(16)

Eqs. (14)–(16) are inventory balance constraints. Since blood components are perishable, tracking the production time of blood components’ inventories is inevitable (Kallrath, 2005). In this way, an index which determines the period during which the product is produced (r), is applied for inventory variables (Invjprt ). In fact, inventory level of each blood type at each period only consists of those products that are not perished. Eq. (14) determines the inventory level at the end of the first period. Owing to the assumption that there is no initial inventory, it deals with blood component production through the first period and the amount of products transported to MC during the first period. Transshipments between MCs are also considered in the formulation. Eq. (15) specifies the inventory levels for those periods smaller than the shelf life (t < up ). During the aforementioned period, it is still possible to have inventory of blood components produced at period r. Hence, all blood components’ inventories, regardless of production times, are t−1 calculated as (∑r = 0 Invjprt − 1). Eq. (16) ensures that for periods (t ⩾ up ), perished products are not figured in inventory calculation. Thus, only those blood components that are produced from period t + 1−up until the period t −1 can be regarded as the inventory of period t −1 (i.e., t−1 ∑r = t + 1 − up Invjprt − 1).

∑ ∑ n∈N

∑ ∑

wjhpnt ⩽

∑ ∑

∀ j ∈ J , p ∈ P , t < up ,

ejhprt ,

h∈H r=0

∑ ∑

wjhpnt ⩽

∑

∑ ∑

um jj′ pt ⩽

j ́ ∈J

vm jj′ prt ,

∑

∑

rhlvt −

l∈M

p ∈ P , j ∈ J , t < up ,

j ́ ∈J

j′ ∈ J r = t + 1 − up

vmjj ′ prt ,

rlhvt = 0,

∀ h ∈ M, v ∈ V , t ∈ T, (26)

l∈M

∑

rhlvt +

r jlvt −ajh ⩽ 1,

∀ h ∈ M −{j}, v ∈ V , t ∈ T , j ∈ J

l∈M

(27)

(19)

∑

(25)

Eq. (27) shows that a hospital or an MC can be assigned to a specified MC in each period, only if there is a route passed by the hospital or MC and originated from that specified MC.

∑

t

∑

∀ t ∈ T , h, l ∈ M , v ∈ V ,

Eq. (26) indicates that a vehicle should enter and leave a node at the same period.

(17)

(18)

j′ ∈ J r = 0

um jj′ pt ⩽

(24)

Eq. (25) is the sub-tour elimination constraint, which guarantees that each route is originated from an MC and terminated at the same MC.

qjpvt =

v∈V

∑

∀ h ∈ M, t ∈ T,

∀ j ∈ J , p ∈ P , t ⩾ up ,

h ∈ H r = t + 1 − up

∑ ∑

∀ h ∈ H , p ∈ P, t ∈ T , (23)

Mhvt −Mlvt + (B × rhlvt ) ⩽ B−1,

t

∑

rhlvt = 1,

l∈M

ejhprt ,

dhnpt ,

n∈N

Eq. (24) ensures that each hospital is visited by only one vehicle at each period.

∑

∑

∑

wjhpnt ⩽

j∈J

l∈M v∈V

t n∈N h∈H

(22)

Constraint (23) states that the number of each type of blood components transported from MCs to each hospital should be lower than or equal to the demand of each product.

t n∈N h∈H

∀ r, t ∈ T , p

Eqs. (21) and (22) determine blood components’ inventory levels based on their production times. For periods equal to production time t = r , the blood component’s inventory level is equal to the production level minus the transported level to hospitals and the tradeoffs (i.e. net incoming transshipment flows) between MCs (see Eq. (21)). If inventory level is calculated for those periods different from the production time r, it will depend on the inventory of the previous period and the amount of blood components delivered to hospitals and transshipped between MCs before they perish (t < up + r ).

j′

∈ H , r ∈ T , t ⩾ up ,

vm j′ jprt ,

j′

∈ P , j ∈ J , t −r < up,

⎛ ∑ xijt + yjt ⎞⎟− ∑ ∑ wjhpnt ⎜ ⎠ n∈N h∈H ⎝ j∈J um j′ jpt , ∀ j ∈ J , p ∈ P , h

Invjprt − 1 + lup

− ∑ um jj′ pt +

∑

∑ ∑

wjhpnt ,

∀ t ∈ T , j ∈ J , p ∈ P, (28)

n∈N h∈M

Eq. (28) states that the amount of product type p loaded on all vehicles should be equal to the amount of that product which is delivered to all hospitals.

∀ p ∈ P , j ∈ J , t ⩾ up , (20)

Eqs. (17) and (18) track the production time of blood components, which are delivered to hospitals. Similarly, Eqs. (19) and (20) track the production time of blood components, which are transshipped between MCs. All products, which are transported to hospitals and/or transshipped between MCs at period t must be equal to the sum of blood components produced in the previous periods which have not exceeded their shelf lives and transported in the same period. Similar to the inventory balance constraints, these four constraints are period-dependent. In fact, Eqs. (17) and (19) are valid for those periods smaller than the product shelf life and Eqs. (18) and (20) are valid for those periods greater than or equal to products’ shelf lives. Noteworthy, these constraints can be easily merged with inventory balance constraints, however in this paper they are appeared separately for the sake of clearance.

∑

qjpvt ⩽ vcpv,

∀ p ∈ P, v ∈ V , t ∈ T , (29)

j∈J

Eq. (29) ensures that amount of product type p loaded on a vehicle does not exceed the relative capacity.

∑ ∑ ∑ n∈N

wjhpnt =

j∈J h∈M

∑ ∑ ∑

′ rhlvt , qhlpvt

∀ p ∈ P, t ∈ T ,

h∈M l∈M v∈V

(30)

Eq. (30) guarantees that the amounts of products delivered to different demand points are equal to the amounts of products transported from the MCs to the hospitals.

qjpvt =

∑

∑

h∈M−j l∈M−j

′ rhlvt , qhlpvt

∀ j ∈ J , p ∈ P, v ∈ V , t ∈ T , (31)

Eq. (31) ensures that the amounts of products delivered to different 6

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

single-period environments, MSSP approach was introduced to deal with the stochastic data in a dynamic (i.e., multi-period) environment. Through MSSP, data uncertainty is presented through a scenario tree and the objective function is to minimize the total imposed risk to the sequence of decisions (Zahiri et al., 2017). Fig. 2 shows a schematic view of a 4-stage scenario tree with 40 nodes, which further denotes a discrete set of scenarios with their associated probabilities in a dynamic planning horizon. Just for simplifying the presentation and the model formulation, we proceed with the compact form of the presented model as follows:

hospitals are equal to the amounts of products loaded on that vehicle.

∑ ∑

′ rhlvt ⩽ qhlpvt

l∈M v∈V

∑

dhnpt + um jj′ pt ,

∀ p ∈ P , t ∈ T , h ∈ M , j , j′

n∈N

(32)

∈ J,

Eq. (32) states that the amounts of products delivered to each node should be less than or equal to the corresponding demand.

∑ n∈N

wjhpnt ⩽

∑

dhnpt + um jj′ pt ,

∀ p ∈ P , t ∈ T , h ∈ M , j , j′ ∈ J ,

n∈N

(33)

min

Eq. (33) ensures that the amount of blood products transported to each node does not exceed the corresponding demand (including the hospital demands and the amount of transshipment between MCs).

rhlvt , ajh ∈ {0, 1}.

(Cx t + Fyt ), max

t

∑

Gx t

t

s.t. Axt ⩽ L, Bxt = M , Eyt = 1, Nx t ⩾ Hx t , Ix t ⩽ Kx t , Lx t ⩽ Pyt , x t ⩾

′ ⩾ x ijpt , yjpt , wjhpnt , Invjprt , vm jj′ prt , um jj′ pt , ejhprt , qjpvt , qhlpvt 0; Mhvt free variable,

∑

0, yt ∈ {0, 1}. (34)

(36)

(35)

where C, F, G, A, L, B, M, E, N, H, I, K, L and P are the coefficient matrices, and x t and yt are time-dependent continuous and binary variables, respectively. Let us assume that L is a scenario-dependent parameter with a predefined scenario tree by which its estimated values under various scenarios are shown. In this case, the model is transformed as follows:

Finally, Eqs. (34) and (35) show the type of decision variables. In order to derive a valid solution for the problem, we must ensure to find the global optima for the model (8)–(35). Note that the objective function (8) and the constraint sets (30) and (31) are nonlinear. Consequently, to guarantee finding the global optimal solution by an MINLP solver (like BARON in GAMS), convexity of the objective function (8) as well as the constraint sets (30) and (31) are approved in the Supplementary Material file, Section A.

min

∑ ∑

ρ (n)(Cx t (n) + Fyt (n)), max

n ∈ Tree t ∈ tn

∑ ∑

Gρ (n) x t

n ∈ Tree t ∈ tn

s.t. Ax t (n) ⩽ L (n), Bx t (n) = M (n), Eyt (n) = 1, Nx t (n) ⩾

4. The proposed solution approach

Hx t (n), Ix t (n) ⩽ Kx t (n), Lx t (n) ⩽ Pyt (n), x t (n) ⩾ 0, yt (n) ∈ {0, 1}.

To solve the proposed multi-period MINLP problem under data uncertainty, a multi-stage stochastic programming (MSSP) approach is applied in this section.

(37)

where n denotes a node on the scenario tree, and ρ (n) is the occurrence probability of node n (i.e. ∑n ∈ Λ ρ (n) = 1, where set Λ shows a branch of the scenario tree involving node n). Now, assume that M is another scenario-dependent parameter with a different scenario tree to L. In order to have a single stochastic planning model that considers the uncertainty of L and M concurrently, those scenarios related to M are integrated with the scenario tree of L, which leads to a combined scenario tree. Fig. 3 shows a schematic view of a 3-stage combined scenario diagram, where scenarios of L are depicted using full lines, while those pertaining to M are shown as blue dotted lines. At each node of the L scenario tree, different states for M can take place, resulting to have 4! = 24 total scenarios for the combined tree. The model (37) with a combined scenario tree can be rewritten as below:

4.1. Multi-stage stochastic programming approach In the decision problems with dynamic structure, where time and uncertainty play the major roles, the DM should be able to adopt a decision policy that can respond to events as they unfold. These decisions mainly correspond to the time when the new information or adjustments (i.e. recourses) are available for the DM. Such decisions can respond to realized outcomes that are unknown a priori, and are modeled under MSSP approach (Birge & Louveaux, 2011; Zahiri, Torabi, & Tavakkoli-Moghaddam, 2017). Despite the two-stage stochastic programming (TSSP) approach which is merely adopted in

Fig. 2. A schematic view of a scenario tree with 27 scenarios.

7

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

main phase of SDE is started, which includes searching the solution space using the self-adaptive mutation operators and VNS algorithm. In the main phase of the proposed SDE, the first step is to generate a population of random solutions with size PopSize . The individuals are then ranked and the non-dominated ones are archived by applying the fuzzy dominance sorting (Kundu et al., 2011). 5.2. The proposed MSDV After archiving the non-dominated solutions, a neighborhood search structure (NSS) is performed to improve the quality of non-dominated solutions. In this paper, the NSS is applied by utilizing the variable neighborhood search (i.e. VNS) algorithm (Mladenović & Hansen, 1997). In recent years, VNS has increasingly gained lots of attention and a large number of successful applications have been reported (Hansen &

Fig. 3. A schematic view of a combined scenario tree.

min

∑ ∑

ρi (n)(Cx ti (n) + Fyti (n)), max

n ∈ Tree t ∈ tn

∑ ∑

Gρi (n) xti

n ∈ Tree t ∈ tn

s.t. Ax ti (n) ⩽ L (n), Bx ti (n) = M (n), Eyti (n) = 1, Nx ti (n) ⩾ Hxti (n), Ixti (n) ⩽ Kx ti (n), Lx ti (n) ⩽ Pyti (n), x ti (n) ⩾ 0, yti (n) ∈ {0, 1}.

(38)

Mladenović, 2001). Using two or more neighborhoods, instead of one, in its structure has differentiated the VNS from the most other local search heuristics. Particularly, it is based on the principle of systematic change of neighborhood during the search (Behnamian, Zandieh, & Ghomi, 2009). In addition, to avoid costing too much computational time, the best number of neighborhoods is often three (Mohammadi, Jolai, & Tavakkoli-Moghaddam, 2013), which is followed by our algorithm. The three neighborhoods employed in our algorithm are defined below:

where i denotes the scenarios of M. 5. The hybrid solution algorithm In this section, we propose a new multi-objective meta-heuristic algorithm based on a self-adaptive differential evolution (SDE) algorithm and variable neighborhood search (VNS), which is briefly called MSDV algorithm, wherein a fuzzy dominance sorting technique is employed to obtain the non-dominated solutions.

I. Swap: By this operator, the places of two randomly selected bits are exchanged. II. Reversion: Via this operator, a random part of a solution is selected and its permutation is reversed. III. Inversion: By this operator, one bit is chosen randomly and its value is replaced with a new random value.

5.1. The proposed SDE SDE algorithm is a new variant of DE algorithm, wherein a selfadaptive version of mutation operator is applied. In the literature, several mutation operators have been introduced while using all of them simultaneously in an algorithm extremely increases the computational time. Table 3 contains four different mutation operators wherein e v s are randomly selected solutions from the population and ebest is the best solution in the population. Therefore, many studies have employed only one mutation operator to search the solution space. However, in the proposed SDE algorithm, several mutation operators are simultaneously applied without increasing the computational time. SDE algorithm includes two phase as Initialization phase and Main phase. In the initialization phase, each mutation operator obtains a score if it could create a better solution at each iteration. After a predefined number of iterations in the initialization phase, a selection probability (SP) metric is calculated for each operator by dividing the obtained score by the number of iterations. It should be noted that the sum of all SPs is equal to 1. The Pseudo code of the initialization phase of the proposed SDE algorithm is shown in Fig. I in Section B of the Supplementary Material. When the initialization phase is finished, the

It should be noted that the VNS is performed on a portion of whole population (PE ). Furthermore, the number of outer and inner iterations of the proposed VNS are equal to ItVNSmax and nRepeat , respectively. This procedure is iterated until the population size exceeds a predefined upper limit (i.e. Popmax ). The population is again ranked in order to select the best Popmax individual solutions. This procedure is continued until the stopping criterion is met. The number of function calls-NFC (Črepinšek, Liu, & Mernik, 2012) has been considered as the stopping criterion in the proposed MSDV. Pseudo code of the SADE’s initialization phase is provided in Fig. 4. 6. Computational experiments In this section, applicability of the formulated model and the developed solution approach are validated. We also conduct some sensitivity analyses and provide managerial insights from the numerical results. In doing so, we first present the results of a number of randomly generated problem instances in different sizes to show the validity of our results as well as demonstrating the performance of the introduced solution algorithm compared to Multi-Objective Imperialist Competitive Algorithm (MOICA) and Non Dominated Sorting Genetic Algorithm (NSGA-II), and then focus on our case study in the next subsection. Tables I and II in Section C of the Supplementary Material show the

Table 3 Mutation operators. Mutation Strategy

Operator

MO MO MO MO

Mi Mi Mi Mi

#1 #2 #3 #4

= = = =

ebest + F (ev2−ev3) ev1 + F (ev2−ev3) ebest + F (ev1 + ev2−ev3−ev 4 ) ei + F (ev1 + ev2−ev3−ev 4 )

8

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

Fig. 4. The Pseudo code of fuzzy dominance sorting.

6.2. Discussion

ranges of randomly generated input parameters. Suppose that the demand scenario tree has three branches with values of 1300 with probability of 0.3, 1400 with possibility of 0.4 and 1700 with possibility of 0.3 (see Fig. II in Section C of the Supplementary Material).

The performance of the proposed MSDV algorithm compared to the exact results (i.e., GAMS outputs - can be accessed on https://www. gams.com/) as well as the results of two competing algorithms (i.e. MOICA and NSGA-II) under different metrics are evaluated in Table 4 and Tables III and IV in Section D of the Supplementary Material. As we can see, the results of GAMS are only available for the small- and medium-size instances, which further shows the NP-hard nature of the problem and limited memory of the GAMS software. QM values equal to 1 for the proposed MSDV, indicate that the Pareto solutions of MSDV completely dominate the solutions obtained by two other algorithms. They also show the superior performance of the proposed algorithm compared to the other two algorithms considering the higher values in DM and lower values in CM, SM and MID. Higher value of DM and lower value of SM show larger Pareto front of the proposed MSDV

6.1. Performance metrics The five commonly used metrics in the evolutionary multi-objective optimization literature are: (1) quality metric (QM), (2) convergence metric (CM), (3) divergence metric (DM), (4) spacing metric (SM), and (5) mean ideal distance (MID) metric, which represent both quantitative and qualitative comparisons with MOEAs (see Deb, 2001 for further information).

Table 4 Results of small sized instances. Problem No

1 2 3 4 5 6 7 8 9 10

MSDV

NSGA-II

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

1 1 1 1 1 1 1 1 1 1

2.024 2.112 1.913 1.798 2.009 1.937 2.193 2.059 2.051 2.191

1.173 1.148 1.246 1.06 1.359 1.183 1.332 1.533 1.3 1.659

0.14 0.252 0.251 0.212 0.226 0.266 0.367 0.218 0.227 0.432

0.526 0.443 0.387 0.387 0.493 0.359 0.32 0.394 0.243 0.347

43 61 88 93 97 98 105 105 110 111

0 0 0 0 0 0 0 0 0 0

2.857 3.268 2.399 2.887 3.26 3.104 3.553 3.146 2.831 3.319

0.678 0.769 0.887 0.949 1.152 1.041 1.069 0.847 1.215 1.324

0.286 0.442 0.432 0.332 0.473 0.457 0.663 0.332 0.398 0.628

0.676 0.651 0.741 0.596 0.702 0.525 0.504 0.682 0.556 0.655

163 187 206 244 269 298 313 348 360 397

MOICA

1 2 3 4 5 6 7 8 9 10

GAMS

QM

CM

DM

SM

MID

CPU (s)

QM

CM

DM

SM

MID

CPU (s)

0 0 0 0 0 0 0 0 0 0

2.803 3.537 2.245 2.829 3.112 2.659 3.321 3.428 3.651 3.82

0.785 0.474 0.793 1.01 1.523 0.842 0.926 0.87 0.999 0.943

0.508 0.615 0.527 0.581 0.853 0.583 0.769 0.486 0.744 0.61

0.73 0.82 0.778 1 1.026 0.878 0.61 0.932 0.827 0.893

173 177 188 191 203 210 253 271 309 321

– – – – – – – – – –

– – – – – – – – – –

0.423 0.665 0.395 0.501 0.735 0.632 0.716 0.514 0.58 0.738

0.655 0.503 0.524 0.563 0.604 0.745 0.544 0.606 0.585 0.641

0.831 0.809 0.762 0.858 0.656 0.853 0.884 0.802 0.677 0.873

73 116 289 501 778 1134 1898 3340 4133 7760

9

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

budget. II. Most of the delivered products were not fresh (most especially PLTs), which resulted in a large amount of waste. As discussed before, among different derivations of WB, PLT has the shortest lifespan (i.e. normally just 5–7 days), which makes the concerned inventory and distribution planning quite difficult. Based on the annual reports of the concerned hospitals, the average remaining lifetime of platelet is 2.5 days, whose discarded cost includes 4% of the total hospitals’ daily operational costs. III. Although the amount of products sufficed to meet the total demand, significant waiting times of recipients at the hospitals, and ignoring the priority level of recipients have highly decreased the customer satisfaction level. SBDC is currently serving 12 hospitals, and has 5 active BMs, by which the potential donations in different parts of the city and suburbs are covered. Since both the number of donors (i.e. supply size) and demand are uncertain and variable from one stage to another, only forecasting of donations and demands for the next two stages (which corresponds to 60 days in our setting) is reliable. The time unit is considered as daily. As discussed earlier, among the different product types which are derived in SBDC, PLT because of its shortest life time and higher demand is the most concerned product in the institution. In this respect, SBDC is assumed to have an unlimited storage capacity. This is a reasonable assumption since as previously mentioned, owing to the small shelf life and high production costs, it is not possible to build a large amount of inventory for PLT at any particular time unit. Consequently, the storage capacity is not an active limitation to the problem. In the beginning of each stage, the scenarios associated with the demand size are defined by the directors based on their observations and past experiences (see Tables V–VII in Section E of the Supplementary Material for the data corresponding the considered case study). At the first stage, three different levels for each scenario dependent parameter as the most optimistic, the most probable, and the most pessimistic values are set. Accordingly, at the next stage, three other values based on the previous branch are set each of which associating with a probability of occurrence. Furthermore, historical data in SBDC reveals that donation size has a stationary probability distribution. Therefore, even though it is set for three levels (i.e., the most pessimistic, the most optimistic, and the most likely values), the values remain the same during the planning horizon. It should be noted that since the variations in both demand and donation sizes within a month are negligible, the daily rates for donations and demands are computed as total demands / donations in the stage .

Fig. 5. CPU times vs. problem instances for competing algorithms.

wherein solutions are distributed uniformly across the Pareto front. This provides more flexibility for the decision maker. On the other hand, lower values for CM and MID indicate that the Pareto front is as close as the optimal solutions. Furthermore, the efficiency of the proposed algorithm can be detected considering its lower computational times to achieve nondominated fronts (NFs). The CPU time varies from 163 to 3502 s, 173 to 3382 from the smallest scale to the largest for NSGA-II and MOICA, respectively, while this value for MSDV is 43 to 1553, which is quite a significant improvement in the runtime. Fig. 5 further shows the schematic view of such comparison. Notably, all computations have been done on an operator with Core-i5, 2.5 GHz of CPU and 4 GB of RAM processor. Hereafter, the quality performance of the proposed algorithm is schematically compared with the other two algorithms and the exact Pareto fronts (PFs) of GAMS under the test problem #6 in Fig. 6. As the results show, the NFs of the MSDV dominate those of other algorithms. 6.3. Implementation and evaluation 6.3.1. Case study: specifications and challenges Founded in 1979, Sari Blood Donation Center (SBDC) is one of the main blood donation centers in Mazandaran province, Iran. Statistics in the year 2014 reveal that approximately 85% of donation attempts were successful (in terms of acceptable age of donors, density of the certain concentrations, etc.), which resulted in collecting more than 135,000 units of donated blood. Although the amount of collected blood sufficiently met the total demand, some critical issues still exist:

total number of days in that stage

The scenario trees of both demand (i.e. Fig. 7(a)) and donation amounts (i.e., Fig. 7(b)) are prepared by the SBDC directors. Each arc corresponds to each stage, which is associated with the forecasted value on that period and its occurrence probability. The combined scenario tree is provided in Fig. 8, which involves 27 scenarios. In order to see the behavior of the objective functions at their extreme points, we first assume that only the first objective function is important to the DM, and the model is thus transformed into a singleobjective one (we name this model as (i)). The optimal configuration of delivery-routing for t = 1 is illustrated in Fig. 9(a). With the total inventory level of 4320 and OFV1 = 3691092 , the delivery-route passes two routes 3−6−1−5−4 and 2−9−11−10−12−7−8. Now, let us assume that only the maximum freshness of products at the hospitals is important for the DM (this model is called as (ii)). Reconsidering the model with only the second objective function, the delivery-routing for t = 1 can be derived as Fig. 9(b). The total inventory level for this configuration with OFV2 = 11705.5 and the optimal route of 6−1, 5−4 , 3−9−2−8, and 7−11−10−12, is 7389. Based on the obtained results, the following conclusions can be made;

I. The current network is not cost efficient and in many cases, the annual cost reached or even exceeded the maximum available

Fig. 6. PF vs. NFs. 10

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

Fig. 7. Scenario tree of the demand size (a) and donation size (b).

Fig. 8. Combined scenario tree for demand and donation sizes.

Fig. 9. The optimal configuration for the case (i) (a) case (ii) (b), and The optimal configuration for the middle point of NF (c).

• Since the maximum freshness is the only concern in the case (ii), • •

To show how the second objective function changes the solution configuration of the first objective function, to pursue maximum freshness, the solution of the middle point of NF is derived. Fig. 9(c) illustrates the derived solution whose objective function values are OFV1 = 3781265 and OFV2 = 11313.2 . We can see a tradeoff between the two objectives, which results in three delivery-routes and higher levels of inventory levels.

compared to the derived configuration in the case (i), we see higher values in inventory levels and the number of passed routes. For any situations between these two extreme points, nodes 1−6, 4−5, and 11−10−12 are passed within the same route, which can provide valuable managerial insights for transportation/delivery organizations, whose major concerns are on their available budget/ resources. The maximum number of routes derived from the case (ii) is 4, which can be considered as the maximum number of routes, since it is independent of delivery cost.

6.4. Sensitivity analyses In this section, some sensitivity analyses based on the provided data 11

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

Fig. 10. OFV1 (a) and OFV2 (b) vs. alteration of Cap(i).

Fig. 11. OFV1 (a) and OFV2 (b) vs. alteration of tvc.

remains constant. It is worth noting that for Cap (i) ⩾ 1000 , the model is transformed to its uncapacitated counterpart, since both OFV1 and OFV2 values remain constant. The indicated state of SBDC (marked as triangle) shows the acceptable level of capacities. The same behaviors for the OFVs can be seen, when the capacity of vehicles for distribution of products are altered. Since the capacity levels of available vehicles are different, we define a new variable as follows;

tvc =

∑ v

vc v

(39)

where tvc denotes the total available capacity of vehicles. As Fig. 11(a) shows, with increasing the tvc, OFV1 decreases as the delivery can be done using lesser vehicles, while for lower levels of capacity, more vehicles and consequently more routes should be passed. Fig. 11(b) illustrates the increasing behavior of OFV2 as tvc increases. Based on the indicated current state of SBDC, it is revealed that the same delivery amount can be done with lower capacities (or lesser amount of vehicles) since for tvc ⩾ 80 , we see a constant value for both objective functions. Consequently, SBDC can significantly reduce the vehicles’ amortization/service cost. As mentioned before, about 30% of donations to SBDC are occurred directly to the MC. Now, let us call this ratio δ . Fig. 12 shows the sensitivity analysis on different values of δ . As can be seen, with increasing this value, the OFV1 decreases, since transporting products from BMs to SBDC is shortened and SBDC can fulfill its demand by direct donations. Now, if SBDC tries to increase the ratio by 10%, it can save approximately 5.26% in its total cost. In other words, the maximum cost

Fig. 12. OFV1 vs. alteration of δ.

sets are conducted on both model’s critical parameters and the parameters of the solution approach. The behaviors of the first and second objective functions under changes in capacity levels of BMs are demonstrated in Fig. 10(a) and (b), respectively. As can be seen, with increasing the capacity levels, the first objective function is decreased. However, for the values Cap (i) ⩽ 100 , we see a significant shift in the OFV. In this interval, based on the capacity limitations, the model is faced with a significant unmet demand, whose penalty cost increases the total cost of the network. For the second objective function, as the Cap (i) increases, the OFV2 increases as well until the point Cap (i) = 1000 ; thereafter, it 12

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

Fig. 13. OFV1 (a) OFV2 (b) and vs. alteration of u.

further research. Furthermore, in terms of solution methodology, introducing other uncertainty programming methods such as mixed fuzzy stochastic programming and comparing the efficiency of the results with the current approach is recommended. Applying other meta-heuristic solution techniques may also elevate the accuracy of the results.

that SBCD can be imposed under advertisement or other activities of the sort to raise the ratio by 10% is 194,477 (i.e. the difference between the current total cost under δ = 0.3 and the cost with δ = 0.4 ). It should be noted that alteration of δ does not affect the OFV2. Another parameter, which significantly affects the network’s storage policy and consequently OFVs is the shelf life of products (in our case, PLT). Therefore, a sensitivity analysis on this parameter can provide valuable managerial implications. As stated before, to make our network reliable, the considered shelf life for platelet is 5 days (this life time varies between 5 and 7 days, 5 ⩽ u ⩽ 7 according to the environmental situation). In the next sensitivity analyses, the behavior of the OFVs under different values of u (3 ⩽ u ⩽ 10 ) is tracked. As Fig. 13(a) shows, with increasing the value of u, the first objective function decreases, as the feasible region for optimality increases. The same behavior (i.e. approaching to the optimality as u increases) is seen for OFV2 in Fig. 13(b). It should be noted that the saving cost of the network in the case of increasing the life time of PLT for 1 day (i.e. u = 6) is about 0.015% of the current total cost. In other words, this amount is the maximum investment cost in technology and laboratory environment in order to increase the life time of PLT for 1 day.

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cie.2018.05.041. References Abdulwahab, U., & Wahab, M. (2014). Approximate dynamic programming modeling for a typical blood platelet bank. Computers & Industrial Engineering, 78, 259–270. Behnamian, J., Zandieh, M., & Ghomi, S. F. (2009). Parallel-machine scheduling problems with sequence-dependent setup times using an ACO, SA and VNS hybrid algorithm. Expert Systems with Applications, 36, 9637–9644. Beliën, J., & Forcé, H. (2012). Supply chain management of blood products: A literature review. European Journal of Operational Research, 217, 1–16. Birge, J. R., & Louveaux, F. (2011). Introduction to stochastic programming. Springer Science & Business Media. Blake, J. T. (2009). On the use of Operational Research for managing platelet inventory and ordering. Transfusion, 49, 396–401. Cetin, E., & Sarul, L. S. (2009). A blood bank location model: A multi-objective approach. European Journal of Pure and Applied Mathematics, 2, 112–124. Chapman, J. (2007). Unlocking the essentials of effective blood inventory management. Transfusion, 47, 190S–196S. Črepinšek, M., Liu, S.-H., & Mernik, L. (2012). A note on teaching–learning-based optimization algorithm. Information Sciences, 212, 79–93. Deb, K. (2001). Multi-objective optimization using evolutionary algorithms, Vol. 16. John Wiley & Sons. Drackley, A., Newbold, K. B., Paez, A., & Heddle, N. (2012). Forecasting Ontario's blood supply and demand. Transfusion, 52, 366–374. Eyles, N. H. J., Webert, K., Arnold, E., & McCurdy, B. (2011). Do expert assessments converge? An exploratory case study of evaluating and managing a blood supply risk. BMC Public Health, 11. Ghandforoush, P., & Sen, T. K. (2010). A DSS to manage platelet production supply chain for regional blood centers. Decision Support Systems, 50, 32–42. Gunpinar, S., & Centeno, G. (2015). Stochastic integer programming models for reducing wastages and shortages of blood products at hospitals. Computers & Operations Research, 54, 129–141. Haijema, R., van der Wal, J., & van Dijk, N. M. (2007). Blood platelet production: Optimization by dynamic programming and simulation. Computers & Operations Research, 34, 760–779. Haijema, R., van Dijk, N., van der Wal, J., & Smit Sibinga, C. (2009). Blood platelet production with breaks: Optimization by SDP and simulation. International Journal of Production Economics, 121, 464–473. Hansen, P., & Mladenović, N. (2001). Variable neighborhood search: Principles and applications. European Journal of Operational Research, 130, 449–467. Heddle, N. M., Liu, Y., Barty, R., Webert, K. E., Whittaker, S., Gagliardi, K., et al. (2009). Factors affecting the frequency of red blood cell outdates: An approach to establish benchmarking targets. Transfusion, 49, 219–226. Hemmelmayr, V., Doerner, K. F., Hartl, R. F., & Savelsbergh, M. W. P. (2010). Vendor managed inventory for environments with stochastic product usage. European Journal of Operational Research, 202, 686–695. Hussein, E., & Teruya, J. (2012). Evaluation of blood supply operation and infectious disease markers in blood donors during the Egyptian revolution. Transfusion, 52, 2321–2328.

7. Conclusion This paper develops a bi-objective mixed integer programming model for an integrated planning of main processes related to blood products. The proposed model addresses conflicting measures (i.e. reliability, efficiency, and efficacy) in a blood main center. In order to cope with the uncertain nature of the supply and demand, a multi-stage stochastic programming approach with a combined scenario tree is utilized. A main limitation to the presented model and the applied multi-stage stochastic programming approach is the high complexity of the model which requires developing heuristic/metaheuristic algorithms to efficiently solve the model at large scale instances. In doing so, a novel metaheuristic algorithm is also developed whose results show its superiority compared to the other competing algorithms. Although this paper is focused on the blood supply chain network design, the general framework presented here with some minor modifications can be applied to many other healthcare systems/supply chains with the same structure more especially those industries with major concerns regarding perishability and high degree of uncertainty in supply and demand data. There are several research directions that can be pursued based on the current work. The inventory decisions can be extended throughout the supply chain by which the limited inventory is also forced at demand points, along with the main centers. In addition, considering strategic decisions (e.g. facility location of MCs) and selecting reliable facility locations by taking into account disruption probabilities and/or backorder shortage in coping with demands are other worthwhile directions for 13

Computers & Industrial Engineering 122 (2018) 1–14

B. Zahiri et al.

Rytilä, J. S., & Spens, K. M. (2006). Using simulation to increase efficiency in blood supply chains. Management Research News, 29, 801–819. Şahin, G., Süral, H., & Meral, S. (2007). Locational analysis for regionalization of Turkish Red Crescent blood services. Computers & Operations Research, 34, 692–704. Şahinyazan, F. G., Kara, B. Y., & Taner, M. R. (2015). Selective vehicle routing for a mobile blood donation system. European Journal of Operational Research, 245, 22–34. Seifried, E., Klueter, H., Weidmann, C., Staudenmaier, T., Schrezenmeier, H., Henschler, R., et al. (2011). How much blood is needed? Vox Sanguinis, 100, 10–21. Sha, Y., & Huang, J. (2012). The multi-period location-allocation problem of engineering emergency blood supply systems. Systems Engineering Procedia, 5, 21–28. Stanger, S. H. W., Yates, N., Wilding, R., & Cotton, S. (2012). Blood inventory management: Hospital best practice. Transfusion Medicine Reviews, 26, 153–163. The World Bank (2015). < http://data.worldbank.org/indicator/SH.XPD.TOTL.ZS > . Accessed 01.05.16. Torabi, S. A., & Moghaddam, M. (2012). Multi-site integrated production-distribution planning with trans-shipment: A fuzzy goal programming approach. International Journal of Production Research, 50, 1726–1748. Van Der Bom, J. G., Heckbert, S. R., Lumley, T., Holmes, C. E., Cushman, M., Folsom, A. R., et al. (2009). Platelet count and the risk for thrombosis and death in the elderly. Journal of Thrombosis and Haemostasis, 7, 399–405. Van Dijk, N., Haijema, R., Van Der Wal, J., & Sibinga, C. S. (2009). Blood platelet production: A novel approach for practical optimization. Transfusion, 49, 411–420. Vera Hemmelmayr, K. F. D., Hartl, R. F., & Savelsbergh, M. W. P. (2009). Delivery strategies for blood products supplies. OR Spectrum, 31, 707–725. Zahiri, B., & Pishvaee, M. S. (2017). Blood supply chain network design considering blood group compatibility under uncertainty. International Journal of Production Research, 1–21. Zahiri, B., Tavakkoli-Moghaddam, R., Mohammadi, M., & Jula, P. (2014). Multi-objective design of an organ transplant network under uncertainty. Transportation Research Part E: Logistics and Transportation Review, 72, 101–124. Zahiri, B., Torabi, S. A., Mousazadeh, M., & Mansouri, S. (2015). Blood collection management: Methodology and application. Applied Mathematical Modelling, 39, 7680–7696. Zahiri, B., Torabi, S. A., & Tavakkoli-Moghaddam, R. (2017). A novel multi-stage possibilistic stochastic programming approach (with an application in relief distribution planning). Information Sciences, 385–386, 225–249. Zhou, D., Leung, L. C., & Pierskalla, W. P. (2011). Inventory management of platelets in hospitals: Optimal inventory policy for perishable products with regular and optional expedited replenishments. Manufacturing & Service Operations Management, 13, 420–438.

Jabbarzadeh, A., Fahimnia, B., & Seuring, S. (2014). Dynamic supply chain network design for the supply of blood in disasters: A robust model with real world application. Transportation Research Part E: Logistics and Transportation Review, 70, 225–244. Kaiser Health News (2015). < http://www.kaiserhealthnews.org/Daily-Reports/2012/ June/13/health-care-costs.aspx > . Accessed 01.05.16. Kallrath, J. (2005). Solving planning and design problems in the process industry using mixed integer and global optimization. Annals of Operations Research, 140, 339–373. Kamp, C., Heiden, M., Henseler, O., & Seitz, R. (2010). Management of blood supplies during an influenza pandemic. Transfusion, 50, 231–239. Katsaliaki, K. (2008). Cost-effective practices in the blood service sector. Health Policy, 86, 276–287. Kundu, D., Suresh, K., Ghosh, S., Das, S., Panigrahi, B. K., & Das, S. (2011). Multi-objective optimization with artificial weed colonies. Information Sciences, 181, 2441–2454. Madden, E., Murphy, E. L., & Custer, B. (2007). Modeling red cell procurement with both double-red-cell and whole-blood collection and the impact of European travel deferral on units available for transfusion. Transfusion, 47, 2025–2037. Margaret, F. S., Brandeau, L., & Pierskalla, W. P. (2004). Operations research and health care a handbook of methods and applications. Mladenović, N., & Hansen, P. (1997). Variable neighborhood search. Computers & Operations Research, 24, 1097–1100. Mohammadi, M., Jolai, F., & Tavakkoli-Moghaddam, R. (2013). Solving a new stochastic multi-mode p-hub covering location problem considering risk by a novel multi-objective algorithm. Applied Mathematical Modelling, 37, 10053–10073. Mousazadeh, M., Torabi, S. A., & Zahiri, B. (2015). A robust possibilistic programming approach for pharmaceutical supply chain network design. Computers & Chemical Engineering. Mustafee, N., Taylor, S. J., Katsaliaki, K., & Brailsford, S. (2009). Facilitating the analysis of a UK national blood service supply chain using distributed simulation. Simulation, 85, 113–128. Nagurney, A., Masoumi, A., & Yu, M. (2012). Supply chain network operations management of a blood banking system with cost and risk minimization. Computational Management Science, 9, 205–231. Osorio, A. F., Brailsford, S. C., & Smith, H. K. (2015). A structured review of quantitative models in the blood supply chain: A taxonomic framework for decision-making. International Journal of Production Research, 1–22. Perera, G., Hyam, C., Taylor, C., & Chapman, J. F. (2009). Hospital Blood Inventory Practice: The factors affecting stock level and wastage. Transfusion Medicine, 19, 99–104. Prastacos, G. P., & Brodheim, E. (1980). PBDS: A decision support system for regional blood management. Management Science, 26, 451–463.

14