Introduction

Blood is an essential element of human life and required for different treatments including anemia treatments, cancer, organ transplants, and major surgeries [1, 2]. Nearly 30 million units of blood products are employed in blood transfusion each year, while this amount is expected to increase because of an increase in life expectancy and advances in medical procedures requiring blood transfusion [3, 4]. Hence, improving productivity of blood management systems is a key concern. In this regard, supply, inventory, and routing for blood products are important problems for healthcare managers, since such decisions have significant consequences for society and patients’ health [5, 6]. This highlights the need for further collaboration across the whole supply chain network and integrated decision-making at both strategic and operational levels to reduce delivery time and increase customer service level [7] while reducing product waste and environmental, social, and economic consequences of perishability.

The integrated supply and distribution management of blood is more sophisticated than the other products. This includes making joint decisions on the number of products to be supplied and distributed to blood centers, the level of inventory to be held, and the vehicle routing [8]. However, the efficient coordination of these decisions is a challenging task due to product shortage and waste, especially when it comes to uncertain situations [9,10,11]. Therefore, supply and inventory decisions need to be taken in such a way that demand is fully satisfied with the minimum waste [12, 13]. In spite of the significance of this problem, it has not been adequately addressed in the literature and there is a lack of studies on integrated inventory-routing of blood products to hospitals under demand uncertainty. The available literature has not taken into account transshipment among blood centers to investigate how it may improve inventory management of blood products [14]. To fill these gaps, in this paper, a mathematical modeling approach is adopted. A two-echelon supply chain consisting of a single supplier and a set of blood centers is considered, where the supplier has the job of supplying blood to the blood centers. The objective is to determine the optimal amounts of supplying and inventory, delivery plan and time, and vehicle routing over a multi-period planning horizon. While the supplier seeks to optimize quantity of supply, there may be shortage or excess inventory at the blood centers due to uncertainty of customer demand. Thus, transshipment of blood from the blood centers having excess inventory to those with the shortage is allowed to decrease total costs including shortage, waste, and transportation costs. The problem is addressed using a mathematical modeling approach and robust optimization is adopted to deal with uncertainty in demand. Since the proposed model is NP-hard, a heuristic solution algorithm is developed. Numerical experiments and data from a real case of blood products are presented to illustrate the efficiency of the proposed model and solution algorithm.

This paper is organized as follows: Sect. “Literature review” presents a review of the relevant literature. Section “Proposed mathematical model” develops the mathematical model for deterministic and uncertain situations. The proposed solution algorithm is presented in Sect.“ Proposed solution algorithm”. Section “Computational experiments” discusses numerical results based on available data in the literature and the data from a real case of blood products. Finally, research conclusions and some recommendations for future research are provided in Sect. Conclusions.

Literature review

Most previous studies on production and distribution management were conducted with non-perishable items [2, 15]. However, recently, researchers have become more interested in studying the production and distribution management of perishable goods due to the complexity inherent in the problem [16]. Li et al. [17] investigated a novel integrated planning model for intelligent food logistics systems with the objectives of minimizing total production, inventory, and transportation costs and maximizing average food quality. The results indicated that their proposed approach could reduce the total costs by 10.77% compared to an existing three-phase heuristic method. In another research, Marandi and Zegordi [18] considered costs and assessed how the quality of blood products could be improved by shortening the time interval between production and distribution. The researchers dealt with a variation of Integrated Production Distribution System (IPDS) that contains a product with short shelf-life.

In addition to the mentioned studies, Li et al. [19] formulated the production-inventory-routing planning with an integrated Mixed Integer Linear Programming (MILP) model, where the food quality level is explicitly traced throughout the supply chain. The results showed that the problem is so complex that it is hard to solve even in a small case with 15 retailers, three periods, and two vehicles.

Some other authors conducted research on Integrated Production-Inventory-Routing Problem (IPIRP) by focusing on different aspects. Agra et al. [20] considered a stochastic single item IPIRP with a single producer and multiple clients. Clients were allowed to ask for a penalty cost as a backlog when demands were considered uncertain. Preliminary tests based on random instances were conducted illustrating that the hybrid heuristic had a better performance than the classical sample approximation algorithm for hard instances.

Qiu et al. [21] presented a generalized IPIRP model with perishable inventory. The authors analyzed the optimal integrated decisions of when and how much to deliver and sell products with varying manufacturing periods. The researchers discussed main inventory management policies to demonstrate the applicability of the model in real-world applications for Production-Routing Problems (PRPs) with perishable inventory. Widyadana and Irohara. [22] extended the IPIRP by considering recovery and remanufacturing of End-Of-Life (EOL) products as stable factors. They aimed at using a mathematical programming approach to determine the optimal amounts of goods (i.e., lot-sizes) required to be produced, reproduced, and stored while meeting customers' requests (i.e., deliveries and pick-ups) with the least total costs as a result of the integrated operations. The results of the study demonstrated a need for development of approximate resolution methods in large instances.

Neves-Moreira et al. [23] utilized a new three-phase methodology to solve a large Production-Routing Problem (PRP) by combining realistic features. The large-size instances caused intractable problems, which needed to be resolved via efficient solution methods. In 2018, Agra et al. [20] carried out a study on a single item IPIRP with a single producer/supplier and multiple retailers. Inventory management constraints were taken into account both at the producer and at the retailers following a vendor-managed inventory approach and decision on the replenishment policy for each retailer.

