Designing efficient medical cold supply chain networks: A multi-objective P-hub location and routing model under perishability and capacity constraints

Abstract

This study proposes an optimization model to strategically locate blood hubs and route blood product flows. The model considers real-world factors, including depreciation costs at hubs, blood wastage, product decomposition at lab centers, and handling of multiple blood products. It employs both single and multiple allocation strategies to mitigate shortages and meet all demand. The proposed multi-objective model aims to: (i) minimize the total cost, such as transportation costs, depreciation costs, and hub installation costs, and (ii) minimize the maximum time that products stay in the network. Computational results show that NSGA-II produces diverse Pareto-efficient solutions and achieves small optimality gaps for small-sized instances when compared with the ε-constraint method, while remaining effective for larger instances where exact methods become computationally intractable. The ε-constraint method is used to derive Pareto optimal solutions, and NSGA-II is applied to solve the proposed model. Sensitivity analysis indicates that blood demand and depreciation costs are the most influential parameters, and that optimal solutions occur at moderate hub capacities. Managerial implications reveal that the optimal number and location of hubs are determined not solely by geography but also by demand distribution, donor availability, and product compatibility.

1. Introduction

According to the World Health Organization, approximately 117 million blood donations are collected each year globally, yet a substantial shortage persists in many countries. ¹ Timely and equitable access to safe blood is a persistent challenge in global healthcare.^2,3 The availability of blood and its components remains highly constrained and subject to uncertainty. ⁴ Blood is an indispensable medical resource that requires careful stewardship to ensure timely and equitable access for patients in need. ⁵

In low and middle-income nations, insufficient access to safe blood continues to be a leading contributor to preventable deaths affecting both healthcare systems. ⁶ Blood supplies are sourced from voluntary donations, which face significant challenges, including limited donor participation rates, logistical delays in screening and testing and the inherently limited shelf lives of blood products.^7,8 Given these challenges, optimizing the blood supply chain is paramount. Developing a resilient and responsive network for blood collection, processing, distribution, and transfusion is essential for safeguarding public health and improving population-level health outcomes.^3,9

The design and optimization of hub networks in blood transfusion supply chains must account for specific economic pressures and service requirements unique to healthcare logistics.^7,10 An important challenge in hub location problems (HLPs) is to design models that reduce costs while enhancing service efficiency.^11–13 This necessitates a comprehensive modeling approach that accounts for capital and operating expenses. ¹⁴ Capital expenses in this context might include investments in specialized storage facilities, vehicles equipped for handling blood products, and the infrastructure required to support a sterile and efficient network. ¹⁴ Operating costs typically encompass expenses related to personnel, such as medical technicians, as well as ongoing costs for fuel and vehicle maintenance. Importantly, the models should incorporate detailed time management strategies to handle the perishable nature of blood products, ensuring that the maximum time blood spends in transit is minimized. ¹⁵

While the transportation of blood products from donors to recipients may prima facie resemble the logistics networks of other perishable products, it involves unique complexities inherent to blood transfusion supply chains. ¹⁵ Managing the blood distribution network is a challenge that necessitates planning and optimization to guarantee the reliability and availability of blood units. ¹⁰ Although hub location is generally a strategic decision and routing is an operational task, in the context of blood supply chains these two decisions are inherently interconnected. The perishability of blood products and strict compatibility constraints mean that routing considerations such as travel time and allocation flexibility directly influence the suitability of hub locations. Conversely, the strategic placement of hubs shapes feasible and cost-effective routing options. Addressing both simultaneously allows the model to more accurately represent the real-world trade-offs between cost and service level.

Recent studies highlight the growing use of multi-objective optimization and metaheuristic methods in healthcare and cold supply chain applications, offering improved resilience and efficiency under uncertainty. ¹⁶ For instance, Altunoglu and Batur Sir ¹⁷ optimized fixed and mobile blood donation centers using the ε-constraint method to balance cost and timeliness. Tunga, Dhar ¹⁸ applied an evolutionary neural approach to manage uncertainty efficiently. Guo, Chen ¹⁹ developed a multi-period, multi-objective closed-loop blood supply chain model that accounts for disruptions and uncertainties, emphasizing the role of mitigation strategies and multi-period design in enhancing flexibility and resilience in healthcare logistics. Abdolazimi, Pishvaee ²⁰ designed a tri-objective blood supply chain model under pandemic conditions using a Benders decomposition–based heuristic, improving computational efficiency by over 80% while addressing uncertainty and perishability in COVID-19 disruptions. Mansur, Wangsa ²¹ formulated a multi-echelon mixed-integer model integrating facility location, inventory, and production decisions to minimize total cost and effectively manage outdated blood products through coordinated planning across supply chain levels. However, despite this progress, few studies comprehensively integrate biological perishability, capacity constraints, and routing within a unified hub location framework, particularly for blood transfusion systems.

Despite the critical importance of blood transfusion networks, there has been a notable absence of research addressing the p-hub center problem within this specific context.^22,23 Existing studies often focus on either facility location or routing but rarely integrate them while accounting for biological constraints, depreciation costs, and multi-product flows.^21,24–27 This lack of comprehensive, multi-objective optimization models for blood transfusion networks leaves a gap in both academic literature and practical applications. Furthermore, biological factors such as perishability, depreciation, and blood-type compatibility are often simplified or treated separately, limiting the applicability of such models to real healthcare systems. Capacity restrictions and allocation flexibility across hub facilities have also received limited attention, despite their critical role in maintaining service continuity and minimizing waste. Prior approaches have emphasized algorithmic efficiency, while the managerial interpretation of cost–time trade-offs has been less explored. These gaps highlight the need for an integrated multi-objective framework that addresses both operational and biological dimensions of blood supply systems, enabling more accurate and practical decision-support tools for healthcare planners. Accordingly, the objectives of this study are as follows:

• To develop a multi-objective p-hub center and routing model that integrates hub location and routing decisions under perishability and capacity constraints.

• To minimize both total operational costs and the maximum circulation time of blood products within the network.

• To apply the ε-constraint method and NSGA-II to generate Pareto-optimal solutions and provide managerial insights for enhancing the efficiency and resilience of healthcare blood supply systems.

To address this gap, this paper develops a capacitated p-hub center location and routing model tailored to blood supply chains. The model simultaneously minimizes total costs including transportation, depreciation, and hub installation and the maximum time blood products remain in the network. It incorporates perishability, product decomposition, wastage rates, and multi-product flows while using both single- and multi-allocation strategies to ensure demand fulfillment and minimize transit times. These elements are crucial as they address several complex issues previously unexplored, such as managing depreciation costs, blood wastage, and product decomposition in lab centers while also accommodating the diverse needs of components like donor centers, processing centers, and distribution centers.

The proposed model utilizes both single and multi-allocation methods to effectively manage blood product flows across the network, ensuring that all customer demands are met and minimizing the time blood products spend in transit. To solve the inherently NP-hard problem, we employ a meta-heuristic algorithm, the non-dominated genetic algorithm (NSGA-II), complemented by the ϵ-constraint method to achieve Pareto optimal solutions. The efficacy of the model is validated through several test problems, with results benchmarked against GAMS solutions, highlighting its practical applicability and robustness in handling the real-world complexities of blood transfusion logistics.

The proposed model contributes to the literature in a number of ways. First, it applies the HLP to the blood distribution network, addressing a complex and vital system. Second, the model accounts for the cost of depreciation related to blood products, which is critical to an efficient distribution network. Third, it considers the blood products in the testing and processing centers, which affect the availability and quality of blood products. Fourth, the model incorporates both single allocation and multi-allocation strategies in the network, allowing for a balanced trade-off concerning costs and service level. Fifth, the model is designed for multiple objectives and utilizes the ϵ-constraint method along with the NSGA-II to determine Pareto-optimal solutions. The main contributions of this research are as follows:

• Extends the p-hub center model with biological and logistical constraints unique to blood supply chains.

• Applies NSGA-II and the ε-constraint method in a novel way to solve a complex, real-world blood network problem.

• Addresses a gap in healthcare logistics by delivering adaptable solutions and a decision-support tool for blood transfusion networks.

The remainder of the paper is organized as follows. First, the relevant literature on multi-objective p-hub center problems is reviewed. Next, the problem description, model assumptions, and mathematical formulation are presented. The NSGA-II algorithm, including its techniques and parameter settings, is then explained. This is followed by the computational results, including data generation, the ε-constraint method, sensitivity analysis, and managerial implications. The paper then discusses the main findings and implications before concluding with potential directions for future research.

