Blood inventory-routing problem under uncertainty

Abstract

The problem of delivering blood products from community blood centers to the demand points including hospital blood banks falls within the context of perishable inventory-routing problems (PIRP). This is due to the fact that the delivery should be made on the right time with the right delivery quantity at the right place such that the total possible perished items as well as routing and inventory costs are minimized. However, some unique characteristics of blood logistics including assigned and unassigned inventories, crossmatch release period, transfusion to crossmatch ratio and older-first policy have made the problem more difficult than the routine PIRPs and thus proposing a new modeling of the problem is required. In this paper, we first develop a mixed integer programming formulation for blood inventory-routing problem. Then, to cope with uncertainties, a novel robust possibilistic programming (RPP) approach is proposed. Afterward, a novel iterative branch-and-cut is developed to solve a number of numerical examples to optimality. Finally, by implementing a test scenario on the data inspired from a real Iranian blood supply chain, the significance and applicability of the proposed model and RPP approach is proven.

Keywords

Healthcare management blood inventory-routing perishable products robust possibilistic programming branch-and-cut algorithm

1 Introduction

Inventory-routing problem (IRP) simultaneously determines when, how, and how much products should be delivered to the customers by integrating inventory and routing decisions. The work of Bell et al. [1] is the first research conducted on IRP, 33 years ago. Since then the concept has been evolved and many of its variants are introduced and applied to many real applications such as bottled water supply chain [2], crude oil transportation network [3], fuel delivery network [4] and waste vegetable oil collection problem for biodiesel production [5]. However, the application of IRP to blood distribution networks has never been tackled.

Blood inventory-routing problem (Blood-IRP) which is introduced in this paper has several specific features. First of all, blood products are perishable [6]. It means they become spoiled after a known period of time. Although this is also the case for radioactive and chemical products [7] and food industries such as dairy products, vegetables and fruits [8, 9], the significance of considering perishability of blood products is magnified due to its fatal consequences.

Recently, some authors have developed IRP models to take into account the perishability of products [10 –13] and called this class as PIRP. But these advancements in the context of IRP is not yet enough to address the other major challenges of Blood-IRPs.

Secondly, blood can be produced only by human beings and no other product can be used as an alternative [14]. In Iran the whole blood units are collected through a donation-based process which implies how uncertain the amount of collected units is. Community blood centers (CBC) are responsible for collecting whole blood units, decomposing them into their components and delivering the components to the demand points including hospital blood banks (HBB) and other medical centers and laboratories. Moreover, when these products are received at HBBs, a number of units are crossmatched with the patients’ bloods and kept reserved somewhere called assigned inventory. The crossmatched units remain in the assigned inventory for a maximum period of time called crossmatch release period (CRP) until they are either transfused to patients or turned back to the unassigned inventory for future utilizations. The ratio of number of transfused units to crossmatched units is called transfusion to crossmatch ratio [15]. In this process, the units that become spoiled are considered as wastage and are discarded.

The above-mentioned features are unique to blood logistics that are not yet tackled in the current literature. Thus, proposing a new IRP model capable of handling these specific challenges for blood supply chains is imperative. Furthermore, because of the indefinite emergency demands of blood units, employing a robust or fuzzy approach to cope with Blood-IRP is necessary.

In this paper, a mixed integer linear programming (MILP) formulation for modeling Blood-IRP under uncertainty is developed for the first time. All of the distinctive characteristics of the blood logistics are considered in the proposed formulation. In this model, older-first policy (OF) is adopted as the preferred issuing policy in HBBs and order-up-to level (OU) policy is adopted as the dominant policy of the blood inventory management. To cope with the uncertainty of parameters, a novel robust possibilistic programming approach is utilized. Moreover, a novel iterative Branch-and-cut (B&C) is proposed to solve a number of numerical examples to optimality. Finally, by conducting sensitivity analysis on the data inspired from a real Iranian blood supply chain the implication and applicability of the model and the solution approach is proven.

The rest of the paper is organized as follows. In Section 2 Blood-IRP is officially described and its notations are provided. In Section 3 the proposed mathematical formulation along with its underlying assumptions is elaborated. The robust possibilistic programming approach is explained in Section 4 and the robust possibilistic version of the proposed formulation is presented. Section 5 describes the developed B&C method. In Section 6 the data inspired from an Iranian blood supply chain are illustrated and numerical results and discussions are provided. Finally, in Section 7 conclusion remarks are drawn and research directions for future studiesare given.

2 Problem description

The description of the problem is formally provided here.

Consider G = (V, E) a graph,

Where V ={ 0, …, n } is the vertex set and E ={ (i, j) : i, j ∈ V, i < j } is the arc set while a routing cost rc_ij is associated with each arc (i, j) ∈ E. Vertex 0 stands for the CBC and the other vertices V′ = V ∖ {0} are symbols of HBBs. Consider set T = {1, …, lp} as the planning horizon of length lp and set K = {1, …, nv} as the set of vehicles available at CBC at each time period. Each vehicle has a limited capacity of Qu to be respected. Also, to each unit of blood product an age group g is assigned where g ∈ SL = {0, … sl} and sl is the shelf life of blood units. When one unit of blood becomes outdated a fuzzy wastage cost of $\tilde{oc}$ is incurred.

At the beginning of the planning horizon the initial inventories of HBBs are known to CBC. CBC is provided with $\tilde{r^{t}}$ quantity of fresh blood units at the start of each period. Demand of each HBB at each period $\tilde{{de}_{i}^{t}}$ is an uncertain parameter where its range is determined by experts and estimations regarding both elective and emergency demands. When HBB i is visited by vehicle k in period t, a quantity of $\sum_{g \in SL} q_{i}^{gkt}$ blood units is delivered such that its inventory level reaches a predefined target inventory level TIL_i. At each period t, $\tilde{{de}_{i}^{t}}$ amount of blood units are crossmatched and kept in assigned inventory. A ratio p of these crossmatched units are transfused to patients and the remaining units that are not yet outdated will be turned back to the unassigned inventory after crossmatch release period R for future usage. For each unit of age g remains in inventory at the end of each period, $\tilde{{hc}_{i}^{g}}$ Unit holding cost isincurred.

The problem is to decide on the time to visit each HBB, the route of the visited HBBs and the quantity to be delivered to each of them such that the mentioned restrictions are met while the inventory, transportation and wastage costs are minimized.

In the following the notation used in the proposed MILP formulation is presented.

Indices:

i, j∈ V = { 0, …, n }

Indices of locations where 0 represents CBC and V′ = V∖ { 0 } is set of HBBs

g∈ SL = { 0, …, sl }

Index of age groups that belongs to the discrete set SL (shelf life) where 0 represents fresh blood units

t∈ T = { 1, …, lp }

Index of time period where lp is the length of planning horizon

k∈ K = { 1, …, nv }

Index of vehicles. A set K of vehicles are available.

Parameters:

Total number of HBBs

Shelf life of blood units after which the blood units are outdated and become unsuitable for transfusion

The length of planning horizon

Total number of available vehicles

\tilde{{hc}_{i}^{g}}

Unit holding cost of blood of age g for HBB i

\tilde{oc}

Unit wastage cost of blood

rc _ij