Overall, to identify our contributions, the literature review reveals that the integrated decision-making on supply, inventory, and routing of blood products has not been widely studied. In addition, most of the previous studies on supply and distribution management of blood products have presented a replenishment policy for blood centers, i.e., the place, where blood products are produced. Nevertheless, routing to hospitals has not been taken into account in previous works, whereas hospitals can be viewed as the retailers having the job of satisfying patients’ demand, i.e., customers. Additionally, in previous studies, only the supplier was assumed to be responsible for meeting blood centers’ demand and transshipment flows among retailers were not allowed. However, in practical situations in the real world, transshipment can be applied as an effective strategy to reduce the overall shortage and excess inventory under uncertainty. While some researchers have investigated replenishment policy for blood products in hospitals, only a few researchers such as Hemmelmayr et al. [24] have addressed the routing problem in the distribution of blood products to hospitals, which is indicative of a wide gap in the available literature. The literature review also reveals that integrated inventory-routing of blood products has not been investigated in uncertain conditions. The lack of information and research in this regard has been stressed by Adulyasak et al. [25] who argued that no studies have been conducted on robust optimization of the integrated production-inventory problem. Moreover, Coelho et al. [26] provide a comprehensive review of this literature, based on a new classification of the problem. The authors categorize IRPs with respect to their structural variants and the availability of information on customer demand.

Reviewing past research works as well as realizing the mentioned gaps in the literature motivated the researchers to investigate the problem under more realistic conditions. The current paper contributes to the available literature by developing a model for integrating supply, inventory, and routing of blood with uncertain demand. The model proposed in this study allows for transshipping products among blood centers to help improve the management of shortage and excess inventories. The robust optimization approach is adopted to deal with uncertainty and a heuristic algorithm is developed to obtain high-quality solutions. The data from blood distribution services are used to further validate the proposed model.

Proposed mathematical model

We suppose a two-echelon supply chain consisting of one supplier of a blood product, i.e., blood, and a set of blood centers over a finite time horizon. The supplier has a limited capacity in each period and decides on supplying quantity at the beginning of each period. There is a single-vehicle with limited capacity that traverses a route to deliver products from the supplier to a subset of blood centers at the beginning of each period. Transshipment among blood centers and direct delivery from supplier to the blood centers is allowed during each period. Figure 1 shows a graphical illustration of the network with four blood centers over three periods. In the first period, blood products are shipped from the supplier to blood centers 1 and 2, respectively, and there is a transshipment from blood center 2 to blood center 3. In the second period, blood is directly delivered to blood center 4 and then to blood center 3, while we have blood transshipment from blood center 3 to blood center 2. Finally, in period 3, there is direct shipment from the supplier to blood center 4, while both blood centers 2 and 1 receive blood products via transshipment.

Fig. 1
figure 1

Graphical illustration of the problem with four blood centers over three time-periods

Full size image

Other assumptions are made as follows:

  • The inventory level of products in each age group at both the supplier and the blood center is known at the beginning of the planning horizon [27, 28].

  • The blood centers have limited capacities [29, 30].

  • The supplying capacity in each period is known [31, 32].

  • The blood centers follow the Fresher First (FF) policy to satisfy customers’ demand. Adopting this policy, rather than the relaxed model, can decrease total costs including shortage, waste, and transportation costs of the whole supply chain. Fresh First (FF) policy by which the retailer always sells the fresher items first. This policy ensures a longer shelf life and increases utility for the customers but, at the same time, yields a higher spoilage rate [13]. In addition, it ensures a longer shelf life and increases goods utility for the customers but, at the same time, it yields a higher spoilage rate [13].

The model proposed by Coelho et al. [33] was considered as the base model and extended by incorporating product age, FF policy, and supplying planning decisions. To do so, we first investigated the most relevant models presented in previous works, including Qiu et al. [34], Elalouf et al. [35], and Hosseinifard and Abbasi [36]. Choosing these models was because td with optimizing integrated production and distribution decisions making similar or neighboring assumptions and objectives. The mathematical model is then formulated using the initial set of parameters and decision variables extracted from these models and further lopment othe objective function and constraints by taking into account the unique characteristics of the addressed problem and observations from the real-case of the blood distribution system. The sets and indices, parameters, and decision variables used in this study are defined as follows:

Sets and indices:

\(I\)

Set of the supplier and blood centers \(i \in I\)

\(J\)

Set of blood centers \(j \in J\)

\(G\)

Set of the ages of blood products \(g \in G\)

\(S\)

Set of scenarios \( s \in S\)

\(T\)

Set of time periods \(t \in T\)

Parameters:

\(b_{ij}\)

Transshipment cost between node \(i\) and \(j\)

\(b^t\)

Supply capacity of the supplier at period \(t\)

\(c_{ij}\)

Transportation cost of vehicle between nodes \(i\) and \(j\)

\(C_i\)

Capacity of holding inventory at the \(i{\text{th}}\) blood center

\(a_i^t\)

Demand of \(i{\text{th}}\) blood center at period \(t\)

\(a_i^{t,s}\)

Realization of \(a_i^{t,s}\) under scenario \(s\)

\(F^t\)

Total supply at period \(t\)

\(h_i\)

Cost of holding inventory at the supplier and blood center

\(a_i^{0,g}\)

Initial inventory of age \(g\) at the supplier and \(i{\text{th}}\) blood center at the beginning of the first period

\(M\)

Maximum allowable loss of product

\(p_s\)

Probability of scenario \(s\)

\(P\)

Supplying cost of each unit of blood product

\(Q\)

Capacity of vehicle

\(\lambda\)

A weight factor for the variance of the objective function

\(\omega\)

A weight factor for balancing between solution robustness and model robustness

Decision variables:

\({d}_{i}^{t,g}\)

Portion of demand in \(i{\text{th}}\) blood center which is satisfied by the products of age \(g\) at period \(t\)

\({d}_{i}^{t,g,s}\)

Realization of \({d}_{i}^{t,g}\) under scenario \(s\)

\({l}_{i}^{t,g}\)

An auxiliary binary variable for adopting the FF policy; 1 if a product of age \(g\) is sold by blood center \(i\) at period \(t\), and 0 otherwise