2. Literature review

2.1. Hub location problems

An important strategic decision in supply chain management is determining the optimal location of operational nodes.^28–33 This process involves selecting certain nodes to serve as hubs and allocating non-hub nodes to them.^34–36 Hub location problems (HLPs) arise in multi-node systems where goods, people, or data must be transported between locations.^37,38 In order to simplify the movement of goods, HLPs entail choosing locations for hubs and assigning non-hub nodes to these hubs.^38,39 The hub location–allocation problem is typically classified into three subcategories: p-hub median, p-hub center, and p-hub covering. ⁴⁰ When hubs number is predetermined, HLPs are commonly termed p-hub location problems. ⁴¹

Prior research has extensively explored HLPs for efficient network design across different industries.^35,37 Theoretical and practical approaches to HLPs have evolved, incorporating advanced mathematical models, simulation techniques, and, more recently, AI-driven algorithms to address the complexities of modern network designs. ⁴² For extensive surveys of HLPs, interested readers may refer to Basallo-Triana, Vidal-Holguín ⁴³ and Alumur, Campbell. ³⁷ While the application of hub location strategies has evolved considerably, their integration into healthcare logistics, especially the blood supply chain, remains relatively underexplored, presenting a worthwhile opportunity for research and development.^8,24,44,45

2.2. Hub location problems in medical cold supply chains

This section reviews key studies on HLPs within the blood transfusion network. We highlight the gaps and limitations of the existing models and methods and motivate the need for a multi-objective approach that considers both operational efficiency and service quality. The p-hub center problem is particularly relevant for distribution networks requiring rapid handling of time-sensitive or perishable items, such as critical medical supplies and express courier services.^12,41 Its operational principles align closely with the requirements of blood supply chain design. ⁴⁶ An important decision in designing the blood distribution network is determining where to place blood processing and distribution centers that collect donations from collection sites and deliver them to hospitals. ⁸ Strategically locating hubs and optimizing operational routes are interconnected decisions in blood supply chains. Integrating both decision levels allows the model to support cost-effective blood distribution and flexible responses to demand variability and operational uncertainties. Despite the critical importance of the hub location choices, the literature on HLPs in blood transfusion networks remains limited, addressing different aspects of the problem. ⁷

A number of research address hub location or allocation as the primary decision layer.^7,8,47 For example, Sharma, Ramkumar ⁸ developed a hybrid tool combining Tabu search heuristics and Bayesian belief networks to optimally establish emergency blood centers during and after disasters, ensuring quick response times for hospitals. Wemelsfelder, den Hertog ⁴⁷ outlined a mixed-integer linear programming model to determine the blood distribution centers locations, striving to maximize access to hospitals and lower transportation costs. Chaiwuttisak, Smith ⁷ introduced a novel location-allocation model for establishing economical blood donation and supply in Thailand, advancing from traditional approaches by incorporating budgetary constraints and a dual facility type system.

A growing body of research integrates hub location with other operational layers, such as routing and inventory control. For instance, Rostami, Sadjadi ⁴⁸ developed a heterogeneous fleet open vehicle routing problem model with dynamic hub location to optimize blood sample collection in clinical laboratory networks, demonstrating significant transportation cost reductions compared to classical OVRP variants. Gilani, Sahebi ⁴⁹ proposed a multi-objective robust possibilistic programming model for a sustainable-resilient blood supply chain network, integrating facility location, inventory, and routing decisions under uncertainty to minimize cost, environmental impact, and shortage risk. Alikhani, Dezfoulian ⁵⁰ designed a multi-objective metaheuristic framework combining NSGA-II and particle swarm optimization to optimize location, allocation, and routing in a perishable blood supply chain, achieving improved cost-efficiency and reduced wastage. Fu, Ma ³ developed a multi-objective optimization framework for home health care routing and scheduling with shared services, combining a mixed integer programming model and a problem-specific knowledge-based artificial bee colony algorithm to simultaneously minimize operational costs and service tardiness.

Kaya and Ozkok ⁵¹ presented an integrated approach to designing a blood distribution network that optimizes facility location, inventory levels, and routing decisions. They proposed an innovative approach where selected hospitals act as local blood distribution centers, reducing inventory levels at each hospital by pooling resources. Hosseini-Motlagh, Samani ⁵² introduced a dynamic optimization approach for the procurement of blood from various potential donors, adopting a sequential scheduling method to address uncertainty. Their approach dynamically updates strategies based on daily fluctuations in donor behavior and urban conditions, significantly improving blood collection efficiency in a blood network. Dehaghani, Nawaz ⁵³ introduced a multi-objective mathematical approach aimed at optimizing the blood distribution process, applying the m/m/m/k queue model to solve the location-allocation problem. Rekabi, Garjan ²³ discussed the complexities of blood distribution, highlighting the importance of strategic allocation and precise positioning of network elements such as donor centers and blood processing facilities, while also addressing the inherent complexities associated with the varying shelf lives of blood products.

Across these studies, two trends emerge. First, there is a movement toward integrated, multi-objective models that address multiple decision layers. Second, there is growing recognition of perishability-aware optimization, acknowledging that ignoring time-sensitive decay in blood products or other medical supplies can undermine both cost efficiency and service quality. As shown in Table 1, previous research has made substantial progress in integrating facility location, inventory, and routing decisions across various healthcare logistics contexts. Most studies either simplify biological constraints or address perishability indirectly through inventory control. Few explicitly combine hub-based location and routing decisions under perishability, depreciation, and capacity constraints within a unified multi-objective model. However, persistent gaps remain: very few models simultaneously address multi-product flows, multi-allocation strategies, and capacity constraints while capturing the temporal dynamics of blood deterioration. Moreover, many approaches rely on deterministic assumptions, limiting their robustness in real-world conditions marked by fluctuating demand and supply.

Table 1.

Comparative synthesis of prior work on blood supply chain optimization and hub-based designs.

Authors	Domain/setting	Decision layers covered	Key findings	Reported gaps/limitations
Altunoglu and Batur Sir 17	Blood donation center network	Location (fixed/mobile); limited routing	Demonstrates trade-offs between cost and timeliness in siting decisions	Routing/perishability not integrated; capacity treatment simplified
Tunga Tunga, Dhar 18	Blood supply management under uncertainty	Network design focus; limited routing	Improves responsiveness under uncertain demand/supply	Limited biological constraints; capacity/perishability simplified
Guo, Chen 19	Closed-loop blood SC (multi-period)	Configuration, disruption mitigation	Flexibility via multi-period design and disruption handling	Routing integration limited; perishability often abstracted
Abdolazimi, Pishvaee 20	Pandemic blood SC	Location–inventory under disruption	Efficient solutions during pandemic disruptions	Deterministic sub-structures; routing/perishability simplified
Mansur, Wangsa 21	Multi-echelon medical network	Location–inventory–production	Coordinated planning reduces obsolescence	No explicit routing integration; perishability modeled indirectly
Kaya and Ozkok 51	Blood distribution network	Location–inventory–routing	Pooling at selected hospitals reduces inventory	Biological decay simplified; hub interdependencies underexplored
Rostami, Sadjadi 48	Sample collection network	Dynamic hub + OVRP routing	Dynamic hubs reduce cost vs. classical OVRP	No multi-product perishability; limited capacity interplay
Alikhani, Dezfoulian 50	Perishable blood SC	Location–allocation–routing	Lower wastage via integrated decisions	Depreciation/time-in-network not explicit; limited p-hub design
Dehaghani, Nawaz 53	Blood distribution process	Location–allocation	Optimized location–allocation improves service	Inventory/routing decisions excluded
Rekabi, Garjan 23	Blood distribution systems	Allocation/network positioning	Emphasized strategic allocation and shelf-life complexity	Lacked a quantitative optimization framework
Hosseini-Motlagh, Samani 52	Blood donor procurement	Scheduling under uncertainty	Adaptive donor scheduling improves efficiency	Limited integration with network routing

To address these gaps, this study formulates a capacitated p-hub center problem for a blood transfusion network, which considers multiple types of blood products, multiple objectives, and multiple layers of facilities. We use the epsilon-constraint technique together with the NSGA-II metaheuristic to solve the model and generate Pareto-optimal solutions. We implement our model and methodology in a synthetically generated case study data focused on the design of a blood transfusion network, examining the outcomes and evaluating the sensitivity of the model’s parameters.