Routing cost associated with edge (i, j) ∈ E where E ={ (i, j) : i, j ∈ V, i < j }

TIL _i

Target inventory level determined for each HBB

\tilde{r^{t}}

Number of fresh blood units CBC receives at each time period

\tilde{{de}_{i}^{t}}

Demand of each HBB i for each time period t

Vehicle capacity (a homogeneous fleet of vehicles are considered)

Transfusion to crossmatch ratio

Crossmatch release period

Variables:

${Inv}_{i}^{gt}$

Inventory level of blood units of age g at the end of time t at each HBB

u_{i}^{t}

Number of outdated blood units at each HBB at the end of period t

x_{ij}^{kt}

Number of times edge (i, j) is traversed by vehicle k in period t

q_{i}^{gkt}

The quantity of blood units of age g delivered by vehicle k in period t to HBB i

y_{i}^{kt}

Binary variable equal to one if and only if HBB i is visited by vehicle k in period t

z_{i}^{t}

Binary variable equal to one if and only if HBB i is visited in period t

β_{i}^{gt}

Number of blood units of age g returned from assigned inventory of HBB i to its unassigned inventory at the beginning of time t

d_{i}^{gt}

Number of blood units of age g used to satisfy the demand of HBB i in period t

L_{i}^{gt}

Binary variable equal to one if and only if blood units of age g are used to satisfy demand in time period t

α_{i}^{gt}

An auxiliary variable associated with age group g in time t that captures the number of blood units in an age group left to be utilized for the subsequent period if all available blood in this age group is not completely used to meet the demand of present period

To provide a better understanding of the problem, Fig. 1 shows the inventory flow of product of age g at HBB i in period t.

3 Mathematical formulation

The proposed mathematical formulation is based on the assumptions provided below.

The developed Blood-IRP model belongs to PIRP class of problems. It means that blood units are perishable and have a defined shelf life.

This model has uncertain parameters and hence is not deterministic. Demand of each HBB at each period, the quantity of fresh blood units made available at CBC at each time period, unit holding cost of blood at each HBB and unit wastage cost of blood are the ambiguous parameters of the model. To cope with these imprecise parameters, trapezoidal possibility distributions are adopted.

Stock-out is not permitted.

The OU policy is considered as the dominant inventory management policy at HBBs.

The fleet of vehicles is homogeneous.

The time horizon is finite.

The issuing policy at both levels of the supply chain is set to be Older-first policy. It means that while a product of older age is available, the use of the fresher units is not permitted.

Assigned and unassigned inventories, crossmatch release period and transfusion to crossmatch ratio are considered in the developed formulation.

Partial delivery is not acceptable.

3.1 Mixed integer programming formulation