\({d}_{i}^{t,g,s}\)

Realization of \({l}_{i}^{t,g}\) under scenario \(s\)

\({I}_{i}^{t,g}\)

Inventory level of age \(g\) at the supplier and blood centers at the end of period \(t\)

\({I}_{i}^{t,g,s}\)

Realization of \({I}_{i}^{t,g}\) under scenario \(s\)

\({v}_{i}^{t}\)

An auxiliary variable for sub-tour elimination

\({x}_{i,j}^{t}\)

A binary variable; 1 if node \(i\) is visited after node \(j\) in period \(t\) and 0 otherwise

\({q}_{i}^{t,g}\)

Quantity of blood product of age \(g\) delivered to blood center \(i\) at period \(t\)

\({w}_{ij}^{t,g}\)

Quantity of blood product with age \(g\) transshipped from node \(i\) to j at period \(t\)

\({w}_{ij}^{t,g,s}\)

Realization of \({w}_{ij}^{t,g}\) under scenario \(s\)

\({\delta }_{i}^{t,s}\)

Realization of blood product demand (shortage) in \(i{\text{th}}\) blood center under scenario \(s\)

Based on the notations provided above, the deterministic mathematical model is formulated as follows:

$$ {\text{Min}} \sum \limits_{i \in V} \sum \limits_{t \in T} \sum \limits_{g \in G} h_i I_{i,t}^g + \sum \limits_{i \in V} \sum \limits_{j \in V} \sum \limits_{t \in T} c_{ij} x_{ij}^t + \sum \limits_{i \in V} \sum \limits_{j \in V^{\prime}} \sum \limits_{t \in T} \sum \limits_{g \in G} b_{ij} w_{i,j}^{t,g} + P \sum \limits_{t \in T} F^t . $$
(1)

Subject to:

$$ I_0^{t,g} = I_0^{t - 1,g - 1} - \sum \limits_{i \in V^{\prime}} q_i^{t,g} - \sum \limits_{i \in V^{\prime}} w_{0,i}^{t,g} t \in T, g = 2,3, \ldots H, $$
(2)
$$ I_0^{t,g} = F^t - \sum \limits_{i \in V^{\prime}} q_i^{t,g} - \sum \limits_{i \in V^{\prime}} w_{0,i}^{t,g} t \in T, g = 1 $$
(3)
$$ I_0^{t,g} = I_i^{t - 1,g - 1} + q_i^{t,g} + \sum \limits_{j \in V} w_{i,j}^{t,g} - \sum \limits_{j \in V^{\prime}} w_{i,j}^{t,g} - d_i^{t,g} t \in T, g = 2,3, \ldots H, i \in I $$
(4)
$$ I_i^{t,g} = q_i^{t,g} + \sum \limits_{j \in V} w_{j,i}^{t,g} - \sum \limits_{j \in V^{\prime}} w_{i,j}^{t,g} - d_i^{t,g} t \in T, g = 1, i \in I $$
(5)
$$ \sum \limits_{t \in T} \sum \limits_{i \in V} \left( {l_i^{t,g} } \right) \le M \quad g = H $$
(6)
$$ F^t \le B^t \quad t \in T $$
(7)
$$ a_i^t = \sum \limits_{g \in G} d_i^{t,g} \quad t \in T, i \in I, $$
(8)
$$ \left( {a_i^{t,g} } \right) \le C_i l_i^{t,g} \quad t \in T, g \in G, i \in I $$
(9)
$$ \left( {l_i^{t,g - 1} } \right) \le l_i^{t,g} \quad t \in T, g = 2,3, \ldots ,H, i \in I $$
(10)
$$ C_i \left( {1 - l_i^{t,g - 1} } \right) \ge \sum \limits_{g^{\prime} = g} \left( {I_i^{t - 1,g^{\prime} - 1} + q_i^{t,g^{\prime}} + \sum \limits_{j \in V} w_{j,i}^{t,g^{\prime}} - \sum \limits_{j \in V} w_{i,j}^{t,g^{\prime}} } \right) - a_i^t + 1 \quad t \in T, g = 2,3, \ldots ,H, i \in I $$
(11)
$$ \sum \limits_{g \in G} I_i^{t,g} \le C_i \quad t \in T, i \in I $$
(12)
$$ q_i^{t,g} \le C_i \sum \limits_{j \in V} x_{i,j}^t \quad t \in T, i \in I, g \in G $$
(13)
$$ \sum \limits_{g \in G} q_i^{t,g} \le C_i - \mathop \sum \limits_{g = 1}^{H - 1} I_i^{t - 1,g} \quad t \in T, i \in I $$
(14)
$$ \sum \limits_{g \in G} \sum \limits_{i \in V^{\prime}} q_i^{t,g} \le Q \quad t \in T $$
(15)
$$ \sum \limits_{j \in J} x_{i,j}^t = \sum \limits_{j \in J} x_{ji}^t \quad t \in T, i \in I $$
(16)
$$ v_i^t - v_j^t + Qx_{ij}^t \le Q - \sum \limits_{g \in G} q_j^{t,g} \quad t \in T, j \in J $$
(17)
$$ \sum \limits_{g \in G} q_j^{t,g} \le v_i^t \le Q\quad i \in I $$
(18)
$$ \sum \limits_{i \in v} x_{i0}^t \le 1 $$
(19)
$$ x_{ij}^t \in \left\{ {0,1} \right\}\quad t \in T, i \in I,j \in J,i \ne j $$
(20)
$$ I_i^{t,g} \ge 0, w_{i,j}^{t,g} \ge 0 d_{i,j}^{t,g} \ge 0,v_i^t \ge 0,F^t \ge 0,q_i^{t,g} \ge 0\quad t \in T, g \in G, i \in I. $$
(21)