Our study advances this literature in four ways. First, it integrates p-hub center location and routing within a single formulation while explicitly modeling perishability, depreciation, product decomposition, and capacity limits, moving beyond models that treat facility location or routing in isolation. Second, it supports both single- and multi-allocation, showing when multi-allocation yields meaningful time reductions for modest cost increases—evidence conveyed through decision-ready Pareto frontiers. Third, we employ a hybrid ε-constraint NSGA-II approach with parameter tuning to efficiently approximate trade-offs for medium-to-large instances where exact methods struggle. Fourth, we translate results into managerial guidance (hub count/location, capacity right-sizing, and allocation policy) to connect algorithmic outcomes with real planning choices in medical cold chains. The following section offers a detailed description of our model.

3. Problem description

This section introduces a mixed-integer model for designing a multi-objective p-hub location problem in a blood transfusion network. The network layers consist of collection sites, where blood is donated, testing and processing centers, where blood is screened and prepared, distribution centers, where blood is stored and delivered, and demand centers, such as hospitals, research laboratories, and others, that use blood. It should be noticed that in this blood transfusion network, processing centers and distribution centers are hubs. Figure 1 illustrates the visual representation of this design. The blood transfusion network is modeled as a multi-product network. The model considers products such as red blood cells, whole blood, plasma and platelets which vary in perishability rates and have different uses.

Figure 1.

Schematic illustration of the blood transfusion network.

The blood collected from donors at different sites is thoroughly tested for diseases, such as HIV, HBV, HCV and syphilis. Then, the safe blood is separated into three products according to the demand requirements: plasma, platelets, and red blood cells. This process takes place at testing and processing centers, which are the Type I hubs. From a 450 ml sample of whole blood, platelets, plasma, and red blood cells are derived, with any remaining whole blood reserved for transfusions. These products are transported to central distribution centers (Type II hubs) and then distributed to hospitals as per their demand (demand centers).

This paper addresses a multi-objective problem involving making strategic decisions for placing p hubs in the network and allocating non-hub nodes to hubs with single-multi allocation and determining the amount of each blood product type to be moved between nodes within the blood transfusion network.

3.1. Explicit assumptions regarding depreciation and perishability

Before presenting the mathematical formulation of the proposed model, we first introduce some assumptions made to simplify the problem. These assumptions are based on the characteristics of the blood transfusion network and the objectives of the decision makers. The main assumptions include the following: A one-to-one allocation method is adopted to assign each origin node to only one Type I hub (testing and processing centers), while a multi allocation strategy is adopted to assign destination nodes to more than one Type II hub (blood distribution centers). The hubs are predefined in number, and each has a restricted capacity. The demand is deterministic and predefined and the wastage rate at Type I hubs is known. The hub network is fully interconnected, and each hub has a unique opening fixed cost. Transparency of the model relies on clearly defined depreciation, perishability, time horizon, units of measurement, and consistency assumptions.

- Perishability of blood products

• The model assumes that blood products are perishable and must be transferred and delivered within their valid shelf lives to ensure safety and efficacy.

• The perishability constraints are expressed as a maximum allowable circulation time ( ${\rho }^f$ ) (e.g., hours or days) for each blood product type (f).

• Time Horizon: The perishable window is measured from the moment of donation or processing at the origin until delivery at the hospital.

• Units of Measurement: All time-related parameters, including demand fulfillment times, transit times, and shelf lives, are consistently expressed in hours or days throughout the model.

• This consistency allows the model to accurately enforce temporal constraints and limit transit durations related to each blood product’s perishability.

- Depreciation of blood and costs

• Depreciation in this context refers to the cost associated with the aging of blood products during transit and storage, reflecting value deterioration and increased wastage risk.

• The depreciation cost rate ( $C_d$ ) (e.g., cost per kilometer or per hour) captures the incremental value loss due to blood’s age during transportation.

• Time Horizon for Depreciation: It is assumed that the cost accumulates proportionally over the transit and handling time for each blood unit, measured in hours or days.

• The model explicitly factors in the distance traveled and time spent in transit to estimate depreciation costs for each flow, ensuring these costs align with the perishability constraints.

• Units of Measurement: All cost calculations related to depreciation are expressed in monetary units per kilometer or per hour, maintaining consistency across the model.

- Consistency and application

• The time units for perishability ( ${\rho }^f$ ), transit times ( $T_{ik}$ ), and depreciation costs ( $C_d$ ) are uniform and carefully defined to avoid ambiguity.

• All temporal parameters are modeled in the same unit (either hours or days), and this unit is clearly stated in the model formulation.

• The perishability constraints enforce that the total transit and handling time for each unit does not exceed the maximum valid shelf lives, ensuring product safety.

• The depreciation costs are calculated based on these same time measures, providing a realistic valuation of the deterioration of blood products during transit.

3.2. Sets

i: Sets of origin nodes {i = 1, 2, . . ., I}

j: Sets of destination nodes {j = 1, 2, . . ., J}

k: Sets of Type I hubs {k = 1, 2, . . ., K}

l: Sets of Type II hubs {l = 1, 2, . . ., L}

f: Sets of different blood products {f =1, 2, . . ., F}

3.3. Indices and parameters

Table 2 provides and defines the indices, parameters and variables employed in the model.

Table 2.

Summary of parameters, variables, and symbols used in the mathematical model.

Type	Symbol	Description
Indices
​	i	Origin nodes (collection sites)
​	j	Destination nodes (demand centers)
Parameters
​	K	Number of Type I hubs
​	L	Number of Type II hubs.
​	f	Blood product type
​	f1	Whole blood
​	f2	Red blood cell
​	f3	Platelet
​	f4	Plasma
​	$D_j^f$	The demand of demand center j for blood product f.
​	${cl}_k$	Testing and processing center k capacity.
​	${cb}_l$	Blood distribution center l capacity.
​	$C_{ik}$	Transportation cost for each unit of flow over the distance separating nodes i and k.
​	$C_{kl}$	Transportation cost for each unit of flow over the distance separating nodes k and l
​	$C_{lj}$	Transportation cost for each unit of flow over the distance separating nodes l and j.
​	$F_k$	The fixed cost associated with establishing a hub at node k.
​	$F_l$	The fixed cost associated with establishing a hub at node l.
​	$T_{ik}$	The travel time from node i to node k.
​	$T_{kl}$	The travel time from node k to node l.
​	$T_{lj}$	The travel time from node l to node j.
​	α	Percentage of blood wastage at hubs.
​	$w_i$	Processing time at the collection center i.
​	$w_k^f$	Time of processing blood in the testing and processing center k for blood product f.
​	$w_l^f$	Time of blood storage in blood distribution center l for blood product f.
​	ß	The cost reduction factor applied to routes linking hubs.
​	$C_d$	Depreciation cost for each unit of flow per kilometer.
​	$d_{ik}$	The length of the route between nodes i and k.
​	$d_{kl}$	The length of the route between nodes k and l.
​	$d_{lj}$	The length of the route between nodes l and j.
​	M	A large number.
Decision variables
​	$z_{ik\hspace{0.17em}}^f$	The amount of blood product travelling from collection center i to hub k.
​	$z_{kl\hspace{0.17em}}^f$	The amount of blood product f travelling from hub k (processing center k) to hub $l$ (blood distribution center l).
​	$z_{lj\hspace{0.17em}}^f$	The amount of blood product $f$ travelling from hub $l$ (blood distribution center) to demand center $j$ (hospital j).
​	$x_k$	Equals 1 if node k hosts a Type I hub; 0 otherwise
​	$x_l$	Equals 1 if node $l$ hosts a Type II hub; 0 otherwise.
​	$x_{ik\hspace{0.17em}}$	Takes a value of 1 if node i is assigned to hub at node k; otherwise, it is 0.
​	$x_{kl}$	Takes a value of 1 if hub k is assigned to hub at node l; otherwise, it is 0.
​	$x_{lj}$	Takes a value of 1 if node j is assigned to hub at node l; otherwise, it is 0.
​	${\rho }^f$	Indicates the maximum allowable time blood product f can circulate in the network before it deteriorates.

3.4. Model formulation

The mathematical model incorporates two objective functions. The first objective function minimizes: (i) the costs of establishing all Type I hubs, (ii) the costs of establishing all Type II hubs, (iii) the transportation cost from non-hub nodes to Type I hubs, between all hubs, and from non-hub nodes to Type II hubs, and (iv) the total depreciation cost for each unit of flow per kilometers in the network. The second objective aims to minimize the maximum duration that blood products remain within the network before being utilized.