The objective function of the proposed formulation is as follows. $\begin{array}{l} \min \sum_{i \in V} \sum_{g \in S L} \sum_{t \in T} \tilde{{hc}_{i}^{g}} . I n v_{i}^{g t} + \sum_{i \in V^{'}} \sum_{t \in T} \tilde{o c} . u_{i}^{t} \\ + \sum_{i \in V} \sum_{\begin{matrix} j \in V \\ i < j \end{matrix}} \sum_{k \in K} \sum_{t \in T} r c_{i j} . x_{i j}^{k t} \end{array}$ (1)

Subjected to the following constraints: $\begin{matrix} q_{i}^{gkt} \leq {TIL}_{i} . y_{i}^{kt} & \begin{matrix} i \in V^{'}, g \in SL \\ k \in K, t \in T \end{matrix} \end{matrix}$ (2) $\begin{matrix} \begin{matrix} \sum_{g \in SL} \sum_{k \in K} q_{i}^{gkt} + \sum_{g \in SL} {Inv}_{i}^{g - 1, t - 1} + \\ \sum_{g \in SL} β_{i}^{gt} \leq {TIL}_{i} \end{matrix} & i \in V^{'}, t \in T \end{matrix}$ (3) $\begin{matrix} \begin{matrix} \sum_{g \in SL} \sum_{k \in K} q_{i}^{gkt} \geq {TIL}_{i} . z_{i}^{t} - \\ \sum_{g \in SL} {Inv}_{i}^{g - 1, t - 1} - \sum_{g \in SL} β_{i}^{gt} \end{matrix} & i \in V^{'}, t \in T \end{matrix}$ (4) ${Inv}_{0}^{gt} = {Inv}_{0}^{g - 1, t - 1} - \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{gkt} g \in SL ∖ {0}, t \in T$ (5) ${Inv}_{0}^{0 t} = \tilde{r^{t}} - \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{0 kt} t \in T$ (6) $\begin{matrix} \begin{matrix} {Inv}_{i}^{gt} = {Inv}_{i}^{g - 1, t - 1} + \sum_{k \in K} q_{i}^{gkt} \\ + β_{i}^{gt} - d_{i}^{gt} \end{matrix} & \begin{matrix} i \in V^{'}, g \in SL ∖ {0}, \\ t \in T \end{matrix} \end{matrix}$ (7) ${Inv}_{i}^{0 t} = \sum_{k \in K} q_{i}^{0 kt} - d_{i}^{0 t} i \in V^{'}, t \in T$ (8) $\begin{matrix} d_{i}^{gt} \leq {TIL}_{i} . L_{i}^{gt} & \begin{matrix} i \in V^{'}, g \in SL, \\ t \in T \end{matrix} \end{matrix}$ (9) $\begin{matrix} L_{i}^{g - 1, t} \leq L_{i}^{gt} & \begin{matrix} i \in V^{'}, t \in T \\ g \in SL ∖ {0} \end{matrix} \end{matrix}$ (10) $\begin{matrix} \begin{matrix} {TIL}_{i} . (1 - L_{i}^{g - 1, t}) \geq \sum_{j = g}^{sl} {Inv}_{i}^{j - 1, t - 1} + \\ \sum_{j = g}^{sl} \sum_{k \in K} q_{i}^{jkt} + \sum_{j = g}^{sl} β_{i}^{jt} - \tilde{{de}_{i}^{t}} + 1 \end{matrix} & \begin{matrix} i \in V^{'}, t \in T \\ g \in SL ∖ {0} \end{matrix} \end{matrix}$ (11) $\sum_{g \in SL} {Inv}_{i}^{gt} \leq {TIL}_{i} i \in V^{'}, t \in T$ (12) $\tilde{{de}_{i}^{t}} = \sum_{g \in SL} d_{i}^{gt} i \in V^{'}, t \in T$ (13) $\sum_{i \in V^{'}} \sum_{g \in SL} q_{i}^{gkt} \tilde{\leq} Qu . y_{0}^{kt} k \in K, t \in T$ (14) $\begin{matrix} α_{i}^{gt} \leq (L_{i}^{gt} - L_{i}^{g - 1, t}) . i \in V^{'}, t \in T \\ ({Inv}_{i}^{g - 1, t - 1} + \sum_{k \in K} q_{i}^{gkt} + β_{i}^{gt}) g \in SL ∖ {0} \end{matrix}$ (15) $\begin{matrix} \begin{matrix} \begin{matrix} {Inv}_{i}^{gt} = (L_{i}^{gt} - L_{i}^{g - 1, t}) . α_{i}^{gt} + \\ (1 - L_{i}^{gt}) . ({Inv}_{i}^{g - 1, t - 1} + \sum_{k \in K} q_{i}^{gkt} + β_{i}^{gt}) \end{matrix} & \begin{matrix} i \in V^{'}, t \in T \\ g \in SL ∖ {0} \end{matrix} \end{matrix} \end{matrix}$ (16)

$\begin{matrix} \begin{matrix} β_{i}^{gt} = ⌊ \begin{matrix} [\begin{matrix} (\begin{matrix} {Inv}_{i}^{g - R - 1, t - R - 1} + \\ \sum_{k \in K} q_{i}^{g - R, k, t - R} + β_{i}^{g - R, t - R} \end{matrix}) . \\ L_{i}^{g - R, t - R} - α_{i}^{g - R, t - R} \end{matrix}] . \\ (1 - p) \end{matrix} ⌋ \end{matrix} & \begin{matrix} g = 1 + R, \dots, sl + R, \\ t = R + 1, \dots, T, i \in V^{'} \end{matrix} \end{matrix}$ (17) $\begin{matrix} \sum_{j \in V, i < j} x_{ij}^{kt} + \sum_{j \in V, j < i} x_{ji}^{kt} = 2 . y_{i}^{kt}, & \begin{matrix} i \in V, k \in K \\ t \in T \end{matrix} \end{matrix}$ (18) $\begin{matrix} \sum_{i \in S} \sum_{j \in S, i < j} x_{ij}^{kt} \leq \sum_{i \in S} y_{i}^{kt} - y_{m}^{kt} & \begin{matrix} S \subseteq V^{'}, k \in K, \\ t \in T, m \in S \end{matrix} \end{matrix}$ (19) $\begin{matrix} \sum_{k \in K} y_{i}^{kt} \leq 1 & i \in V^{'}, t \in T \end{matrix}$ (20) $\begin{matrix} {Inv}_{i}^{gt} = β_{i}^{gt} & \begin{matrix} g = sl + 1, \dots, sl + R, \\ i \in V^{'}, t \in T \end{matrix} \end{matrix}$ (21) $\begin{matrix} β_{i}^{gt} = 0, & \begin{matrix} g = 1 + R, \dots, sl + R \\ t = 1, \dots R, i \in V^{'} \end{matrix} \end{matrix}$ (22) $\begin{matrix} β_{i}^{gt} = 0 & \begin{matrix} g = 0, \dots, R, \\ i \in V^{'}, t \in T \end{matrix} \end{matrix}$ (23) $\begin{matrix} u_{i}^{t} = \sum_{j = 0}^{R} {Inv}_{i}^{sl + j, t} & i \in V^{'}, t \in T \end{matrix}$ (24) $\begin{matrix} z_{i}^{t} \leq \sum_{k \in K} y_{i}^{kt} & i \in V^{'}, t \in T \end{matrix}$ (25) $\begin{matrix} z_{i}^{t} \geq \frac{1}{nv} . \sum_{k \in K} y_{i}^{kt} & i \in V^{'}, t \in T \end{matrix}$ (26) $\begin{matrix} \begin{matrix} {Inv}_{i}^{gt}, d_{i}^{gt}, q_{i}^{gkt}, α_{i}^{gt}, β_{i}^{gt}, \\ u_{i}^{t} \in ℤ^{+} \end{matrix} & \begin{matrix} i \in V^{'}, g \in SL, \\ t \in T, k \in K \end{matrix} \end{matrix}$ (27) $\begin{matrix} x_{0 i}^{kt} \in {0, 1, 2}, & \begin{matrix} i \in V^{'}, k \in K \\ t \in T \end{matrix} \end{matrix}$ (28) $\begin{matrix} x_{ij}^{kt} \in {0, 1}, & \begin{matrix} i, j \in V^{'}, k \in K \\ t \in T \end{matrix} \end{matrix}$ (29) $\begin{matrix} L_{i}^{gt} \in {0, 1}, & \begin{matrix} i \in V^{'}, g \in SL \\ t \in T \end{matrix} \end{matrix}$ (30) $\begin{matrix} y_{i}^{kt} \in {0, 1}, & \begin{matrix} i \in V, k \in K \\ t \in T \end{matrix} \end{matrix}$ (31)

Total costs of the supply chain are minimized using objective function (1). In Equation (1), routing costs are deterministic because they are computed by the product of three deterministic parameters; The length of route (in kilometers), the gasoline used (liter per kilometer) and the price of each liter of gasoline (Currency per liter). The length of the route between each pair of demand points is a fixed value. The gasoline used per kilometer is also a vehicle specific feature and is considered as a fixed value. Finally, in Iran, the price of each liter of gasoline is a constant value that is determined by the government annually and hence can be regarded as a deterministic parameter. Thus, the product of these three constant values (i.e., routing costs) is also deterministic in nature. The Order-up-to level inventory policy is imposed by constraints (2) – (4). Constraint (5) indicates the inventory balance of CBC for each period and each age group and constraint (6) denotes CBC is provided with $\tilde{r^{t}}$ quantity of fresh blood units at the start of each period. Since the supply of fresh blood units from collection points is uncertain due to its donation-based nature, parameter $\tilde{r^{t}}$ becomes uncertain accordingly. The inventory balance of HBBs for non-fresh and fresh blood units is shown by constraints (7) and (8), respectively. Constraints (9) – (11) guarantee the OF issuing policy. Constraint (12) ensures that capacity constraints of HBBs are respected. Constraint (13) is obviously established by definition. Constraint (14) enforces vehicle capacity constraints. Constraints (15) – (17) together with (22) and (23) are incorporated in order to compute the quantity of blood units turned back from assigned inventory to unassigned inventory at the beginning of each period. Constraints (18) and (19) are degree constraints and subtour elimination constraints, respectively. Partial deliveries are avoided by constraint (20). Constraint (21) implies that HBBs should not use or keep the expired blood units. The number of wasted units is calculated using constraint (24). Inequalities (25) and (26) establish the relation between z and y variables. Integrality and non-negativity conditions on variables are imposed by constraints(27) – (31).

3.2 Linearization

The mixed integer programming formulation presented in the previous subsection includes nonlinear terms through inequalities (15) – (17) due to the interactions between binary and integer variables and the presence of floor function. Nevertheless, there exist well-known approaches to cope with such types of nonlinearities and convert them to their linear counterparts. After implementing the linearization techniques, a linear model will be developed that we refer to as mixed integer linear programming (MILP) model. To keep coherence and brevity, the details of the linearization techniques are not delineatedhere.

3.3 Valid inequalities

In order to solve the proposed model to optimality, using a branch-and-cut algorithm, it is necessary to strengthen the model through adding several valid inequalities: $\begin{matrix} x_{0 i}^{kt} \leq 2 . y_{i}^{kt}, & \begin{matrix} i \in V, k \in K \\ t \in T \end{matrix} \end{matrix}$ (32) $\begin{matrix} x_{ij}^{kt} \leq y_{i}^{kt} & \begin{matrix} i, j \in V^{'}, \\ k \in K, t \in T \end{matrix} \end{matrix}$ (33) $\begin{matrix} y_{i}^{kt} \leq y_{0}^{kt} & \begin{matrix} i \in V^{'}, k \in K, \\ t \in T \end{matrix} \end{matrix}$ (34) $\begin{matrix} y_{0}^{kt} \leq y_{0}^{k - 1, t} & \begin{matrix} k \in K ∖ {1}, \\ t \in T \end{matrix} \end{matrix}$ (35) $\begin{matrix} y_{i}^{kt} \leq \sum_{j < i} y_{j}^{k - 1, t} & \begin{matrix} k \in K ∖ {1}, \\ t \in T, i \in V \end{matrix} \end{matrix}$ (36)

Since the fleet of vehicles is homogeneous, vehicle symmetry breaking constraints (35) and (36) are valid for HBB vertices.

4 The proposed robust possibilistic programming approach

In some cases, as in blood supply chain, we encounter profound uncertainties where no information is available on the parameter probability distributions [16]. For example, as illustrated in the work of Pierskalla [15], daily O+ crossmatches for a large Chicago hospital ranges from a low of less than 10 units to a high of over 140 units and with considerable daily fluctuations throughout a year. An appropriate way of solving stochastic IRPs in such cases is the use of robust optimization. In this approach, the problem is optimized so that the problem remains feasible for all possible realizations of bounded parameters. For a comprehensive review on recent advances in robust optimization, the reader is referred to the work of Gabrel, Murat, and Thiele [17].

Pishvaee, Razmi and Torabi [18] extended the theory of robust optimization into the possibilistic programming framework and called their new approach as robust possibilistic programming (RPP). As is known, possibilistic programming copes with ambiguous coefficients of objective functions and constraints that are usually modeled considering available subjective and objective data [19]. Thus, RPP benefits from the capabilities of both robust optimization and possibilistic programming approaches and hence is a fruitful approach for our proposed uncertain mathematical formulation.

However, RPP proposed by Pishvaee, Razmi and Torabi [18] does not take into account equality constraints with necessity measures (constraints (6) and (13)) and fuzzy inequalities (constraint (14)). Therefore, in this study, a novel RPP is developed which is able to address the mentioned challenges. The new approach is as follows.

4.1 The proposed RPP approach

Step 1: Convert all the fuzzy inequalities into their crisp counterparts using Verdegey’s method [20]. The crisp version of constraint (14) is presented below: $\sum_{i \in V^{'}} \sum_{g \in SL} q_{i}^{gkt} \leq (Qu + (1 - ω) \tilde{q}) . y_{0}^{kt} k \in K, t \in T$ (37)

Where ω is a real number between 0 and 1 denoting ω-cut concept in fuzzy mathematical programming and $\tilde{q}$ is fuzzy tolerance (permitted violation) of the corresponding constraint.

Step 2: Obtain the possibilistic chance constrained programming (PCCP) version of the proposed model. For equality constraints with necessity measures (here constraints (6) and (13)) use the method of Zahiri, Tavakkoli-Moghaddam and Pishvaee [21]. The PCCP is as follows: $\begin{array}{l} \min \sum_{i \in V} \sum_{g \in S L} \sum_{t \in T} E [\tilde{{hc}_{i}^{g}}] . I n v_{i}^{g t} \\ + \sum_{i \in V^{'}} \sum_{t \in T} E [\tilde{o c}] . u_{i}^{t} + \sum_{i \in V} \sum_{\begin{matrix} j \in V \\ i < j \end{matrix}} \sum_{k \in K} \sum_{t \in T} r c_{i j} . x_{i j}^{k t} \end{array}$ (38)

Subjected to the following constraints: $\begin{matrix} Nec {{Inv}_{0}^{0 t} = \tilde{r^{t}} - \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{0 kt}} \geq ω_{1} & t \in T \end{matrix}$ (39) $\begin{matrix} Nec {\begin{matrix} {TIL}_{i} . (1 - L_{i}^{g - 1, t}) \geq \sum_{j = g}^{sl} {Inv}_{i}^{j - 1, t - 1} + \\ \sum_{j = g}^{sl} \sum_{k \in K} q_{i}^{jkt} + \sum_{j = g}^{sl} β_{i}^{jt} - {\tilde{de}}_{i}^{t} + 1 \end{matrix}} \geq ω_{2} & \begin{matrix} i \in V^{'}, \\ t \in T, \\ g \in SL ∖ {0} \end{matrix} \end{matrix}$ (40) $\begin{matrix} Nec {{\tilde{de}}_{i}^{t} = \sum_{g \in SL} d_{i}^{gt}} \geq ω_{3} & \begin{matrix} i \in V^{'}, \\ t \in T \end{matrix} \end{matrix}$ (41) $\begin{matrix} Nec {\begin{matrix} \sum_{i \in V^{'}} \sum_{g \in SL} q_{i}^{gkt} \leq \\ [Qu + (1 - ω) \tilde{q}] . y_{0}^{kt} \end{matrix}} \geq ω_{4} & \begin{matrix} k \in K, \\ t \in T \end{matrix} \end{matrix}$ (42)

And constraints (2)–(5), (7)–(10), (12), linear counterparts of constraints (15)–(17), (18)–(36) where ω_i indicates the minimum satisfaction level of the corresponding constraint.

Step 3: Obtain the crisp counterpart of PCCP according to the work of Pishvaee, Razmi and Torabi [18] assuming trapezoidal possibility distribution for modeling imprecise parameters. We know that triangular possibility distribution is the most common tool for modeling the imprecise nature of the ambiguous parameters due to its computational efficiency and simplicity in data acquisition [22]. However, since triangular possibility distribution is a special case of trapezoidal possibility distribution, to make the approach more general, the latter is adopted in this paper.

Noteworthy, a trapezoidal possibility distribution is defined by its four prominent points, e.g., $\tilde{ψ} = (ψ_{(1)}, ψ_{(2)}, ψ_{(3)}, ψ_{(4)})$ . The crisp counterpart of PCCP (CCPCCP) is presented below. $\begin{array}{l} \min E [O F] = \sum_{i \in V} \sum_{g \in S L} \sum_{t \in T} (\frac{h c_{i (1)}^{g} + h c_{i (2)}^{g} + h c_{i (3)}^{g} + h c_{i (4)}^{g}}{4}) . I n v_{i}^{g t} \\ + \sum_{i \in V^{'}} \sum_{t \in T} (\frac{o c_{(1)} + o c_{(2)} + o c_{(3)} + o c_{(4)}}{4}) . u_{i}^{t} \\ + \sum_{i \in V} \sum_{\begin{matrix} j \in V \\ i < j \end{matrix}} \sum_{k \in K} \sum_{t \in T} r c_{i j} . x_{i j}^{k t} \end{array}$ (43)

Subjected to: $\begin{matrix} \begin{matrix} {Inv}_{0}^{0 t} \geq ω_{1} r_{(2)}^{t} + (1 - ω_{1}) r_{(1)}^{t} - \\ \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{0 kt} \end{matrix} & t \in T \end{matrix}$ (44) $\begin{matrix} \begin{matrix} {Inv}_{0}^{0 t} \leq ω_{1} r_{(3)}^{t} + (1 - ω_{1}) r_{(4)}^{t} - \\ \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{0 kt} \end{matrix} & t \in T \end{matrix}$ (45) $\begin{matrix} \begin{matrix} {TIL}_{i} . (1 - L_{i}^{g - 1, t}) \\ \geq \sum_{j = g}^{sl} {Inv}_{i}^{j - 1, t - 1} + \sum_{j = g}^{sl} \sum_{k \in K} q_{i}^{jkt} + \sum_{j = g}^{sl} β_{i}^{jt} \\ - ({de}_{i (2)}^{t} + ω_{2} ({de}_{i (1)}^{t} - {de}_{i (2)}^{t})) + 1 \end{matrix} & \begin{matrix} i \in V^{'}, t \in T \\ g \in SL ∖ {0} \end{matrix} \end{matrix}$ (46) $\begin{matrix} \sum_{g \in SL} d_{i}^{gt} \geq ω_{3} {de}_{i (2)}^{t} + (1 - ω_{3}) {de}_{i (1)}^{t} & i \in V^{'}, t \in T \end{matrix}$ (47) $\begin{matrix} \sum_{g \in SL} d_{i}^{gt} \leq ω_{3} {de}_{i (3)}^{t} + (1 - ω_{3}) {de}_{i (4)}^{t} & i \in V^{'}, t \in T \end{matrix}$ (48) $\begin{matrix} \begin{matrix} \sum_{i \in V^{'}} \sum_{g \in SL} q_{i}^{gkt} \\ \leq [\begin{matrix} Qu + \\ (1 - ω) . (ω_{4} q_{(1)} + (1 - ω_{4}) q_{(2)}) \end{matrix}] . y_{0}^{kt} \end{matrix} & k \in K, t \in T \end{matrix}$ (59)

And constraints (2)–(5), (7)–(10), (12), linear counterparts of constraints (15)–(17), (18)–(36) where OF is the abbreviation of objective function.

Remark. From decision maker’s (DM) viewpoint, two main problems confront the CCPCCP model provided above. Firstly, the model only minimizes the expected value of objective function and does not take deviations from this expected value into account. As a result, the DM can be exposed to high levels of risk in some real situations. Secondly, to find the best combination of values of ω _i , DM should set the values of ω _i using an exhaustive interactive sensitivity analysis approach. Indeed, there is no guarantee for the obtained values to be the best possible choice since the values are chosen in a somehow subjective manner. To cope with the mentioned problems, in Step 4, a robust optimization approach isapplied.

Step 4: Convert the model developed in Step 3 to the form of RPP-II model proposed by Pishvaee, Razmi and Torabi [18]. RPP-II not only minimizes the maximum deviation from the expected value but also determines the optimum ω_i values in a proactive manner. This is due to the fact that in RPP-II, ω_i are considered as variables and the model tries to find their optimum values. After conversion, remove nonlinear terms in the model, if needed.

The obtained RPP-II model is presented below. $\begin{matrix} min E [OF] + θ_{1} ({OF}_{max} - E [OF]) \\ + θ_{2} \sum_{t \in T} [\begin{matrix} r_{(2)}^{t} - ω_{1} r_{(2)}^{t} - (1 - ω_{1}) r_{(1)}^{t} \\ + ω_{1} r_{(3)}^{t} + (1 - ω_{1}) r_{(4)}^{t} - r_{(3)}^{t} \end{matrix}] \\ + θ_{3} \sum_{i \in V^{'}} \sum_{t \in T} [{de}_{i (2)}^{t} + ω_{2} ({de}_{i (1)}^{t} - {de}_{i (2)}^{t}) - {de}_{i (1)}^{t}] \\ + θ_{4} \sum_{i \in V^{'}} \sum_{t \in T} [\begin{matrix} {de}_{i (2)}^{t} - ω_{3} {de}_{i (2)}^{t} - (1 - ω_{3}) {de}_{i (1)}^{t} \\ + ω_{3} {de}_{i (3)}^{t} + (1 - ω_{3}) {de}_{i (4)}^{t} - {de}_{i (3)}^{t} \end{matrix}] \\ + θ_{5} \sum_{k \in K} \sum_{t \in T} [(1 - ω) . (ω_{4} q_{(1)} + (1 - ω_{4}) q_{(2)} - q_{(1)})] . y_{0}^{kt} \end{matrix}$ (50)

Subjected to the following constraints: $\begin{matrix} \begin{matrix} {Inv}_{0}^{0 t} \geq ω_{1} r_{(2)}^{t} + (1 - ω_{1}) r_{(1)}^{t} - \\ \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{0 kt} \end{matrix} & t \in T \end{matrix}$ (51) $\begin{matrix} \begin{matrix} {Inv}_{0}^{0 t} \leq ω_{1} r_{(3)}^{t} + (1 - ω_{1}) r_{(4)}^{t} - \\ \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{0 kt} \end{matrix} & t \in T \end{matrix}$ (52) $\begin{matrix} \begin{matrix} {TIL}_{i} . (1 - L_{i}^{g - 1, t}) \\ \geq \sum_{j = g}^{sl} {Inv}_{i}^{j - 1, t - 1} + \sum_{j = g}^{sl} \sum_{k \in K} q_{i}^{jkt} + \sum_{j = g}^{sl} β_{i}^{jt} \\ - ({de}_{i (2)}^{t} + ω_{2} ({de}_{i (1)}^{t} - {de}_{i (2)}^{t})) + 1 \end{matrix} & \begin{matrix} i \in V^{'}, t \in T \\ g \in SL ∖ {0} \end{matrix} \end{matrix}$ (53) $\begin{matrix} \sum_{g \in SL} d_{i}^{gt} \geq ω_{3} {de}_{i (2)}^{t} + (1 - ω_{3}) {de}_{i (1)}^{t} & i \in V^{'}, t \in T \end{matrix}$ (54) $\begin{matrix} \sum_{g \in SL} d_{i}^{gt} \leq ω_{3} {de}_{i (3)}^{t} + (1 - ω_{3}) {de}_{i (4)}^{t} & i \in V^{'}, t \in T \end{matrix}$ (55) $\begin{matrix} \begin{matrix} \sum_{i \in V^{'}} \sum_{g \in SL} q_{i}^{gkt} \\ \leq [\begin{matrix} Qu + \\ (1 - ω) . (ω_{4} q_{(1)} + (1 - ω_{4}) q_{(2)}) \end{matrix}] . y_{0}^{kt} \end{matrix} & k \in K, t \in T \end{matrix}$ (56) $\begin{matrix} 0.5 < ω_{i} \leq 1 & i \in {1, 2, 3, 4} \end{matrix}$ (57)

And constraints (2)–(5), (7)–(10), (12), linear counterparts of constraints (15)–(17), (18)–(36) where OF_max is computed using Equation (58) as follows: $\begin{array}{l} O F_{\max} \\ = \sum_{i \in V} \sum_{g \in S L} \sum_{t \in T} h c_{i (4)}^{g} . I n v_{i}^{g t} + \sum_{i \in V^{'}} \sum_{t \in T} o c_{(4)} . u_{i}^{t} \\ + \sum_{i \in V} \sum_{\begin{matrix} j \in V \\ i < j \end{matrix}} \sum_{k \in K} \sum_{t \in T} r c_{i j} . x_{i j}^{k t} \end{array}$ (58)

It should be pointed out that both the last term in objective function (50) and the right hand side of constraint (56) include nonlinear terms ( $ω_{4} . y_{0}^{kt}$ ). To eliminate the nonlinear term, an auxiliary variable $y_{0}^{kt}$ is defined such that: $\begin{matrix} y_{0}^{' kt} = ω_{4} . y_{0}^{kt} & k \in K, t \in T \end{matrix}$ (59)

Subjected to the following conditions: $\begin{matrix} y_{0}^{' kt} \leq M . y_{0}^{kt} & k \in K, t \in T \end{matrix}$ (60) $\begin{matrix} y_{0}^{' kt} \geq M . (y_{0}^{kt} - 1) + ω_{4} & k \in K, t \in T \end{matrix}$ (61) $\begin{matrix} y_{0}^{' kt} \leq ω_{4} & k \in K, t \in T \end{matrix}$ (62) $\begin{matrix} y_{0}^{' kt} \geq 0 & k \in K, t \in T \end{matrix}$ (63)

Where M is a predefined large number. Therefore, the final linearized RPP formulation (Final-RPP) is as follows: $\begin{matrix} min E [OF] + θ_{1} ({OF}_{max} - E [OF]) \\ + θ_{2} \sum_{t \in T} [\begin{matrix} r_{(2)}^{t} - ω_{1} r_{(2)}^{t} - (1 - ω_{1}) r_{(1)}^{t} \\ + ω_{1} r_{(3)}^{t} + (1 - ω_{1}) r_{(4)}^{t} - r_{(3)}^{t} \end{matrix}] \\ + θ_{3} \sum_{i \in V^{'}} \sum_{t \in T} [{de}_{i (2)}^{t} + ω_{2} ({de}_{i (1)}^{t} - {de}_{i (2)}^{t}) - {de}_{i (1)}^{t}] \\ + θ_{4} \sum_{i \in V^{'}} \sum_{t \in T} [\begin{matrix} {de}_{i (2)}^{t} - ω_{3} {de}_{i (2)}^{t} - (1 - ω_{3}) {de}_{i (1)}^{t} \\ + ω_{3} {de}_{i (3)}^{t} + (1 - ω_{3}) {de}_{i (4)}^{t} - {de}_{i (3)}^{t} \end{matrix}] \\ + θ_{5} \sum_{k \in K} \sum_{t \in T} [(1 - ω) . [q_{(1)} (y_{0}^{' kt} - y_{0}^{kt}) + q_{(2)} (y_{0}^{kt} - y_{0}^{' kt})]] \end{matrix}$ (64)

Subjected to: $\begin{matrix} \begin{matrix} \sum_{i \in V^{'}} \sum_{g \in SL} q_{i}^{gkt} \leq Qu . y_{0}^{kt} \\ + (1 - ω) . [q_{(1)} y_{0}^{' kt} + q_{(2)} (y_{0}^{kt} - y_{0}^{' kt})] \end{matrix} & k \in K, t \in T \end{matrix}$ (65)

And constraints (2)–(5), (7)–(10), (12), linear counterparts of constraints (15)–(17), (18)–(36), (51)–(55), (57), (58), (60)–(63).

In Final-RPP model, the first term of the objective function (64) minimizes the expected value of OF. The second term minimizes the maximum deviation of OF over its expected value while parameter θ₁ indicates the importance of this term compared to the other terms. In fact the second term controls the optimality robustness of the solution vector. The other terms in the objective function control the feasibility robustness of the solution vector by minimizing the confidence level of each constraint where θ_i (2 ≤ i ≤ 5) is the penalty unit of possible violation of each constraint that has imprecise parameters. Thus, in Final-RPP the values of parameters θ_i and ω are determined based on the DM’s preferences while the values of variables ω_i are determined by the model itself in a proactive manner.

5 Branch-and-cut algorithm

It is known that IRP is a NP-hard problem. Thus, our proposed model becomes intractable for large-scaled instances. But, for modest instance sizes it can be handled by efficient B&C algorithms [23 –26].

In this paper, a novel B&C is developed to solve the problem. B&Cs already proposed in the literature drop subtour elimination constraints (SEC) (19) from the formulation and during the solution procedure add them to any integer solution that violates them. However, in this approach after dropping constraints (19) we keep solving the formulation until reaching an optimum solution. Then, the violated constraints are identified and added to the problem and the problem is then reoptimized. The reason for employing this new approach is twofold. Firstly, a small number of cuts in the form of constraints (19) are generated in this approach because it only generates the active ones that can affect the optimum solution of the original problem. Secondly, this approach can be easily employed within the body of matheuristic algorithms that contain repairing procedures for eliminating violated subtour constraints. The three-step B&C algorithm is provided below.

5.1 The proposed Branch-and-cut algorithm

Step 1: Drop subtour elimination constraints (19) from the proposed mathematical formulation and solve the rest of the formulation using CPLEX software.

Step 2: Suppose z^* is the optimum solution found by the previous step and f (z^*) is its objective function value. Check if z^* has any violated subtour elimination constraints (VSEC).

If the answer is no,

z^* is also the optimum solution of the original problem and f (z^*) is the optimum objective function value and hence terminate the algorithm.

else, Use the subtour identification procedure (explained in subsection 5.2.) to find all of the connected graphs of each vehicle in each period. Then, add the corresponding SECs into the model and go to Step 3.

Step 3: Set z* as the initial infeasible solution of the new model. Solve the new model using CPLEX software with z* as its initial infeasible solution and then go to Step 2.

In order to make the best use of CPLEX capabilities in Steps 1 and 3 some remarks should be made. Firstly, it should be kept in mind that valid inequalities (32–36) should also be added to the model. These inequalities are very useful in helping CPLEX generate new cuts. Secondly, it is recommended that branching priorities be defined for the problem. It is clear by definition that a higher priority should be assigned to y variables rather than x variables.

In Step 3 when z* is entered as an initial solution to the CPLEX software, CPLEX starts several heuristics to repair the solution and find a feasible solution. Since z* is the optimal solution of the problem in the absence of SECs, the found feasible solution is a heuristic initial solution of good quality.

5.2 Subtour identification procedure

The primary role of Subtour identification procedure is to discover the subtours of an on hand solution and generate the corresponding VSECs.

In this procedure, a tour is initiated for each vehicle k in each time period t. The tour starts from the supplier and tries to meet all the nodes with $y_{i}^{kt} = 1$ only once and in a consecutive way and finishes at the supplier again. If such a tour exists it is concluded that the route of vehicle k in period t contains no subtours. If not, after reaching to the starting point some nodes will remain unmet. At this point, it is clear that a minimum of two SECs are violated. Then again a new tour is initiated from one arbitrary unmet node and navigates through its connected nodes until reaching to its starting point. This procedure is repeated until having all of the nodes visited. At this point, a number of VSECs equal to the number of subtours of vehicle k in period t are found. For all vehicles this procedure is repeated and all the VSECs are generated.

6 Implementation and evaluation

To evaluate the performance and applicability of the proposed CCPCCP and Final-RPP models, data derived from a real blood supply chain in Iran is used. In the considered supply chain, Sari (a big city in province of Mazandaran in northern Iran) is the supplier that provides 8 demand points including Galugah, Behshahr, Neka, Juybar, Babol, Qaemshahr, Zirab and Sari itself with blood units. Figure 2 shows the location of these cities on the map taken from Google Map free service.

The possibility distribution of unit holding costs is randomly generated. However, it has been estimated using the knowledge and data of field experts for the other imprecise parameters. The value of the salient parameters related to demand of each city at each period, number of fresh blood units becomes available at CBC (located in Sari) at each time period and distance between each pair of cities are illustrated through Tables 1 to 3, respectively. Besides crossmatch release period is set equal to 1 period, transfusion to crossmatch ratio is set to 0.5 and shelf life of blood units is considered equal to 2 periods. The values of other parameters are not provided here due to space limitations. However, they can be provided upon request.

In order to assess the robustness of solutions obtained by models CCPCCP and Final-RPP a test scenario is implemented. In this scenario, CCPCCP model for three different values of ω_i(0.7, 0.8 and 0.9) and Final-RPP model are firstly solved using nominal data provided in Tables 1–3. The outcomes of this phase are shown in Table 4. It is worth noting that all computations were done on a personal computer with two cores of 2.2 GHz and a 2 GB RAM and no parallel computing technology was accessible for the authors.

As can be seen from Table 4, the performance of Final-RPP under nominal data is worse than the other three models. Indeed, it is the cost that is incurred because of obtaining more robust solutions.

In the next step of the test scenario, to compare the desirability and robustness of solutions obtained in the previous step, 15 random realizations are generated uniformly and then the performance of the obtained solutions is tested under each realization. To do so, suppose that $\tilde{ψ} = (ψ_{(1)}, ψ_{(2)}, ψ_{(3)}, ψ_{(4)})$ is an imprecise parameter with trapezoidal possibility distribution function. The realization is generated by producing a random number uniformly from [ψ₍₁₎, ψ₍₄₎] interval. Subsequently, the solutions obtained by the models under nominal data (Inv^*, u^*, x^*, q^*, L^*, β^*, d^*, y^*) will be replaced in the following linear programming formulation: $\begin{array}{l} \min \sum_{i \in V} \sum_{g \in S L} \sum_{t \in T} h c_{r e a l i}^{g} . I n v_{i}^{* g t} + \sum_{i \in V^{'}} \sum_{t \in T} o c_{r e a l} . u_{i}^{* t} \\ + \sum_{i \in V} \sum_{\begin{matrix} j \in V \\ i < j \end{matrix}} \sum_{k \in K} \sum_{t \in T} r c_{i j} . x_{i j}^{* k t} \\ + θ_{2} . (\sum_{t \in T} S V_{1}^{t} + \sum_{t \in T} S V_{2}^{t}) \\ + θ_{3} . (\sum_{i \in V^{'}} \sum_{g \in S L \ {0}} \sum_{t \in T} S V_{3 i}^{g t}) \\ + θ_{4} . (\sum_{i \in V^{'}} \sum_{t \in T} S V_{4 i}^{t} + \sum_{i \in V^{'}} \sum_{t \in T} S V_{5 i}^{t}) \\ + θ_{5} . (\sum_{k \in K} \sum_{t \in T} S V_{6}^{k t}) \end{array}$ (66)

Subjected to the following constraints: ${Inv}_{0}^{* 0 t} + {SV}_{1}^{t} \geq r_{real}^{t} - \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{* 0 kt} t \in T$ (67) ${Inv}_{0}^{* 0 t} \leq r_{real}^{t} - \sum_{i \in V^{'}} \sum_{k \in K} q_{i}^{* 0 kt} + {SV}_{2}^{t} t \in T$ (68) $\begin{matrix} \begin{matrix} {TIL}_{i} . (1 - L_{i}^{* g - 1, t}) \geq \sum_{j = g}^{sl} {Inv}_{i}^{* j - 1, t - 1} \\ + \sum_{j = g}^{sl} \sum_{k \in K} q_{i}^{* jkt} + \sum_{j = g}^{sl} β_{i}^{* jt} - {de}_{real i}^{t} - \\ {SV}_{3 i}^{gt} + 1 \end{matrix} & \begin{matrix} i \in V^{'}, t \in T \\ g \in SL ∖ {0} \end{matrix} \end{matrix}$ (69) $\begin{matrix} \sum_{g \in SL} d_{i}^{* gt} + {SV}_{4 i}^{t} \geq {de}_{real i}^{t} & i \in V^{'}, t \in T \end{matrix}$ (70) $\begin{matrix} + \sum_{g \in SL} d_{i}^{* gt} \leq {de}_{real i}^{t} + {SV}_{5 i}^{t} & i \in V^{'}, t \in T \end{matrix}$ (71) $\begin{matrix} \sum_{i \in V^{'}} \sum_{g \in SL} q_{i}^{* gkt} \leq {Qu}_{real} . y_{0}^{* kt} + {SV}_{6}^{kt} & k \in K, t \in T \end{matrix}$ (72) ${SV}_{1}^{t}, {SV}_{2}^{t} \geq 0 t \in T$ (73) $\begin{matrix} {SV}_{3 i}^{gt} \geq 0 & \begin{matrix} i \in V^{'}, t \in T \\ g \in SL ∖ {0} \end{matrix} \end{matrix}$ (74) ${SV}_{4 i}^{t}, {SV}_{5 i}^{t} \geq 0 i \in V^{'}, t \in T$ (75) ${SV}_{6}^{kt} \geq 0 k \in K, t \in T$ (76)

In the above formulation, SV variables are the only decision variables and determine the violation of chance constraints under random realization. To evaluate the proposed models, the average and standard deviation of objective function values of the above linear formulation is used as performance measures. The obtained results are shown in Table 5.

To infer that Final-RPP outperforms the CCPCCP models from the data provided in Table 5, a series of statistical analyses are conducted in the following.

Firstly, a Kolmogorov-Smirnov test is carried out to ensure that the data provided in Table 5 are approximately normally distributed.

Afterwards, a one-way analysis of variance (ANOVA) is employed to see whether there are any significant differences between the means of the different groups (Averages of CCPCCP (ω_i = 0.7), CCPCCP (ω_i = 0.8), CCPCCP (ω_i = 0.9) and Final-RPP). The obtained results showed that the differences between the means are significant (See Table 6). To discover the source of the differences, post-hoc tests (e.g., the Fisher Least Significant Difference (LSD) method and Tukey test) should be implemented. Thus, a Levene’s test, as a prerequisite for post-hoc tests, is performed to make sure that the homogeneity of variances is satisfactory. The outcomes of Levene’s test showed that the differences are not significant (p_value = 0.736).

Finally, both LSD and Tukey tests are employed. These results are provided in Table 7.

In Table 7, the bold numbers relate to those mean differences that are significant at 0.01 confidence level. One can see that in both LSD and Tukey cases, it is only Final-RPP model that significantly differs from the other groups.

In a nutshell, the following conclusions from the numerical results provided in this section can be derived:

Final-RPP is significantly better than its rivals in terms of quality of solutions.

Final-RPP has the least objective function value among the 4 proposed models (Table 5).

The standard deviation for Final-RPP objective function value is openly shorter than that of CCPCCP models. It shows that Final-RPP provides more robust solutions in uncertain and volatile situations.

From solution time’s perspective, CCPCCP with ω_i = 0.7 is the best model, while CCPCCP with ω_i = 0.8 shows the weakest performance.

Considering the high quality and robust solutions Final-RPP generates, its performance in terms of solution times is satisfactory.

Summarizing the above points, Final-RPP is the most preferred method to cope with the inherent uncertainty of imprecise parameters. Particularly, in blood-supply chains, where we encounter a deep uncertain environment in terms of demand fluctuations at both supply and demand points, unknown costs of wastage units due to its fatal consequences and lack of knowledge about the exact holding costs, the significance of acquiring more robust solutions hence utilizing Final-RPP model is magnified.

7 Conclusions and directions for further research

In this paper, the problem of delivering blood units from a community blood center to a number of hospital blood banks under uncertain environment was tackled for the first time in the literature, its major challenges were discussed and an MILP formulation was developed as its solution approach. Since blood supply chains are overwhelmed with uncertain parameters, there is a need of developing suitable approaches to cope with these uncertainties.

Thus, in this paper, a possibilistic chance constrained programming (PCCP) model and a robust possibilistic programming approach was proposed.

To solve these models to optimality, a novel branch-and-cut method was developed. Afterwards, the data driven from a real blood supply chain was employed to assess the performance of both PCCP and RPP approaches. Also, a number of random realizations were produced and the performances of both approaches were evaluated against them using ANOVA technique and LSD and Tukey tests. The results showed that Final-RPP is significantly better than CCPCCP formulations under different realizations. Moreover, Final-RPP had the least objective function value not only on the average but also on every single one of the 15 realizations. It was also observed that the performance of both approaches was satisfactory from solution times viewpoint.

This paper can be considered as a pioneering work in the context of Blood-IRPs. Nonetheless, area for improving the current research is widely open. For instance, considering partial delivery in the distribution network will be of great interest. Also, blood demands can be categorized in two elective and emergency groups while stock-outs can be permitted for some elective demands. The topic of heterogeneous vehicles is also an interesting area.

References

Bell

W.J.

, Dalberto

L.M.

, Fisher

M.L.

, Greenfield

A.J.

, Jaikumar

, Kedia

, Mack

R.G.

and Prutzman

P.J.

, Improving the distribution of industrial gases with an on-line computerized routing and scheduling optimizer, Interfaces 13(6) (1983), 4–23.

Lmariouh

, Hachemi

N.E.

, Jamali

and Bouami

, The Multi-vehicle Mutli-Product Inventory-Routing Problem: A case study, Paper Presented at the Advanced Logistics and Transport (ICALT), 2014 International Conference on 2014–.

Shen

, Chu

and Chen

, A Lagrangian relaxation approach for a multi-mode inventory routing problem with transshipment in crude oil transportation, Computers & Chemical Engineering 35(10) (2011), 2113–2123.

Vidović

, Popović

and Ratković

, Mixed integer and heuristics model for the inventory routing problem in fuel delivery, International Journal of Production Economics 147(Part C) (2014), 593–604.

Aksen

, Kaya

, Salman

F.S.

and Tüncel

Ö.

, An adaptive large neighborhood search algorithm for a selective and periodic inventory routing problem, European Journal of Operational Research 239(2) (2014), 413–426.

Jagannathan

and Sen

, Storing crossmatched blood: A perishable inventory model with prior allocation, Management Science 37(3) (1991), 251–266.

Aelion

, Castells

and Veroutis

, Life cycle inventory analysis of chemical processes, Environmental Progress 14(3) (1995), 193–200.

Blackburn

and Scudder

, Supply chain strategies for perishable products: The case of fresh produce, Production and Operations Management 18(2) (2009), 129–137.

Farahani

, Grunow

and Günther

H.-O.

, Integrated production and distribution planning for perishable food products, Flexible Services and Manufacturing Journal 24(1) (2012), 28–51.

10.

Coelho

L.C.

and Laporte

, Optimal Joint Replenishment and Delivery of Perishable Products. Report, 2013a.

11.

, Diabat

, Richard

J.-P.

and Yih

, A column generation-based heuristic algorithm for an inventory routing problem with perishable goods, Optimization Letters 7(7) (2013), 1481–1502.

12.

Al Shamsi

, Raisi

A.A.

and Aftab

, Pollution-Inventory Routing Problem with Perishable Goods, In Logistics Operations, Supply Chain Management and Sustainability, Springer, 2014, pp. 585–596.

13.

Coelho

L.C.

and Laporte

, Optimal joint replenishment, delivery and inventory management policies for perishable products, Computers & Operations Research 47 (2014), 42–52.

14.

Gunpinar

, Supply Chain Optimization of Blood Products. University of South Florida, 2013.

15.

Pierskalla

W.P.

, Supply Chain Management of Blood Banks. In Operations Research and Health Care, edited by Brandeau

Margaret L

, Sainfort

François

and Pierskalla

William P

, Springer US, pp. 103–145.

16.

Solyali

, Cordeau

J.-F.

and Laporte

, Robust inventory routing under demand uncertainty, Transportation Science 46(3) (2012), 327–340.

17.

Gabrel

, Murat

and Thiele

, Recent advances in robust optimization: An overview, European Journal of Operational Research 235(3) (2014), 471–483.

18.

Pishvaee

M.S.

, Razmi

and Torabi

S.A.

, Robust possibilistic programming for socially responsible supply chain network design: A new approach, Fuzzy Sets and Systems 206 (2012), 1–20.

19.

Zahiri

, Tavakkoli-Moghaddam

, Mohammadi

and Jula

, Multi-objective design of an organ transplant network under uncertainty, Transportation Research Part E: Logistics and Transportation Review 72 (2014), 101–124.

20.

Verdegay

J.L.

, Fuzzy mathematical programming, Fuzzy Information and Decision Processes 231 (1982), 237.

21.

Zahiri

, Tavakkoli-Moghaddam

and Pishvaee

M.S.

, A robust possibilistic programming approach to multi-period location–allocation of organ transplant centers under uncertainty, Computers & Industrial Engineering 74 (2014), 139–148.

22.

Torabi

S.A.

and Hassini

, An interactive possibilistic programming approach for multiple objective supply chain master planning, Fuzzy Sets and Systems 159(2) (2008), 193–214.

23.

Archetti

, Bertazzi

, Laporte

and Speranza

M.G.

, A branch-and-cut algorithm for a vendor-managed inventory-routing problem, Transportation Science 41(3) (2007), 382–391.

24.

Coelho

L.C.

and Laporte

, A branch-and-cut algorithm for the multi-product multi-vehicle inventory-routing problem, International Journal of Production Research 51(23-24) (2013b), 7156–7169.

25.

Coelho

L.C.

and Laporte

, The exact solution of several classes of inventory-routing problems, Computers & Operations Research 40(2) (2013c), 558–565.

26.

Adulyasak

, Cordeau

J.-F.

and Jans

, Formulations and branch-and-cut algorithms for multivehicle production and inventory routing problems, INFORMS Journal on Computing 26(1) (2014), 103–120.