The problem setting aims at determining the optimal amounts of supply and inventory, and the optimal delivery plan and time, as well as vehicle routing over a multi-period planning horizon. The objective function (1) minimizes total costs defined in four terms. The first term implies total inventory costs comprising of holding costs and initial inventory costs. The second term implies the expected routing cost based on the selected route. The third term presents the amount of blood transshipped in the network and, finally, the last term shows total supplying costs based on the cost of each unit of supply at each period. This formulation is consistent with the functions defined for total cost minimization in previous studies regarding integrated inventory-routing problems. It is also well-suited to meet the requirements of blood distribution system in the real situation.

Taking into account the main characteristics of blood products, constraint (2) updates the inventory level of age two and above at the supplier. In addition, constraint (3) implies the inventory level of age one at the supplier based on the supply rate. Constraints (4) and (5) indicate the inventory level of age two and above as well as the inventory level of age 1 at the blood centers, respectively. Constraint (6) ensures the maximum allowable loss of products and constraint (7) shows the supply capacity. Constraint (8) implies that products with different ages can be used to satisfy demand.

Constraints (9)–(11) are used for adopting the FF policy for blood products. Constraint (9) implies that when demand for a product of age \(g\) is satisfied, the relevant auxiliary variable (\({l}_{i}^{t,g}\)) will be equal to one. Constraint (10) ensures that the inventory of age \(g-1\) is used if the inventory of age \(g\) is not enough for satisfying demand. In addition, constraint (11) ensures that the inventory of age \(g-1\) cannot be used when available inventory of age g and above is adequate for satisfying demand.

Constraint (12) shows the limited capacity of blood centers for holding inventory. Constraint (13) guarantees that the product is delivered to the blood center if it is visited by the vehicle and constraint (14) shows that the amount of product delivered to the blood center should be within its available capacity. Constraint (15) denotes the limited capacity of the vehicle. Constraint (16) constructs the tours and constraints (17) and (18) are the sub-tour elimination constraints. Constraint (19) implies that only a single vehicle is available. Finally, constraints (20) and (21) show the types of decision variables.

Robust optimization model

To account for demand uncertainty in the aforementioned problem, a robust optimization approach is adopted. The robust optimization is a powerful and effective approach for handling optimization problems and decision-making under uncertainty [37, 38]. It has become widely popular and can been applied to various supply chain planning and optimization problems in uncertain conditions. In this paper, we adopt the same robust optimization approach and measure as the one used in several other relevant studies that are similar in nature (e.g., Pishvaee et al. [39], Rabbani et al. [40], Wee et al. [41]; Hendalianpour et al. [42, 43], Wu [44, 45]). According to Pishvaee et al. [39], such an approach is capable of controlling the degree of feasibility and optimality robustness and also establishing a reasonable tradeoff between the robustness, cost of robustness, and other objectives such as improving the average performance of the system. Thus, it provides a flexible approach to decision making under uncertainty and is applicable in most business cases. Suppose the presented deterministic model as follows [39, 41]:

$$ \min \xi = c^{\text{T}} x + d^{\text{T}} $$
(22)
$$ Ax = b $$
(23)
$$ Bx + Cy = e $$
(24)
$$ x,y \ge 0. $$
(25)

The Eq. (22) represents deterministic objective function, where \(x\) and \(y\) denote design and control variables. Constraint (23) shows design constraint in which A and b parameters are known. Constraint 24 implies control constraint in which \(B,C\) and \(e\) parameters are uncertain. Suppose the scenario set \(\left\{ {1,2, \ldots ,{\Omega }} \right\}\) with the probability of \(\sum \limits_{s \in S} p_s = 1\) for each scenario. Then, \(B_s ,C_s ,e_s ,d_s\) parameters denote uncertain coefficients in the above model. The original model will be then converted to the following robust model:

$$ \min \sigma \left( {x,y_1 ,y_2 , \ldots ,y_s } \right) + \omega \rho \left( {\delta_1 ,\delta_2 , \ldots ,\delta_n } \right) $$
(26)
$$ Ax = b $$
(27)
$$ B_s x + C_s y_s + \delta_s = e_s\quad s \in S $$
(28)
$$ x \ge 0, y_s \ge 0\quad s \in S. $$
(29)

\(\delta_s\) in constraint (28) denotes the error vector which measures deviation from control constraint in each scenario. The first part of the objective function (26) indicates solution robustness and tendency of decision maker to achieve a lower cost. The second part indicates model robustness and incurs a penalty to solutions that deviate from control constraints. The \(\omega\) coefficient is used for tradeoff between solution robustness and model robustness. According to El Ghaoui [46], the \(\sigma\) function is constructed as follows:

$$ \sigma \left( {x,y_1 ,y_2 , \ldots ,y_s } \right) = \sum \limits_{s \in S} p_s \xi_s + \lambda \sum \limits_{s \in S} p_s \left( {\left( {\xi_s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} \xi_{s^{\prime}} } \right)} \right)^2 . $$
(30)

In the above function, the first part shows weighted average of the objective function and the second part denotes sum of deviations from the objective function under different scenarios, i.e., variance of the objective function. Yu and Li [47] and Haijema et al. [48] used the following formulation to linearize the robust model.

$$ \sigma \left( {x,y_1 ,y_2 , \ldots ,y_s } \right) = \sum \limits_{s \in S} p_s \xi_s + \lambda \sum \limits_{s \in S} p_s \left| {\left( {\xi_s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} \xi_{s^{\prime}} } \right)} \right| $$
(31)
$$ \min \sum \limits_{s \in S} p_s \xi_s + \lambda \sum \limits_{s \in S} p_s \left( {\left( {\xi_s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} \xi_{s^{\prime}} } \right)} \right) + 2\theta_s $$
(32)
$$ \xi_s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} \xi_{s^{\prime}} + \theta_s \ge 0, \quad s \in \theta_s $$
(33)
$$ \theta_s \ge 0 \quad s \in S. $$
(34)