$\mathit{Min}\hspace{0.5em}\underset{kϵK}{}\hspace{0.17em}(F_k.x_{kk})+\underset{l\in L}{}\hspace{0.17em}(F_l.x_{ll})+\underset{i\in I}{}\underset{k\in K\hspace{0.17em}}{}(C_{ik}.z_{ik}^f)x_{ik}+\underset{f\in F}{}\underset{k\in K}{}\underset{l\in L}{}\hspace{0.17em}({\hspace{0.17em}\beta C}_{kl}.z_{kl}^f)x_{kl}+\underset{f\in F}{}\underset{l\in L\hspace{0.17em}}{}\underset{j\in J}{}\hspace{0.17em}(C_{lj}.z_{lj}^f)x_{lj\hspace{0.17em}}+\underset{i\in I}{}\underset{k\in K}{}\hspace{0.17em}(C_d.d_{ik}.z_{ik}^f)x_{ik}+\underset{f\in F}{}\underset{k\in K}{}\underset{l\in L}{}(C_d.{d_{kl}.z}_{kl}^f)x_{kl}+\underset{f\in F}{}\underset{l\in L\hspace{0.17em}}{}\underset{j\in J\hspace{0.17em}}{}(C_d.d_{lj}.z_{lj}^f)x_{lj\hspace{0.17em}}$

(1)

$\mathit{Min}\hspace{0.17em}\mathit{Max}\hspace{0.17em}\{(\hspace{0.17em}{w}_i+\left(T_{ik}.x_{ik}\right)+\left(T_{kl}{.x}_{kl}\right)+\left(T_{lj}.x_{lj}\right)+\left(w_k^f.x_{kk}\right)+\left(w_l^f.x_{ll}\right))\}$

(2)

Subject to

${\displaystyle \sum }_k\hspace{0.17em}{x}_{ik}=1\hspace{1em}\forall \hspace{0.17em}i\ne k$

(3)

${\displaystyle \sum }_l\hspace{0.17em}{x}_{lj}\ge 1\hspace{1em}\forall \hspace{0.17em}j\ne l$

(4)

${\displaystyle \sum }_k{x}_{kk}=K$

(5)

${\displaystyle \sum }_l{x}_{ll}=L$

(6)

$x_{lj}\le {\hspace{0.17em}x}_{ll}\hspace{1em}\forall \hspace{0.17em}l,j,j\ne l$

(7)

$x_{ik}\le x_{kk}\hspace{1em}\forall \hspace{0.17em}i,k,i\ne k$

(8)

$x_{kk}+x_{ll}-1\le x_{kl}$

(9)

$x_{kl}\le x_{kk}$

(10)

$x_{kl}\le x_{ll}$

(11)

$\left(1-\alpha \right){\displaystyle \sum }_{i\in I}{z}_{ik}\ge {\displaystyle \sum }_{lϵL}{z}_{kl}^f\hspace{1em}\forall f,k$

(12)

${\displaystyle \sum }_{k\in K}{z}_{kl}^f\ge {\displaystyle \sum }_{j\in J}{z}_{lj}^f\hspace{1em}\forall \hspace{0.17em}f,l$

(13)

${\displaystyle \sum }_{f\in F}{\displaystyle \sum }_{l\in L}{z}_{kl}^f\le {cl}_k.x_{kk}\hspace{1em}\forall \hspace{0.17em}k$

(14)

${\displaystyle \sum }_{f\in F}{\displaystyle \sum }_{j\in J}{z}_{lj}^f\le {cb}_l.x_{ll}\hspace{1em}\forall \hspace{0.17em}l$

(15)

${\displaystyle \sum }_l{Z}_{lj}^f\ge D_{j\hspace{0.17em}}^f\hspace{1em}\forall f,j$

(16)

$z_{ik}^f,z_{kl}^f,z_{lj}^f\ge 0\hspace{1em}\forall \hspace{0.17em}\mathrm{f},\mathrm{i},\mathrm{j},\mathrm{l},\mathrm{k}$

(17)

$x_{ik},x_{kl},x_{lj,}{x}_{kk},x_{ll}\in \left\{0,1\right\}$

(18)

The assignment of each origin node (collection site) to exactly one Type I hub is ensured by Equation (3) and constraint (18). Equation (4) incorporates a multi-allocation strategy that allows each destination node to be allocated to more than one Type II hub to ensure the network can meet the maximum possible demand. Ensuring the fulfillment of all potential demand within the blood supply chain is essential, as this has life-or-death implications.

Equations (5) and (6) indicate that the hub network can only establish $K$ (testing and processing centers) and $L$ (distribution centers) hubs. Constraint (7) and constraints (8) mandate that a non-hub node j must be assigned to an open hub at node l, while a non-hub node i must be assigned to an open hub at node k. Constraints (9) to (11) ensure the linearity of the product $x_{kk}$ × $x_{ll}$ . Equation (12) pertains to the decomposition of the blood at Type I hubs. Four primary components—whole blood, platelets, plasma, and red blood cells—are derived from blood, based on the demand for each product. The wastage rate in these equations is represented by the parameter α. Equation (13) outlines the flow equilibrium across testing and processing facilities, distribution centers, and demand points. Equations (14) and (15) impose capacity limitations on testing and processing facilities, blood distribution centers, and demand centers, respectively. The constraints in Equation (16) ensure that the demand is addressed and fulfilled by all blood distribution centers, while the Equation (17) and Equation (18) are positivity and binary constraints respectively.

Figure 2 presents a flowchart of our proposed model. This diagram provides an overview of how different components of the model interact to generate efficient and practical solutions for blood supply chain design. The problem addressed in our study combines capacitated p-hub location decisions with multi-product flows and time-sensitive routing under perishability and depreciation constraints, resulting in a large-scale mixed-integer program that is NP-hard. Even with linearization, the problem size grows rapidly with the number of origins, hubs, destinations, and products, leading to an exponential increase in potential hub configurations and allocation/routing combinations. In practice, exact solvers (e.g., GAMS with ε-constraint) solve small to medium instances within a reasonable time, but larger networks (tens of nodes with multiple products) become intractable. The NSGA-II heuristic, complemented by the ε-constraint method, provides scalable Pareto-frontier solutions, with reported results showing good diversity and near-optimality for medium-to-large instances. For real-time or very large-scale deployments, suggested strategies include problem-specific initialization to seed promising regions, parallelization, and hybridization with exact/decomposition approaches to manage combinatorial growth.

Figure 2.

Flowchart summarizing the methodological framework of the proposed multi-objective blood supply chain optimization model.

3.5. Linearization of the model

The efficiency of a proposed model can be significantly improved by presenting it in a linear form, as linear mathematical models are computed much faster than nonlinear ones. The proposed model is non-linear due to the multiplication of the real positive variables $z_{ik}^f,z_{kl}^f,z_{lj}^f$ with the binary variables $x_{ik},x_{kl},x_{lj}$ in the first objective function. To linearize these terms, we introduce new real positive variables, $R_{ik}{,R}_{kl}{,R}_{lj}$ corresponding to these multiplications. ⁵⁴

$R_{ik}=z_{ik}^f{\hspace{0.17em}*\hspace{0.17em}x}_{ik}$

(19)

$R_{kl}=z_{kl}^f\hspace{0.17em}{*\hspace{0.17em}x}_{kl}$

(20)

$R_{lj}=z_{lj}^f*x_{lj}$

(21)

Note the subsequent logical conditions:

${if\hspace{0.17em}x}_{ik}=0\to R_{ik}=0\hspace{1em};else\hspace{1em}{R}_{ik}=z_{ik}^f\hspace{1em}\forall i,k$

(22)

$\text{if}\hspace{0.17em}{x}_{kl}=0\to R_{kl}=0\hspace{1em};else\hspace{1em}{R}_{ik}=z_{kl}^f\hspace{1em}\forall i,k$

(23)

$\text{if}\hspace{0.17em}{x}_{lj}=0\to R_{lj}=0\hspace{1em};else\hspace{1em}{R}_{lj}=z_{lj}^f\hspace{1em}\forall i,k$

(24)

We can also add the following constraints for new variables $\left(R_{ik},{\hspace{0.17em}R}_{kl},{\hspace{0.17em}R}_{lj}\right)$ :

${\hspace{0.17em}R}_{ik}-M{x}_{ik}\le \mathbf{0}\hspace{0.17em}\forall i,k$

(25)

$-z_{ik}+R_{ik}\le \mathbf{0}\hspace{0.5em}\forall i,k$

(26)

${{\hspace{0.17em}z}_{ik}-R}_{ik}+M{x}_{ik}\le M\hspace{0.17em}\forall i,k$

(27)

$R_{kl}-M{x}_{kl}\le \mathbf{0}\hspace{0.17em}\forall i,k$

(28)

$-z_{kl}^f+R_{kl}\le \mathbf{0}\hspace{0.17em}\forall i,k$

(29)

${z_{kl}^f-R}_{kl}+M{x}_{kl}\le M\hspace{0.17em}\forall i,k$

(30)

$R_{lj}-M{x}_{lj}\le \mathbf{0}\hspace{0.17em}\forall i,k$

(31)

$-z_{lj}^f+R_{lj}\le \mathbf{0}\hspace{0.17em}\forall i,k$

(32)

${z_{lj}^f-R}_{lj}+M{x}_{lj}\le M\hspace{0.17em}\forall i,k$

(33)

Subsequently, the second objective function can be expressed in linear form as detailed below:

$\text{Minimize}\hspace{0.17em}{\hspace{0.17em}\rho }^f$

(34)

Then, the following constraint is introduced into the model:

$\left\{\left(w_i+\left(T_{ik}.x_{ik}\right)+\left(T_{kl}{.x}_{kl}\right)+\left(T_{lj}.x_{lj}\right)+\left(w_k^f.x_{kk}\right)+\left(w_l^f.x_{ll}\right)\right)\right\}{\le \rho }^f$

(35)

Constraint (35) underlines that the maximum duration of each blood product type inside the network does not surpass its associated expiration time $\left({\rho }^f\right)$ .

4. Non-dominated sorting genetic algorithm II (NSGA-II)

We employed the multi-objective algorithm (NSGA-II) to solve the p-hub center problem and generate approximate Pareto-efficient solutions. To illustrate the model’s applicability, several small test problems were solved and their solutions were compared with the GAMS results. The DICOPT solver was used to address the model through a real case study, showcasing the applicability and reliability of the modeling framework. The NSGA-II proposed by Deb, Pratap, ⁵⁴ prioritizes and determines population fronts using the non-dominance and crowding distance techniques. Through crossover and mutation, the algorithm develops new solutions that are merged with the existing population to form the next generation. ⁵⁵ This process follows the concepts of non-dominance and crowding distance. The following sections provide detailed explanations of the non-dominance approach, the calculation of crowding distance, and the crowding selection operator.

In multi-objective network design, hybrid strategies that combine evolutionary search with scalarization or decomposition often enhance convergence and preserve diversity. In this study, NSGA-II provides global exploration and diversity control (via non-dominated sorting and crowding distance), while the ε-constraint mechanism offers structured steering toward extreme trade-off regions. This complementary design helps cover the Pareto set more evenly than using either component alone. For completeness, we note alternative families such as MOEA/D (decomposition-based) and SPEA2 (strength-Pareto based),^56,57 which are well-known baselines for multi-objective optimization; our choice prioritizes transparent trade-off control (ε bounds) together with diverse frontier generation (NSGA-II), which is particularly useful when decision-makers need interpretable cost–time solutions.

The NSGA-II algorithm has demonstrated robust performance in solving complex multi-objective problems across various domains, including transportation and healthcare logistics.^3,33,58 However, its computational efficiency may decrease as the problem size increases, primarily due to the exponential growth of the search space associated with larger networks. ⁵⁹ Specifically, the number of potential hub locations, assignment combinations, and routing possibilities escalates combinatorically with the number of nodes, which can lead to increased computational time and memory requirements.

In our current study, the algorithm effectively generates diverse Pareto frontiers for medium-sized instances within reasonable computational times. For larger-scale problems, potential strategies to enhance scalability include: (i) implementing problem-specific heuristics to generate initial populations close to promising regions, (ii) applying parallel computing techniques to distribute computational load, and (iii) integrating hybrid methods that combine NSGA-II with exact algorithms or decomposition techniques to manage complexity.

While NSGA-II remains a powerful and flexible tool, its application to significantly larger instances should be undertaken with these strategies in mind. Future research can explore these avenues to extend the approach for large-scale, real-world blood supply networks, ensuring that solution quality and computational feasibility are maintained.

4.1. Solution encoding scheme

To apply NSGA-II to the proposed blood supply chain design problem, each chromosome is defined as a structured representation of the key strategic and allocation decisions in the network. As illustrated in Figure 3, the encoding is composed of four linked segments. The first segment represents the opening decisions for Type I hubs, where each gene is binary and indicates whether candidate hub k is opened. The second segment represents the opening decisions for Type II hubs, where each gene similarly indicates whether candidate hub l is opened. The third segment encodes the allocation of origin nodes to Type I hubs, so that each origin node i is assigned to one admissible opened Type I hub. The fourth segment encodes the allocation of destination nodes to Type II hubs, allowing each demand node j to be assigned to one or more admissible opened Type II hubs in accordance with the multi-allocation structure of the model.

Figure 3.

NSGA-II chromosome encoding and decoding process.

The flow variables are not encoded independently as a separate unrestricted chromosome segment. Instead, they are determined through a decoding procedure after the location and allocation decisions are specified. For each chromosome, the opened hubs and allocation pattern define the feasible network structure, and the corresponding blood-product flows are then derived subject to demand satisfaction, capacity limitations, perishability restrictions, and the model constraints. This representation reduces infeasible search regions and ensures that the chromosome directly reflects the core combinatorial decisions of the problem.

During initialization, chromosomes are generated randomly while respecting the required numbers of Type I and Type II hubs and the basic assignment logic of the model. During reproduction, crossover operators are applied to the encoded segments to combine parent hub-location and allocation patterns, while mutation modifies selected genes to explore alternative network structures. After each crossover or mutation step, a repair mechanism is applied when necessary to restore feasibility with respect to hub-count, assignment, and capacity-related rules. The decoded solution is then evaluated according to the two objective functions of the model.

4.2. Non-dominance technique

In the context of the non-dominance technique, assume there are k objective functions. A solution x¹ is said to dominate another solution x² if two conditions are met: first, x¹ is not worse than x² across all objective functions; and second, x¹ is strictly better than x² in at least one of the k objective functions. Solutions in the second front are only outperformed by those in the first front.

4.3. Crowding distance

The crowding distance metric reflects the concentration of solutions near a particular solution. It offers an approximation of the density of surrounding solutions, with a higher value, as determined by equation (36), being more desirable.

${CD}_i={\displaystyle \sum }_{j=1}^k\hspace{0.17em}\frac{f_{j,i+1\hspace{0.17em}}^p-{\hspace{0.17em}f}_{j,i-1}^p}{f_{j,total\hspace{0.17em}}^{p,\mathit{max}\hspace{0.17em}}-f_{j,total}^{p,\mathit{min}}}$

(36)

In this equation, k denotes the number of objective functions. The term $f_{j,i+1\hspace{0.17em}}^p$ corresponds to the value of jth objective function associated with the (i+1)th solution, and $f_{j,i-1\hspace{0.17em}}^p$ corresponds to the value of jth objective function of (i-1)th solution, following the population’s sorting by the jth objective function value. The maximum and minimum value for the jth objective function is represented by the $f_{j,total\hspace{0.17em}}^{p,\mathit{max}\hspace{0.17em}}$ and $f_{j,total\hspace{0.17em}}^{p,\mathit{min}\hspace{0.17em}}$ .

4.4. Tournament selection operator

The binary tournament selection approach was used to choose solutions for crossover and mutation. The process began by selecting two potential population solutions and then determining the superior one using a combination of non-dominated sorting and crowding distance. When two solutions were located on the same front, priority was given to the one with the greater crowding distance.

4.5. Parameters setting

This section focuses on fine-tuning the NSGA-II parameters using the Response Surface Methodology (RSM). RSM is a statistical technique designed to improve processes by analyzing the relationships between various factors and their impact on response variables through regression modeling. First, we identify the parameters with a significant statistical impact on the algorithm results. The study involves two problem sizes, Small-S and Large-L, with each parameter measured at low (L) and high (H) levels to select values that generate superior solutions. The low level (L) is coded as −1, while the high level (H) is coded as +1. The details of the coded variables are presented below:

${\mathit{X}}_{\mathit{i}}=\frac{{\mathit{r}}_{\mathit{i}}-\left(\frac{\mathit{h}+\mathit{l}}{\mathbf{2}}\right)}{\left(\frac{\mathit{h}-\mathit{l}}{\mathbf{2}}\right)}$

(37)