The robust counterpart of the model (26)–(29) is then written in the following way:

$$ \min \sum \limits_{s \in S} p_s \xi_s + \lambda \sum \limits_{s \in S} p_s \left( {\left( {\xi_s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} \xi_{s^{\prime}} } \right) + 2\theta_s } \right) + \omega \sum \limits_{s \in S} p_s \delta_s . $$
(35)

subject to the constraints (33) and (34) and (27)–(29). To develop the robust model, we supposed the components of the objective function as follows:

$$ {\text{PCO}} = P \sum \limits_{t \in T} F^t $$
(36)
$$ {\text{RC}} = \sum \limits_{i \in V} \sum \limits_{t \in T} \sum \limits_{g \in G} c_{ij} x_{ij}^t $$
(37)
$$ Z^s = \sum \limits_{i \in V^{\prime}} \sum \limits_{t \in T} \sum \limits_{g \in G} h_i I_{i,t}^{g,s} + \sum \limits_{i \in V} \sum \limits_{j \in V^{\prime}} \sum \limits_{t \in T} \sum \limits_{g \in G} b_{ij} w_{ij}^{t,g,s} . $$
(38)

The objective function of the robust model is then constructed in the following way:

$$ \min {\text{RC}} + {\text{PCO}} + \sum \limits_{s \in S} p_s Z^s + \lambda \sum \limits_{s \in S} p_s \left( {\left( {Z^s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} Z^{s^{\prime}} } \right) + 2\theta_s } \right) + \omega \sum \limits_{s \in S} \sum \limits_{i \in V} \sum \limits_{t \in T} p_s \delta_i^{t,s} , $$
(39)

subject to constraints (7), (13), (15)–(21) and the following constraints:

$$ I_0^{t,g,s} = I_0^{t - 1,g - 1,s} + \sum \limits_{i \in V^{\prime}} q_i^{t,g,s} - \sum \limits_{i \in V^{\prime}} w_{0,i}^{t,g,s} t \in T, g = 2,3, \ldots H,, s \in S $$
(40)
$$ I_0^{t,g,s} = F^t + \sum \limits_{j \in V^{\prime}} q_i^{t,g,s} - \sum \limits_{j \in V^{\prime}} w_{0,i}^{t,g,s} t \in T, g = 1, s \in S $$
(41)
$$ I_0^{t,g,s} = I_0^{t - 1,g - 1,s} + q_i^{t,g,s} + \sum \limits_{j \in V^{\prime}} w_{i,j}^{t,g,s} - \sum \limits_{j \in V^{\prime}} w_{i,j}^{t,g,s} - d_i^{t,g,s} t \in T, g = 2,3, \ldots H, i \in I, s \in S $$
(42)
$$ I_0^{t,g,s} = q_i^{t,g} + \sum \limits_{j \in V^{\prime}} w_{j,i}^{t,g,s} - \sum \limits_{j \in V^{\prime}} w_{i,j}^{t,g,s} - d_i^{t,g,s} t \in T, g = 1, i \in I, s \in S $$
(43)
$$ \sum \limits_t \sum \limits_i \left( {l_i^{t,g - 1,s} } \right) \le L\quad t \in T, g = H, i \in I, s \in S $$
(44)
$$ a_i^{t,s} = \sum \limits_{g \in G} d_i^{t,g,s} + \delta_i^{t,s}\quad t \in T, i \in I, s \in S $$
(45)
$$ \left( {d_i^{t,g,s} } \right) \le C_i I_i^{t,g,s}\quad t \in T, g \in G, i \in I, s \in S $$
(46)
$$ \left( {l_i^{t,g - 1,s} } \right) \le I_i^{t,g,s} t \in T, g = 2,3, \ldots ,H, i \in I, s \in S $$
(47)
$$ C_i \left( {1 - l_i^{t,g - 1,s} } \right) \ge \sum \limits_{gg = g} \left( {I_i^{t - 1,gg - 1,s} + q_i^{t,gg} + \sum \limits_{j \in v} w_{i,j}^{t,gg,s} } \right) - a_i^{t,s} + 1 t \in T, g = 2,3, \ldots ,H, i \in I, s \in S $$
(48)
$$ \mathop \sum \limits_{{g \in G}} I_{i}^{{t,g,s}} \le C_{i} \quad t \in T,~\,i \in I,~\,s \in S $$
(49)
$$ \sum \limits_{g \in G} q_i^{t,g} \le C_i - \mathop \sum \limits_{g = 1}^{H - 1} I_i^{t - 1,g,s} \quad t \in T, i \in I, s \in S $$
(50)
$$ Z^s - \sum \limits_{s \in S} p_s Z^s + \theta_s \ge 0 \quad s \in S $$
(51)
$$ I_i^{t,g,s} \ge 0, w_{i,j}^{t,g,s} \ge 0 d_{i,j}^{t,g,s} ,\theta_s \ge 0 \quad t \in T, g \in G, i \in I, s \in S. $$
(52)

Constraint (40) updates the inventory level of products of age two and above at the supplier in each scenario. Constraint (41) shows the inventory level of age one at the supplier. Constraints (42) and (43) update the inventory level of blood centers for products of age 2 and above as well as products of age 1. Constraint (44) shows the maximum allowable loss in each scenario. Constraint (45) defines the relationship between demand, the sum of inventory of different ages, and the shortage amount. Constraints (46)–(48) indicate that the FF policy is adopted. Constraint (49) implies the limited capacity of the blood centers, and constraint (50) implies that the quantity delivered by the vehicle cannot exceed the available capacity at the blood centers. Constraint (51) gives the \({\theta }_{s}\) value in each scenario based on constraint (32). Constraint (52) define types of decision variables.

Proposed solution algorithm

The model developed in this paper significantly extends the work by Coelho et al. [33] by incorporating product age, FF policy, and supply planning decisions. The model presented by Coelho and Laporte [13] is NP-hard due to the routing sub-problem, which consists of a vehicle with limited capacity. Thus, the model developed in this paper is also NP-hard and cannot be solved in polynomial time for large-scale problems. In addition, the robust optimization approach increases complexity of the original model and solution time of exact algorithms. Therefore, in this paper, a heuristic solution algorithm is developed for the presented mathematical model.

The proposed solution algorithm begins with an initial solution, which is improved through an iterative cycle including two local search strategies, namely inserting the best and removing the worst solutions [49, 50]. Moreover, a shocking stage is incorporated to prevent falling in local optima. The structure of the proposed solution algorithm is illustrated in Fig. 2. The steps of the proposed algorithm are given as follows:

  • The parameter values, including the maximum number of iterations, allowable computation time, and proportion of deleting stops, as well as that of solution numbers in the solution pool, are identified.

  • An initial solution is generated and put in the solution pool to weigh the initial solution as the best solution in the first step of the algorithm.

  • A solution is chosen from the solution pool and then improved using the neighbourhood structure. If the new solution serves better than the previously-used solutions, then it is replaced.

  • The new obtained solution is replaced with the worst solution in three modes: (1) the completion of the solution pool (2) lack of existence of the current solution in the solution pool, and (3) inefficiency of the available solution in the solution pool compared to the current solution.

  • If the solution pool is not completed and does not have a similar solution, then the current solution is added to the solution pool. The algorithm is implemented until one of the two conditions is satisfied: (1) the maximum number of iterations is reached, or (2) the maximum allowable computation time is established. After stopping the procedure, the best resulting solution is reported.

Fig. 2
figure 2

Structure of the proposed solution algorithm

Full size image

Initial solutions

The initial solutions are generated in two stages: the first stage constructs the routes and the second stage determines the value of other variables based on the fixed routes. Figure 3 illustrates the process of generating the initial solutions. To construct the routes, 75% of the blood centers are chosen randomly in each period. For example, in Fig. 3, blood centers 1, 3, 5, and 6 are chosen from the six blood centers. Then, the chosen blood centers are inserted in the route based on the cheapest insertion rule. Then, the two-opt operator is used to improve the constructed route as much as possible. For example, the constructed route in Fig. 2 is improved by displacing 1–5 and 3–6 arches. Finally, the mathematical model for the flow network problem is solved to determine the value of other decision variables. This is the second stage in generating the initial solutions. The network flow model is formulated as follows:

$$ \min {\text{PCO}} + \sum \limits_{s \in S} p_s Z^s + \lambda \sum \limits_{s \in S} p_s \left( {\left( {Z^s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} Z^{s^{\prime}} } \right) + 2\theta_s } \right) + \omega \sum \limits_{s \in S} \sum \limits_{i \in V} \sum \limits_{t \in T} \sum \limits_{r \in R} p_s \delta_i^{t,r,s} , $$
(53)
Fig. 3
figure 3

Generating initial solutions

Full size image

subject to constraints (7), (21), (40)–(52), and

$$ q_i^{t,g} \le C_i \left( {r_i^t } \right) \quad t \in T, g \in G, i \in I, $$
(54)

where \(r_i^t\) is a 0–1 parameter implying presence of blood center \(i\) in the vehicle route in period \(t\).

Local search

Two local search strategies are applied for improving the solutions. The first strategy, namely inserting the best solution, adds those non-visited blood centers to the route which result in more reduction in total costs. Based on Archetti et al. [29], the following model is formulated to identify these blood centers and the optimal supply quantity, quantity of products delivered by the vehicle, and the quantity of products transshipped among blood centers:

$$ \min \sum \limits_{i \in V} \sum \limits_{t \in T} b_i^t v_i^t + PCO + \sum \limits_{s \in S} p_s Z^s + \lambda \sum \limits_{s \in S} p_s \left( {\left( {Z^s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} Z^{s^{\prime}} } \right) + 2\theta_s } \right) + \omega \sum \limits_{s \in S} \sum \limits_{i \in V} \sum \limits_{t \in T} p_s \delta_i^{t,s} , $$
(55)

subject to constraints (7), (21), (40)–(52), and the followings:

$$ q_i^{t,g} \le C_i \left( {r_i^t - v_i^t } \right)\quad t \in T, g \in G, i \in I, s \in S $$
(56)
$$ v_i^t \le 1 - r_i^t \quad t \in T, i \in I $$
(57)
$$ v_i^t \in \left\{ {0,1} \right\} \quad t \in T, i \in I , $$
(58)

where \(v_i^t\) is a binary variable for presence of blood center \(i\) in the route at period \(t\). \(b_i^t\) denotes the cost of inserting blood center \(i\) at period \(t\) which is obtained from \(b_i^t = c_{ij} + c_{jk} + c_{ik} \) (see Fig. 4). Constraint (56) ensures that the product is only delivered to those blood centers that are present or inserted in the route. Constraint (57) implies that only those blood centers that are not present in the route can be inserted. The second local search aims at removing those visited blood centers that result in more reduction in the objective function compared to other blood centers. A mathematical model is used for determining such blood centers as well as the value of other decision variables:

$$ \min \sum \limits_{i \in V} \sum \limits_{t \in T} a_i^t u_i^t + PCO + \sum \limits_{s \in S} p_s Z^s + \lambda \sum \limits_{s \in S} p_s \left( {\left( {Z^s - \sum \limits_{s^{\prime} \in S} p_{s^{\prime}} Z^{s^{\prime}} } \right) + 2\theta_s } \right) + \omega \sum \limits_{s \in S} \sum \limits_{i \in V} \sum \limits_{t \in T} p_s \delta_i^{t,s} , $$
(59)
Fig. 4
figure 4

Procedure of calculating cost of inserting/removing a specific blood center into/from the route

Full size image

subject to constraints (7), (21), (40)–(52), and

$$ q_{i}^{{t,g}} \le C_{i} \left( {r_{i}^{t} - u_{i}^{t} } \right)\quad t \in T,~g \in G,~i \in I,~s \in S $$
(60)
$$ u_i^t \le r_i^t \quad t \in T, i \in I, $$
(61)
$$ u_i^t \in \left\{ {0,1} \right\} \quad t \in T, i \in I, $$
(62)

\(u_i^t\) is a 0–1 variable implying removal of blood center i from the route at period \(t\). \(a_i^t\) denotes the amount of transportation cost reduction due to removing blood center i from the route which is obtained from \( a_i^t = c_{ik} + c_{jk} - c_{ij}\) (see Fig. 4). Constraint (60) ensures that the products can be only delivered to those blood centers that are present in the route and are not removed. Constraint (68) implies that only the blood centers that are present in the route can be removed.

Shocking stage

Since there is the possibility of falling into local optima, the proposed algorithm includes a shocking stage to change current solutions and further search the solution space. To do so, a time period is chosen randomly and all visited blood centers in that period are removed. The previously-presented flow network model is used to determine the value of other decision variables [51].

Computational experiments

We carried out detailed computational analyses to further investigate the proposed model and evaluate performance of the proposed solution algorithm in terms of solution quality and time. First, the proposed solution algorithm was implemented in a set of Inventory-Routing Problems with Transshipment under the Maximum Level replenishment policy (IRPT-ML) [33]. Results were compared to those obtained from Adaptive Large Neighborhood Search (ALNS) algorithm, which was first proposed by Adulyasak et al. [25] for solving routing problems. In addition, Coelho et al. [33] and Civelek et al. [51] used this algorithm for solving the IRPT-ML problem. In the second stage, data from the real case of distributing blood platelets from Tehran Blood Transfusion Center, Iran to hospitals was used to examine the proposed model and solution algorithm. The algorithm is coded using MATLAB 2014. To obtain the exact solution of the robust model, the flow network model and the local search models were solved using GAMS v24.1 with CPLEX Solver v12.6. All computations were performed on a PC with Core I5, 2.6 GHz CPU and 6 GB RAM.

Application of the proposed solution algorithm on the IRPT-ML problem

The efficiency of the proposed solution algorithm was investigated in a set of IRPT-ML problems. Four clusters of numerical problems were considered with low and high inventory costs with 5–50 nodes over three time-periods and with 5–30 nodes over six time-periods. Five numerical examples were considered for each node and the first one was chosen to be solved using the proposed solution algorithm. The algorithm was implemented with a single product age, a single scenario, a sufficiently large shortage cost, relaxing constraints (40), (47), and (49), and converting constraint (7) to an equality.

The results obtained from the proposed solution algorithm were compared to those obtained from the ALNS algorithm, which was first proposed for solving routing problems. Adulyasak et al. [25] used this algorithm for solving green routing problem with simultaneous pickup and delivery and time window. In addition, Coelho et al. [33] used this algorithm for solving the IRPT-ML problem. Table 1 shows the required input data of all the parameters used in mathematical model. Due to confidentiality of the exact data and proper representation, the input parameter values are scaled down and considered in a uniformly distributed manner [52]. A comparison of the results obtained from the proposed solution algorithm and the ALNS algorithm is provided in Table 2. As it can be seen from the results, most solutions obtained from the ALNS algorithm were improved using the proposed algorithm. In addition, the solution time of the proposed algorithm is considerably less than that of the ALNS algorithm. We assumed the maximum number of seven iterations for solving the examples.

Table 1 Input data values of parameters
Full size table
Table 2 Comparison between results obtained from the proposed solution algorithm and the ALNS algorithm for the IRPT-ML problem
Full size table

Computational results and discussion

In this study, distribution of blood platelets from Tehran Blood Transfusion Center to 70 social security hospitals in Tehran, Iran was investigated using the proposed model. The cost parameters including transportation, holding, and supply costs were identified based on real data. In addition, the amount of demand, inventory level at the beginning of the first period, supply capacity, vehicle capacity, and inventory holding capacity at hospitals were identified based on the number of hospital beds as follows:

  • Vehicle capacity: 35% of sum of the number of beds in hospitals.

  • Capacity of each hospital: 70% of the number of hospital beds.

  • Demand of each hospital: 4% to 8% of the number of hospital beds in the first scenario; 10–14% in the second scenario; 16–20% in the third scenario. The exact values were chosen randomly within the mentioned intervals.

  • Supply capacity: 8–20% of sum of the total number of hospital beds. The exact percentage was chosen randomly within the mentioned interval.

  • Inventory level of blood center at the beginning of the first period: 5–10% of the total number of hospital beds. It was chosen randomly in a way that the inventory of the third age group was 8% of total inventory and the inventory of the second age group was 2% of total inventory.

  • Inventory level of each hospital at the beginning of the first time-period: 5–10% of the number of hospital beds. It was chosen randomly in a way that the inventory level of the third age group was 6% of total inventory and the inventory of the second age group was 4% of total inventory.