Where $X_i$ represents the coded variable and $r_i$ denotes the real variable. The high and low levels of the parameter are denoted by ℎ and l, respectively. In the next section, parameter values and their levels are explained. The parameter tuning process for NSGA-II was conducted using the RSM to identify the optimal combination of crossover rate, mutation probability, and population size. To validate the RSM outputs, each parameter configuration was tested on representative problem instances, and the resulting Pareto fronts were compared using two performance indicators: hypervolume (HV) and spread (Δ). The hypervolume metric evaluates the extent of the objective space covered by the obtained solutions, reflecting convergence and diversity, while the spread indicator measures the evenness of solution distribution along the Pareto front. Higher HV and lower Δ values indicate better solution quality. The RSM results were found to be consistent across multiple replications, confirming the robustness of the tuned parameters.

4.6. Parameter values and assumptions of NSGA-II

In this section, we outline the specific parameter values and underlying assumptions used in the implementation of NSGA-II. The crossover ratio was set to 0.75 and the mutation ratio to 0.3, as per the RSM method. In small-scale problems, the initial population was set to 200, whereas in large-scale problems, it was increased to 300. The NFC-based stopping criteria were established at 40,000 for smaller problems and 120,000 for larger ones. The population was randomly generated, with the crossover operator applied to arbitrarily chosen solutions through one-point, two-point, or weighted combination methods. The mutation process employed a swap operator on a solution at random. For clarity and ease of reference, the main algorithmic settings used in the computational experiments are summarized in Table 3.

Table 3.

Algorithmic settings used in the computational experiments.

Component	Setting/Value	Small instances	Large instances	Remarks
Exact solution method	GAMS + ε-constraint	Used	Not retained as final solution method	Benchmark for small problems
Heuristic solution method	NSGA-II	Used	Used	Main solver for medium/large problems
Crossover rate	0.75	0.75	0.75	RSM-tuned
Mutation rate	0.30	0.30	0.30	RSM-tuned
Initial population size	200/300	200	300	Larger population for larger instances
Stopping criterion	NFC	40,000	120,000	Scaled by problem size
Selection operator	Binary tournament	Same	Same	Based on non-dominance and crowding distance
Crossover operators	One-point/two-point/weighted	Same	Same	Applied during reproduction
Mutation operator	Swap mutation	Same	Same	Applied during reproduction
Parameter tuning method	RSM	Yes	Yes	Validated using HV and spread metrics

5. Results

5.1. Data generation

We generated parameters for the computational experiments using uniform distributions as detailed in Table 4, which shows varied demand across different products. In our analysis, we evaluated and compared the exact solutions with heuristic methods. We employed GAMS and the ε-constraint method for deriving exact solution and the NSGA-II algorithm for the heuristic solutions.

Table 4.

Distributions and units of parameters for computational results.

Parameters	Corresponding random distribution	Unit
$d_{ik},$ $d_{kl},$ $d_{lj}$	Uniform∼(15,20)	Kilometer
$C_{ik},$ $C_{kl},C_{lj}$	Uniform∼(200,300)	Cost unit
$F_k,F_l$	Uniform∼(1000, 2000)	Cost unit
$T_{ik},T_{kl}{,T}_{lj}$	Uniform∼(20,50)	Hours
$D_{\mathrm{j}}^{\mathrm{f}1}$	Uniform∼(5,16)	Whole blood unit
$D_{\mathrm{j}}^{\mathrm{f}2}$	Uniform∼(30,70)	Red blood cell unit
$D_{\mathrm{j}}^{\mathrm{f}3}$	Uniform∼(30,50)	Platelet unit
$D_{\mathrm{j}}^{\mathrm{f}4}$	Uniform∼(20,40)	Plasma unit
${cl}_{\mathrm{k},}{cb}_l$	Uniform∼(200,2000)	Unit
$w_i$	Uniform∼(1,10)	Hours
$w_k^{\hspace{0.17em}f}$	Uniform∼(10,25)	Hours
$w_l^f$	Uniform∼(120,240)	Hours
$C_d$	Uniform∼(20,35)	Cost unit
α	0.1	​
β	0.95	​

The computational study is based on synthetically generated test instances rather than an externally sourced real-world dataset. For each instance, the model parameters were sampled from the distributions reported in Table 4 in order to generate representative blood supply network scenarios with multiple products, capacities, transportation costs, and time parameters. To evaluate the proposed approach across different levels of complexity, the experiments were organized into two groups: small-sized instances, used for exact-versus-heuristic comparison, and larger-sized instances, used to assess the scalability and Pareto-search capability of the NSGA-II algorithm. The problem size is represented by the tuple $i/j/k/l$ , corresponding respectively to the numbers of origin nodes, destination nodes, Type I hubs, and Type II hubs.

5.2. ε-constraint method

We adopt the ε-constraint method with NSGA-II to address the dual objectives of cost minimization and reduced circulation time in blood supply chains, where trade-offs between these goals are critical for ensuring both efficiency and service quality. The rationale for choosing the ε-constraint method combined with NSGA-II in this study is based on several key considerations related to the nature of the problem, the solution requirements, and the advantages of these approaches. The blood supply network design involves conflicting objectives, minimizing costs and minimizing maximum circulation time that require exploring a trade-off space. The ε-constraint method is a rigorous scalarization technique that transforms a multi-objective problem into a series of single-objective problems by constraining one objective while optimizing the others. ⁶⁰ This allows for efficient exploration of Pareto fronts in a controlled manner. NSGA-II is a widely recognized multi-objective evolutionary algorithm known for its ability to produce a diverse set of Pareto optimal solutions across complex search spaces. ³ Its features, such as non-dominated sorting and crowding distance, help maintain solution diversity, which is crucial in decision-making scenarios involving multiple conflicting criteria like cost and time.

Exact methods (e.g., Multi-Objective Integer Programming solvers) may struggle with large-scale, nonlinear, or combinatorial problems due to computational complexity. Metaheuristic approaches like NSGA-II can efficiently generate high-quality solutions within reasonable computational times, especially suitable for large or complex instances typical of real-world transportation and healthcare network problems. The combined approach allows for flexible inclusion of additional constraints or objectives without significant restructuring. The ε-constraint method, coupled with NSGA-II, can be extended easily to handle different problem nuances or additional criteria, making it an adaptable framework. Both the ε-constraint method and NSGA-II have been successfully applied in similar complex logistics and healthcare problems.^59,61 Their robustness, ease of implementation, and ability to handle discrete and nonlinear problems make them attractive choices.

The model is expressed in terms of two objective functions, where Z₁(x) denotes the first objective function, Z₂(x) denotes the second objective function, and X denotes the feasible region defined by the model constraints. The original model can therefore be written as:

$\begin{array}{ll}\min \hspace{0.17em}\left\{Z_1\left(X\right),Z_2\left(X\right)\right\},& \mathrm{x}\in \mathrm{X}\end{array}$

(38)

To explicitly transform the original bi-objective model into a constrained single-objective formulation, the following steps are applied. First, the two objective functions are identified as Z₁(x), and Z₂(x), where Z₁(x) represents the primary cost-related objective and Z₂(x) represents the time-related objective. Second, Z₁(x) is prioritized as the objective to be optimized, because it captures the overall economic performance of the blood supply chain, including transportation, depreciation, and hub establishment costs. Third, Z₂(x) is converted into a constraint by imposing an upper bound ε on its allowable value. Fourth, the value of ε is varied iteratively over its feasible range in order to generate different trade-off solutions. In this way, the original bi-objective model is transformed into a sequence of single-objective optimization problems, each corresponding to a different admissible level of the time objective.

In the ε-constraint approach, one objective function is selected as the primary objective and the other objective function is converted into a bounded constraint. In this study, Z₁(x) is optimized while Z₂(x) is treated through an ε-bound. The resulting formulation is:

$\begin{array}{l}\min \hspace{0.17em}{\mathrm{Z}}_1\left(\mathrm{x}\right)\\ \text{subject}\hspace{0.17em}\text{to}\\ \begin{array}{l}{\mathrm{Z}}_2\left(\mathrm{x}\right)\le \mathrm{\epsilon }\\ \mathrm{x}\in \mathrm{X}\end{array}\end{array}$

(39)

Here, ε denotes the upper bound assigned to the second objective function. By varying the value of ε over its admissible range, different efficient solutions can be generated. However, the conventional ε-constraint method may produce weakly efficient solutions. To address this issue, the augmented ε-constraint method proposed by Mavrotas ⁶² is employed. In the augmented form, a non-negative slack variable is introduced for the constrained objective, and a small augmentation term is incorporated into the primary objective function. The augmented formulation is written as:

$\begin{array}{l}\min \hspace{0.17em}{\mathrm{Z}}_1\left(\mathrm{x}\right)-\delta \left(s_2/r_2\right)\\ \text{subject}\hspace{0.17em}\text{to}\\ \begin{array}{l}{\mathrm{Z}}_2\left(\mathrm{x}\right)+s_2=\mathrm{\epsilon }\\ \mathrm{x}\in \mathrm{X},s_2\in {\mathbb{R}}^{+}\end{array}\end{array}$

(40)

In Equation (40), s₂ denotes the non-negative slack variable associated with the ε-bound imposed on the second objective function, r₂ denotes the normalization range of the second objective function, and δ denotes a sufficiently small positive scalar introduced to avoid weakly efficient solutions without changing the dominance structure of the main objective. In this way, the model first minimizes the primary objective function and, among solutions with the same objective value, gives preference to the solution with smaller slack in the constrained objective.

In the present study, the first objective function is optimized while the second objective function is handled through the augmented ε-constraint mechanism. The value of ε assigned to the second objective function is varied iteratively across its feasible range, and each resulting single-objective model is solved to generate Pareto-efficient solutions. The collection of these solutions is then used to construct the Pareto frontier reported in the computational analysis.

5.3. Computational results

Table 5 presents a comparison between the exact and heuristic solutions. We employed GAMS and the ε-constraint method to obtain the exact solution, and the NSGA-ΙΙ algorithm for the heuristic solutions. The table reports results from the ε-constraint method alongside those from the NSGA-II algorithm for small-sized instances. We compared outcomes across four node sizes: 10, 12, 14, and 15. Our experiments were carried out on a Pentium Core i5, M460 CPU laptop running at 2.53 GHz with 16 GB of RAM. The table highlights the problem sizes, the computational results of the objective functions using both the ε-constraint method and the NSGA-II algorithm, and the gap identified in our experimental tests as well as the Pareto solution for the objective functions.

Table 5.

Comparative analysis of the ε-constraint method and the NSGA-II algorithm.

No. of nodes	Problem size i/j/k/l	Obj 1 (ε-constraint method)	Obj 2 (ε-constraint method)	Obj 1 (heuristic method)	Obj 2 (heuristic method)	Gap obj 1	Gap obj 2
10	3/2/2/3	27925	18609	29498	17204	0.056	0.07
12	3/4/2/3	26528	55251	121217	51066	0.04	0.07
14	5/3/2/4	729996	122148	746409	110226	0.02	0.09
15	6/3/2/4	53818	76872	53135	72264	0.012	0.06

More specifically, the small-sized instances consisted of four test cases with node sizes 10, 12, 14, and 15, corresponding to the configurations 3/2/2/3, 3/4/2/3, 5/3/2/4, and 6/3/2/4, respectively. These instances were used to compare the NSGA-II results with those obtained from the ε-constraint method implemented in GAMS. For larger-sized instances, three configurations were examined, namely 7/4/3/6, 10/6/6/8, and 25/7/5/3, corresponding to node sizes 20, 30, and 40. These larger instances were used to evaluate the robustness and scalability of the NSGA-II algorithm when exact optimization became computationally impractical.

The results show small gaps of less than 10%, supporting the use of heuristic method. The NSGA-ΙΙ algorithm proved effective for small-sized problems with minimal gaps. Table 6 highlights that GAMS could not solve larger-sized problems (i.e., 20, 30, and 40), suggesting that heuristic methods are preferable for these instances. Results for these larger-sized problems using NSGA-ΙΙ algorithm are demonstrated in Table 6.

Table 6.

Pareto solutions for node sizes 20, 30, and 40 using the NSGA-II algorithm.

Parameters		NSGA-II algorithm
No. of nodes	Problem Size i/j/k/l	Obj 1	Obj 2
20	7/4/3/6	62571	251604
​		63253	250270
​		64147	249558
​		65958	249398
​		67786	249046
​		69143	248222
30	10/6/6/8	167984	1439808
​		176604	1426800
​		222408	1424832
40	25/7/5/3	214917	5069220
​		217944	5060265
​		226107	5039470
​		259572	5027055

It should be noted that the results in this section are reported primarily in terms of Pareto-efficient objective values and exact-versus-heuristic comparisons, rather than percentage improvements relative to a separately defined baseline configuration. For small-sized instances, the performance of NSGA-II is evaluated against the ε-constraint method through the reported objective gaps in Table 5. For larger instances, where exact solution by GAMS was not attainable within reasonable computational effort, the NSGA-II algorithm is assessed through the structure and diversity of the Pareto solutions reported in Table 6 and Figures 4–6.

Figure 4.

Pareto frontier for problem size 7/4/3/6 (20-node configuration).

Figure 5.

Pareto frontier for problem size 10/6/6/8 (30-node configuration).

Figure 6.

Pareto frontier for problem size 25/7/5/3 (40-node configuration).

To further visualize the trade-offs between the two objective functions, total cost and the maximum time blood products remain in the network, we present the Pareto frontiers for selected large-scale problem instances solved using the NSGA-II algorithm. Figure 4 illustrates the Pareto frontier for the configuration 7/4/3/6 (20-node problem), showing a smooth trade-off pattern as solutions progress from lower cost to reduced time. Figure 5 presents the Pareto solutions for the 10/6/6/8 instance (30-node problem), where a steeper decline in the time objective is observed as cost increases, indicating diminishing returns in time efficiency at higher costs. Figure 6 depicts the results for the largest test case (25/7/5/3, 40-node problem), where the solutions reveal a broader spread along the Pareto front, demonstrating the scalability and robustness of the NSGA-II approach in handling complex configurations with multiple competing objectives.

The Pareto frontiers presented in Figures 4–6 illustrate the trade-off between total cost and circulation time in the blood supply network. Each point represents an operational scenario that balances efficiency and service quality. The findings show that modest increases in cost can reduce product circulation time, especially for networks with higher perishability sensitivity. This insight helps healthcare planners and policymakers determine the most appropriate operating point based on available budgets and service priorities. For example, in emergency conditions, choosing a slightly higher-cost solution can achieve faster delivery and lower spoilage risk, whereas in stable periods, a cost-efficient solution may be preferred. These trade-offs serve as a decision-support guide for planning hub capacities, route selection, and inventory coordination in real healthcare logistics.

5.4. Sensitivity analysis

We conducted a sensitivity analysis of the main parameters including the depreciation cost per unit of flow per kilometer, demand at the demand centers and hub capacities (testing processing center k capacity and blood distribution center l capacity). We use ANOVA to examine how the optimal solutions vary with different levels of some key parameters that affect the cost. Table 7 summarizes the parameter levels considered in the sensitivity analysis. The results were evaluated using the Taguchi design method facilitated by Minitab 16 software to assess the impact of these parameters on the Pareto frontier and identify most influential factors. It illustrates how these parameters affect the cost objective function, emphasizing that blood demand is the primary factor, as increased demand leads to higher costs. The ANOVA results confirmed that variations in blood demand and depreciation cost per unit of flow per kilometer had statistically significant effects on the total cost objective (p < 0.05), validating their dominant influence on the model’s outcomes. Furthermore, the analysis suggests that the optimal solution occurs at the minimal level of demand and depreciation cost, and at the medium level of hub capacity. From a managerial standpoint, these findings highlight that improving demand forecasting accuracy and reducing depreciation costs through efficient fleet and infrastructure management should be an important priority in planning and operating blood supply networks.

Table 7.

Parameter levels specified for the sensitivity analysis.

Parameters	Level 1	Level 2	Level 3
$D_{\mathrm{j}}^{\mathrm{f}1}$	Uniform∼(1,10)	Uniform∼(2,9)	Uniform∼(5,16)
$D_{\mathrm{j}}^{\mathrm{f}2}$	Uniform∼(20,50)	Uniform∼(30,60)	Uniform∼(30,70)
$D_{\mathrm{j}}^{\mathrm{f}3}$	Uniform∼(15,30)	Uniform∼(20,40)	Uniform∼(30,50)
$D_{\mathrm{j}}^{\mathrm{f}4}$	Uniform∼(10,20)	Uniform∼(15,30)	Uniform∼(20,40)
${cl}_{\begin{array}{l}\mathrm{k}\hspace{0.17em}\end{array}}$	Uniform∼(100,1000)	Uniform∼(200,2000)	Uniform∼(300,5000)
${cb}_l$	Uniform∼(200,300)	Uniform∼(300,500)	Uniform∼(500,1000)
$C_d$	Uniform∼(10,25)	Uniform∼(15,30)	Uniform∼(20,35)

The sensitivity analysis highlights that the impact on the cost objective function is predominantly influenced by fluctuations in blood demand. As illustrated in Figure 7, increased blood demand consistently correlates with heightened costs. Given that higher demand necessitates expanded transportation and handling, this in turn raises operational costs. Minimal levels of demand and depreciation cost alongside medium hub capacities yield the most cost-effective solutions. The depreciation cost per unit of flow per kilometer also significantly impacts overall costs. Lower depreciation costs correlate with reduced total expenses, emphasizing the importance of maintaining efficient and cost-effective transportation infrastructure. This finding underscores the need for investment in durable and reliable transportation assets that help minimize depreciation-related costs over time.

Figure 7.

Impact of key parameter adjustments (depreciation cost, demand at centers, hub capacity) on the objective function.

While increased hub capacities typically suggest a readiness to handle larger volumes, our results indicate an optimal middle ground that avoids the unnecessary overheads associated with maintaining underutilized capacity. This insight is important for operational planning, emphasizing the need for balanced investments in hub capacities tailored to predicted demand levels rather than maximal provisioning. The findings offer recommendations for decision-makers, stressing the necessity of investing in transportation infrastructure to reduce depreciation costs, accurately forecasting demand to avoid over- or under-capacity at hubs, and striking a balance between capacity and demand can lead to substantial cost savings. These strategies collectively can enhance the efficiency and reliability of blood supply chains, ensuring timely and cost-effective delivery of critical blood products.

6. Discussion and implications

Effective management of medical cold supply chains requires decision frameworks that balance cost efficiency with service reliability. ¹⁸ In practice, healthcare planners must determine how much additional investment is justified to improve responsiveness and reduce product deterioration. ²⁷ The results of this study provide guidance for achieving such balance, demonstrating how trade-offs between operational cost and time performance can be systematically evaluated. The insights derived from the model align with best practices in blood supply chain network design, emphasizing the need for integrated decision-making rather than isolated optimization of individual processes.^24,27

We interpret our results through the lens of cost–time trade-offs and allocation/capacity choices that affect perishability outcomes. The findings show that moderate increases in cost can secure meaningful reductions in time-in-network, especially when multi-allocation is enabled. This is consistent with previous research highlighting that multi-objective approaches enhance responsiveness and service quality in blood logistics while maintaining cost efficiency.^17–19,63 Our integrated location–routing formulation with explicit perishability and depreciation complements prior work that emphasizes either facility design or routing in isolation, and it provides decision-ready frontiers that make policy trade-offs transparent to planners. ⁶⁴

Our model provides important insights into managing the blood supply chain. The model assists decision makers in balancing cost and time objectives. We outline some counterintuitive results from our model that can inform decision makers. First, the optimal quantity of hubs is not necessarily equal to the quantity of blood types but depends on the demand distribution and the product compatibility. For instance, when the demand for a specific type of blood is very low, it may be more efficient to supply it from a hub that also serves other blood types, rather than establishing a dedicated hub. Second, the optimal hub location is not always close to the demand centers, it may also be influenced by the availability and cost of collection sites. For example, if a collection site is far away from the demand center with a large and low-cost blood supply, may justify placing a hub nearby and transporting the blood products to the demand center, rather than relying on a closer but more expensive or limited source.

Third, the allocation of non-hubs to hubs is not always best determined by geographical proximity alone, product compatibility and the demand patterns are contributing factors. For example, a non-hub with high demand for a certain blood type might be better served by a more distant hub that has a sufficient supply of that type than by a nearby hub with limited supply availability, thereby avoiding shortages and wastage. Fourth, our study demonstrated that strategic hub location decisions and operational routing considerations are interdependent in perishable medical supply chains. For blood products, routing constraints such as perishability times, compatibility, and delivery urgency can influence the optimal placement of hubs, while hub locations directly shape feasible routing plans. By addressing these two aspects together, the proposed model offers decision-makers a more holistic framework for balancing cost efficiency with service reliability. Lastly, our findings provide practical insights for strategic planning in blood supply chain management. The model supports decisions on the optimal number, location, and capacity of hubs to balance cost and service coverage while preventing redundancy or underutilization. It helps align inventory and replenishment decisions with time-sensitive delivery needs. The integrated optimization of hub placement and routing also enables better coordination among collection, processing, and distribution centers and assists planners in evaluating trade-offs between cost efficiency and product viability. These insights can inform data-driven policies and investment priorities for more resilient and responsive blood distribution systems. These findings help practitioners and decision makers to design and manage blood supply chains more effectively and efficiently by considering the complex interplay of factors that affect the optimal HLP.

This study has a number of theoretical implications for the literature on HLPs and blood distribution networks. First, it enhances existing blood distribution network models by incorporating depreciation cost, blood product decomposition, blood waste and multi-product consumption. These factors reflect the realistic characteristics of blood products and their handling processes. Second, it illustrates the applicability and effectiveness of a hybrid metaheuristic approach that combines a genetic algorithm with ε-constraint to tackle the mixed integer nonlinear programming problem. This algorithm can manage the complexity and nonlinearity of the model and find solutions that approximate the optimal within a feasible timeframe.

7. Conclusion

This paper presents a mixed-integer nonlinear programming formulation to address the hub location problem in a blood distribution network. Four levels are considered in our model: distribution centers, testing and processing centers, collection sites, and demand centers. We developed a method to identify central nodes, allocate peripheral nodes to these central nodes, and establish the flow of transfers within the blood network. The model is multi-objective, taking into account costs such as transportation, depreciation, and hub installation, while also aiming to minimize the maximum time products spend within the network in its second objective. To obtain a Pareto solution, the ϵ-constraint method is utilized. Since HLPs at large scales become NP-hard and time-consuming, NSGA-II is applied to solve the proposed model. Lastly, to demonstrate the problem’s practical relevance, multiple test scenarios are solved, and for small-sized problems, the solutions are evaluated against those from GAMS results.

Our study is not without limitations. Due to data confidentiality policies and restricted access to operational data from blood transfusion centers and hospitals, we relied on synthetically generated data for model evaluation and sensitivity analysis. While the proposed model and algorithms demonstrate promising performance under these controlled conditions, further research is needed to assess their effectiveness in real-world applications. Future work should focus on collaborating with healthcare organizations to access anonymized operational datasets and evaluate the model’s applicability and robustness in practice. Additionally, the current model assumes deterministic demand and supply conditions and does not incorporate uncertainty modeling in parameters such as demand variability, lead times, or transportation disruptions. Future studies could incorporate stochastic demand and supply uncertainty through robust or stochastic programming approaches, improving the model’s resilience under real-world variability.

The computational intensity of the NSGA-II algorithm may limit its applicability in real-time or resource-constrained settings. The model’s reliance on accurate data for costs, demand, and wastage rates poses challenges in environments where data collection is inconsistent or incomplete. To extend the relevance and applicability of this model, future studies could integrate additional operational problems such as incorporating vehicle and equipment maintenance costs and taxes as well as blood compatibility constraints to the blood transfusion network, which could further enhance the model’s utility. Finally, expanding the objective function to reduce the overall duration blood products remain in testing and processing centers for blood analysis and diagnosis of disease in the blood, could streamline operations and reduce the risk of spoilage, ensuring that blood products are used effectively and safely. Another limitation of this study is that it does not factor in the probabilistic variations in the availability and demand for blood. Blood supply and demand are subject to uncertainty and variability due to factors such as emergencies, disasters, epidemics, seasonal fluctuations, donor behavior, and limited shelf lives. Future research could extend the model to incorporate stochastic demand scenarios or perform robustness analysis to improve the model’s practical applicability.

This study provides avenues for further studies on the blood transfusion network optimization model. Future studies could enhance the model’s robustness by considering the dynamic aspects of the blood transfusion network, such as the time-varying blood supply and demand, periodic or random disruptions at blood facilities or transportation links, and the adaptive decision-making of the blood supply chain managers. A dynamic optimization model that incorporates these factors could capture the temporal and spatial variations of blood transfusion networks and provide optimal solutions under diverse scenarios. A worthwhile research direction that could extend the scope and complexity of the blood transfusion network optimization model is to consider the hierarchical and cooperative structure of the blood supply chains, such as regional and national blood centers, local and central blood distribution centers, and inter-organizational and intra-organizational coordination and collaboration. A hierarchical and cooperative optimization model that reflects different levels and roles of blood entities and their interactions and dependencies could enhance blood distribution systems.