Seven clusters of numerical experiments including 10 to 70 hospitals in four time-periods and three scenarios of low, average, and high demand were designed. The shortage cost was set between 10,000 and 100,000 and the variance coefficient \(\lambda \) was set to one for all experiments. The results of solving the problems with 10, 30, and 70 hospitals using the proposed solution algorithm and the branch-and-cut search algorithm are given in Tables 3, 4, and 5. The average results for all examples are also provided in Table 6. The solution time for the branch-and-cut search algorithm in GAMS was limited to 1800s for the problems with 10, 20, and 30 hospitals, while it was limited to 3600 s for the problems having 40 or more hospitals. The fourth column of these tables shows the percentage of improvement in the upper limit with respect to the lower limit. In addition, the last column gives the percentage of improvement of the solutions obtained from the proposed algorithm with respect to the upper limit obtained from the branch-and-cut search algorithm. The results given in Table 6 show that the quality of solutions obtained from the branch-and-cut search algorithm is reduced by increasing the problem dimension. In contrast, the proposed solution algorithm provides better solutions than the branch-and-cut search algorithm.

Table 3 Results of the proposed solution algorithm and the branch-and-cut search algorithm for the presented case study with 10 hospitals
Full size table
Table 4 Results of the proposed solution algorithm and the branch-and-cut search algorithm for the presented case study with 30 hospitals
Full size table
Table 5 Results of the proposed solution algorithm and the branch-and-cut search algorithm for the presented case study with 70 hospitals using
Full size table
Table 6 Average results obtained for the presented case study using the proposed solution algorithm and the branch-and-cut search algorithm
Full size table

Balance between solution robustness and model robustness

Figure 5 shows the balance between solution robustness (total costs were given in Toman) and model robustness (shortage or unsatisfied demand) for the presented case study with 70 hospitals. As illustrated in Fig. 5, by increasing \(\omega \) from 10,000 to 100,000, the total cost was increased up to 818,699,891.5 and the shortage was reduced to zero. In practice, a risk-averse decision-maker chooses higher \(\omega \) values to prevent shortage, but a risk-taking decision-maker chooses lower \(\omega \) values to reduce total cost. In addition, by increasing \(\lambda \), the solution will be more sensitive to changes in the input data. In this study \(\lambda \) is set to one as proposed by Farahani et al. [53].

Fig. 5
figure 5

Balance between model robustness and solution robustness

Full size image

The results of sensitivity analysis of two cases of with and without transshipment

A comparative analysis of total costs, inventory levels and supply rates as well as shortage amounts, and the amounts of supply loss was done in the presented problem with 70 hospitals in two cases: (1) transshipment among hospitals is allowed and (2) transshipment among hospitals is not allowed. The shortage costs were assumed 70,000 Tomans. According to the obtained results, the costs of inventory and supply did not undergo significant changes in cases in which transshipment among blood centers was allowed compared to its counterpart. Nevertheless, as shown in Fig. 6a, b, a considerable decrease in total costs and number of shortages was observed in these two different scenarios, i.e., transshipment allowance and prohibition among blood centers. This is because there is no need for extra supply and holding extra inventory, where transshipment is allowed, so that both supply and inventory costs are both reduced. The results also show that the product loss was the same in two cases and did not exceed the maximum allowable loss. Moreover, the results of the study revealed that the transshipments did not reduce supply loss and that adopting the FF policy as well as the optimal supply policy for product consumption resulted in zero product loss in all periods, except for the first period, in both cases.

Fig. 6
figure 6

Comparison of a total costs, b number of shortages, c production cost, and d inventory levels with and without transshipment among hospitals

Full size image

In the case, where transshipment among hospitals is not allowed, the results show that the supply rate was more than that of the other case (Fig. 6c). This is because a blood center having extra products can give them to a blood center facing shortage which, in turn, results in reduced shortage costs. It also avoids the need for extra supply and reduces supply costs. As shown in Fig. 6d, the comparative analysis also shows higher inventory levels, where transshipment is not allowed, which, in turn, incurs higher inventory holding cost. These results confirm that transshipment among blood centers avoids additional costs and product shortages, which, in turn, reduces total costs of the system.

Overall, the numerical results confirm satisfactory performance of the proposed solution algorithm. The proposed algorithm is benefited from powerful operators and procedures for generating initial solution and improving solution quality that help produce superior results while significantly reducing the solution time. For instance, the solution time of the problem with 30 hospitals is 3.44 s in comparison with the time limit of 1800s in the branch-and-cut algorithm. This is also accompanied with 14.04% improvement in solution quality (Table 6). The improvement in obtained solutions is even more evident in large-sized problem. As an instance, the heuristic algorithm gives a solution with 75.52% improvement in 3.14 s for the problem having 70 hospitals. These results along with those obtained from comparative analyses help establish reliability and validity of the mathematical model and proposed solution algorithm.

Conclusions

In this research, a mathematical model was developed for the integrated inventory-routing of blood goods under demand uncertainty by taking into account the FF policy for consumption of goods. The transshipment among blood centers was incorporated to better handle uncertainty in customer demand. A heuristic algorithm was also proposed for solving the model. The proposed algorithm was tested on a set of numerical examples using the data from the literature showing that it is more efficient with respect to both time and quality of solutions. The proposed model and solution algorithm were also investigated in a case of blood distribution services in Tehran, Iran. A comparison was made between two situations of with and without product transshipment between blood centers and the overall performance of the network with respect to supply, shortage, inventory, and total costs was evaluated. The comparative results indicated that allowing transshipment among blood centers reduces supply quantity at the supplier, the amount of product shortage, and the inventory level. This, in turn, results in the reduced total costs. This research can be extended in several ways. The proposed model can be further investigated by incorporating transshipment of non-consumed products from the blood centers to the supplier. In addition, determining the optimal transshipment routes among the blood centers in each time-period and taking into account time window for delivery of blood products to the blood centers can be the subject of future works. The presented problem can be also investigated in a more complex supply chain network consisting of multiple suppliers with multiple blood products.