1-Introduction

Blood is not a normal product. Blood supply from donors is usually irregular and demands for blood products are stochastic. Still, Supply and Demand Match (SDM) is not effectively solved in different studies. Blood products are perishable that complicates blood issues. Also, in a blood supply chain, all customers are trying to help people and do their best to supply the required blood. This means that customers help each other in blood supply chain and this issue is rarely observed in commercial products; therefore, blood is an unusual and special product [1].
Pierskalla [2] has mentioned various characteristics of designing blood bank. He mentions that, first of all, blood as a product has high level of perishability that is consisted of numerous components with different life spans. Second, since there is high variability in supplying this product, any planning in this regard is difficult. According to Katsaliaki [3] and Schreiber [4], only 5% of blood donors population, donate blood. This issue emphasizes this fact that blood supply can be very unstable. Third, as Belien and Pheresis [1] pointed out, demand for blood products in demand zones is highly random and at least probable. Also, in most of the cases, demands for red blood cells and other blood components are random variables. However, since red blood cells and other blood components are from a common sources (donor), there is high level of dependence in inventory and supply section [2]. Also, the processes of collecting blood samples in standard packages, their storage and distribution require high costs and this is dependent on the number and types of isolated components. Finally, these components have different life spans that complicate supply chain design process.
In a blood supply chain process, whole blood is collected from the donors. The samples are processed at regional laboratories and convert into different blood products. These blood products are transferred to the hospitals and the patients receive them. In a geographical region, a Regional Blood Center (RBC) is responsible to supply blood products for hospitals. Therefore, a schedule should be prepared to collect donors’ bloods a few months earlier. To collect blood from donors, vans and equipment for blood storage are used according to the schedule. The collected bloods are transferred to the processing center where they will be examined for blood group detection, virus, and disease. The obtained components in this stage are stored in blood banks and will be transferred to the hospitals by dedicated transportation system. Then, hospitals decide about how to use these products. If storage of blood samples exceeds from a certain limit, they cannot be used for the patient and must be destroyed. Destroying blood components and products should be specifically investigated. Some of blood components such as platelets can be directly received from donors. This process is called Pheresis. In this process, donor is connected to a machine that is continuously rotating and the blood enters this device. Here, platelet is separated from the blood and the rest returns to donor’s body. This process is more costly compared to platelet extraction from blood ([1], [3], [4]).
Regarding the proposed issue, review of literature is presented in two parts: the first part includes articles and works that have specifically investigated blood as a special products and the second part includes articles that have focused on issues related to designing supply chain, locating facilities, and the related issues.
Blood as a product and inventory issues
Shortage of blood products can cause high costs for the society because it increases death rates [3, 4]. Articles that are reviewed here are related to blood banks where bloods are collected, examined, and identified by blood groups. Blood banks can be organized in different forms that most of them are between centralization and decentralization. This part of review, includes the following products: whole blood, red blood cells, platelets, and blood plasma.
There is not much difference between whole blood and red blood cells and most of the researchers, use “blood” for red blood cells. However, sometimes these two are confused, because red blood cells are the most demanded blood component. However, in those studies that have concentrated on donors, red blood cells were specifically the product of whole blood, because all bloods that are taken from donors, are a whole blood that are divided into other components. Therefore, when the writer does not point to the differences between blood groups, we suppose that he is talking about red blood cells [5]. According to the historical point of view, whole blood and red blood cells have the largest share in demands [6]. Few studies have been conducted on plasma products. Only several studies are concentrated on this issue that even did not exactly concentrate on plasma supply chain that is probably due to this fact that working with plasma availability has less complexities compared to other blood components, because this product has long life span that will not be perished soon [7]. Indeed, plasma is the only non-perishable blood component [6]. Frozen storage of plasma is a suitable way for storage, but since 20 minutes is needed for its thawing that is not normal in terms of emergency, a non-frozen accumulation for 5 days is maintained [8]. Usually, plasma is maintained in frozen state; therefore, those studies that deal with plasma, automatically deal with frozen blood.
Compared to red blood cells and whole blood, fewer studies have concentrated on platelets. As pointed out, the problem of this product is that compared to other blood products, it has smaller life span. This issue causes more complexity and according to the studies on this context, more works are required. Van Dijk et al. (2009) proposed a solution that filled the existing gap in literature. Also, Blake [10] considered their work as an invaluable attempt to solve platelet availability problem management, but rejects the findings of their study. Van Dijk formulized perishability problem by a dynamic planning technique. The problem was that problem solving at real size with dynamic planning is very difficult. They solved this problem by suggesting multiple processes as well as dynamic planning. This solution could be a suitable innovative solution to solve PLT problem, but Blake [10] did not agree with research hypotheses. Van Dijk supposed that it is possible to ignore life span distribution of stored samples while Blake proved that this cannot be ignored. Blake also pointed that the suggested solution for this problem cannot be implemented easily; therefore, he found that their answer in one case is an unjustified example.
Blood supply chain network
Regardless of numerous efforts for locating and designing supply chain problems, blood supply chain designing problem is rarely studied. A recent review that was conducted by Belien and Fores [1], states that regarding strategic decisions about locating facilities, there are many gaps in literature. Pierskalla [2] provided a summary of presented models to allocate donation areas and blood transfusion centers to blood bank centers that aimed to determine number of blood bank centers in each area, locating these facilities, and supply and demand match.
Daskin et al. [11] and Shen et al. [12] presented non-linear integer planning models for a single period locating-inventory issue in blood supply for hospitals. They suggested a comprehensive method to determine numbers and locations of distribution centers and level of inventory for storage in each center. Cetin and Sarul [13] presented a multi-objective model to determine locations of blood banks in hospitals or clinics using a goal programming method. Their model minimizes the expenses of blood bank establishment and the distances between blood banks and hospitals. Nagurney et al. [14] and Nagurney and Masoumi [15] presented models for regional systems network optimization that include blood collection sites, processing facilities, storage facilities, distribution centers, and hospitals. The presented models considered differences such as determination of optimized allocations and waste destruction costs in blood supply chain. Two applied studies on blood supply chain were conducted by Shahin et al. [16] and Sha and Huang [17]. Shahin et al. [16] presented a single period locating-inventory issue model for regional blood services in Turkish Red Crescent. Their issue was considered as a hierarchical issue that included regional blood centers (level one), blood stations (level two), and mobile units (level three). Sha and Huang presented a blood timing issue that focuses on emergency circumstances after earthquake in Beijing. The purpose of this model is to determine location and optimized allocation of blood temporary facilities that minimize operational costs within a certain period of time.
Katsaliaki [3] evaluated the performance of blood supply chain in Britain. His study aimed to discover those policies that lead to more economical production of supply chain that at the same time increases security level. Fontaine et al. [18] conducted a study to improve blood platelets supply chain by emphasis on the importance of close collaboration between blood centers and hospitals. They suggested an interdisciplinary method that was working with supply chain management experts. Rytilä and Spens [5] presented a simulation model from blood supply system. Simulation steps and the results were described in their article. Jabbarzadeh et al. [19] presented a robust model to design blood supply chain during and after natural disasters. A practical optimization model was developed to locate blood centers and decision making regarding allocation for the period after disasters. They provided an analysis to compare robust answers and model through a practical way.
According to the conducted reviews, we find that blood chain can be considered as a system, not only a subsystem of a larger system, and is true for most of supply chain systems. Also, blood supply chain is more attractive in practical terms and theoretical studies can be taken from real events according to perishability problem and disease experiments. Therefore, this study considers the problem of locating facilities and network designing in a uniform way. In fact, this study has presented an Uncertain Capacitated Facility Location-Network Design Problem (UCFLNDP). Since the problems related to blood bank designing are solved for different purposes (such as minimizing the distance and maximizing access), it is better to consider several objective functions at the same time. In this study, the correspondence of blood groups was considered as a binary matrix in model that provides blood products in demand zones using other blood groups. Therefore, the model has the capability for medical services network modelling in blood supply chain under uncertainty of blood groups’ demands that are not available. The supply chain network in this study is consisted of four levels including donation centers or blood collection sites, processing and testing laboratories, blood bank centers or distribution centers, and demand zones that can be hospitals, clinics or research laboratories. Locations of demand zones and laboratories were considered fixed and other facilities are located within the network.
According to the explanations, this study is organized as follow. In sections (2) and (3), the problem is stated and the suggested problem will be formulized using mathematical modelling. In section (4), the model will become linear, then is modelled in a robust state and will become uni-objective using limited Epsilon method. In section (5), the presented model is solved using imperialist competitive algorithm and multi-objective non-dominated genetic sorting algorithm. In sections (6) and (7), sensitivity analysis has been performed on model’s parameters and the management results are presented. Finally, the results are concluded.

2-Problem Definitions

The given network in this study is a multi-product and four-level network where these levels are consisted of blood donation areas, laboratories, blood bank centers, and demand zones. A representation of this network is shown in Figure (1). Also, there is a network in the last layer of the main network that connects some of demand zones to each other. Indeed, consideration of this sub-network is one of the main differences of blood supply chain design with other network designs. By considering this sub-network, hospitals can send their products to each other when there are demands and the blood bank is not capable of providing the required amounts. This seems logical when long hours are needed for blood banks to be capable of meeting the demands; therefore, other hospitals take necessary measures to cover the demands. Since, demand zone is the last part of this network and all products finally reach to this point, therefore, there is not any need for other networks in layers and of course, several hospitals can be advantage points in considering this sub-network at the last layer. Accordingly, this sub-network can be activated earlier, because hospitals are more close to each other compared to blood banks that are scattered. Regarding the performance of this network, blood products are transferred to laboratories and are tested to identify disease and blood groups. Then, these products are transferred to demand zones. These zones include hospitals, research laboratories, clinics, and emergency stations.

Fig 1. Schematic view of blood supply chain network

Research assumptions
The capacity of different suppliers of blood products is determined.
Fixed and operational costs of donation centers are determined.
Fixed and operational costs of blood bank centers are determined.
Transportation costs for every unit of blood product are determined at all network levels.
The interest rate is determined.
The waste rate in laboratory is determined. The potential zones for the establishment of donation centers are determined.
The potential zones for the establishment of blood bank centers are determined.

3-Modeling issues

In this section, mixed integer two-objective linear programming model is studied where locating donation centers and blood bank canters as well as allocation of these two to the donation centers and the amount of transferred products in network are considered. Here, the correspondence of blood groups has been considered. In the following, we will introduce the model and notations.
Notations
Sets and indices
I: candidate points of donation centers i∈{1.2.….I}
J: laboratories j∈{1.2.….J}
K: blood banks k∈{1.2.….K}
L: demand zones l.m∈{1.2.….L}
F: different blood products f∈{1.2.….F}
G: different blood groups g.p∈{1.2.….G}

Parameters:
〖dc〗_( i )^g: capacity of donation center (i) for blood group (g)
〖lc〗_j: capacity of laboratory (j)
〖cbb〗_k: capacity of blood bank (k)
〖hc〗_l: capacity of demand center (i) (e.g. hospital capacity)
d_l^fg: demands of center (i) for blood products (f) and (g)
T_max^f: maximum time that blood product (f) can be used
t_(ij ): travel time on the bow (I, j)
t_(jk ): travel time on the bow (j, k)
t_(kl ): travel time on the bow (k, l)
t_lm: travel time on the bow (l, m)
t_i: blood processing time at donation center (i)
t_( j)^(f ): blood processing time at laboratory (j)
t_( k)^f: storage time at blood bank (k)
t_( l)^f: storage time at hospital (l) for blood product (f)
g_i: fixed costs for the opening of donation center (i)
g_k: fixed costs for the opening of donation center (k)
f_i: operational costs of center (i)
f_k: operational costs of blood bank (k)
c_ij: travel time cost on the bow (I, j) for each blood product unit
c_jk: travel time cost on the bow (j, k) for each blood product unit
c_kl: travel time cost on the bow (k, l) for each blood product unit
c_lm: travel time cost on the bow (I, m) for each blood product unit
p_ij: fixed cost of using bow (I, j)
p_jk: fixed cost of using bow (j, k)
p_kl: fixed cost of using bow (k, l)
p_lm: fixed cost of using bow (I, m)
h_ij: operational cost of using bow (I, j)
h_jk: operational cost of using bow (j, k)
h_kl: operational cost of using bow (k, l)
h_lm: operational cost of using bow (I, m)
M: a large number
λ: interest rate
α: percentage of donated blood that is wasted
〖BC〗^fgp: g and p blood groups correspondence matrix for blood product (f)
Decision Variables
y_i: binary variable that gets one if donation center (i) opens, otherwise, it gets zero.
〖y’〗_k: binary variable that gets one if blood bank (k) opens, otherwise, it gets zero.
x_( ij)^g: the amount of blood group (g) that is transferred from donation center (i) to the laboratory (j)
x_( jk)^fg: the amount of blood group (f) that is transferred from laboratory (j) to the blood bank (k)
x_( kl)^fg: the amount of blood group (f) that is transferred from blood bank (k) to the hospital (i)
x_( lm)^fg: the amount of blood group (f) that is transferred from hospital (i) to the hospital (m)
〖x’〗_kl^( fgp): the amount of blood group (f) and blood group (g) that is transferred from blood bank (k) to the hospital (i) to provide blood group (p) demand.
〖x’〗_lm^( fgp): the amount of blood group (f) and blood group (g) that is transferred from blood bank (i) to the hospital (m) to provide blood group (p) demand.
〖slack〗_l^( fg): the amount of demands for blood group (f) and blood group (g) that are provided in the hospital (i)
z_ij: the binary amount that gets one if bow (I, j) is used, otherwise, it gets zero.
z_jk: the binary amount that gets one if bow ( j, k) is used, otherwise, it gets zero.
z_kl: the binary amount that gets one if bow (k, l) is used, otherwise, it gets zero.
z_lm: the binary amount that gets one if bow (l, m) is used, otherwise, it gets zero.
t^f: the maximum time the blood product (f) can be used before corruption.
Mathematical Model
(1) Min TC=∑_(i∈I)▒〖(g_(i ) 〗+f_i/λ) y_i+∑_(k∈K)▒〖(〖g^’〗_(K ) 〗+〖f^’〗_K/λ) 〖y^’〗_K+∑_(j∈J)▒〖(ρ_(ij ) 〗+h_ij/λ) z_ij+∑_(j∈J)▒〖∑_(k∈K)▒〖(P_(jk ) 〗+h_jk/λ) z_jk+〗 ∑_(k∈K)▒〖∑_(l∈L)▒〖(P_(Kl ) 〗+h_Kl/λ) z_Kl+〗
∑_(l∈L)▒〖∑_(m∈L)▒〖(ρ_(lm ) 〗+h_lm/λ) z_lm 〗+ ∑_(g∈G)▒∑_(i∈I)▒∑_(j∈J)▒〖c_(ij ) X_( ij)^g+∑_(g∈G)▒∑_(f∈F)▒∑_(j∈J)▒∑_(k∈K)▒〖c_(jk ) X_( jk)^fg+∑_(g∈G)▒∑_(f∈F)▒∑_(k∈K)▒∑_(l∈L)▒〖c_(kl ) X_( kl)^fg 〗〗〗+ ∑_(g∈G)▒∑_(f∈F)▒∑_(l∈L)▒∑_(m∈L)▒〖c_(lm ) X_( lm)^fg 〗+M×∑_f▒∑_g▒∑_l▒〖slack〗_l^( fg)
(2) Min T=∑_f▒∑_i▒∑_j▒∑_k▒∑_l▒∑_m▒〖(t_ij 〗 z_ij+t_jk z_jk+t_kl z_kl+t_i y_i+t_j^f+t_k^f y_k+(t_lm+t_l^(f ))×z_(lm))
S.t.
(3) ∀f∈F-{4}.∀g∈G.∀j∈J-G (1-a) ∑_(i∈I)▒X_ij^g ≥∑_█(k@∈K)▒x_jk^4g +max⁡〖{x_jk^fg |f≤3}〗
(4) ∀f∈F.∀g∈G.∀k∈K ∑_(j∈J)▒X_jk^fg ≥∑_(l∈L)▒x_kl^fg
(5) ∀f∈F.∀g∈G.∀l∈L ∑_(k∈K)▒X_KL^fg ≥∑_(m∈L)▒x_lm^fg
(6) ∀i∈I .∀g∈G
∑_(j∈J)▒X_ij^g ≤〖dc〗_i^g×y_(i )
(7) ∀j∈J ∑_(g∈G)▒∑_(f∈F)▒∑_(k∈K)▒〖x_jk^( fg)≤〖lc〗_(j ) 〗
(8) ∀k∈K ∑_(g∈G)▒∑_(f∈F)▒∑_(l∈L)▒〖x_kl^( fg)≤〖cbb〗_k× 〖y’〗_k 〗
(9) ∀l∈L ∑_(g∈G)▒∑_(f∈F)▒∑_(k∈K)▒x_kl^( fg) +∑_(g∈G)▒∑_(f∈F)▒∑_(m∈L)▒x_Ml^( fg) ≤〖hc〗_(l )
(10) ∀f∈F.∀p∈G.∀l∈L ∑_g▒〖〖{BC〗^fgp [∑_(k∈K)▒〖〖x’〗_kl^fgP+∑_(m∈L)▒〖(〖x’〗_ml^fgq 〗-〗〗 〖x’〗_km^fgq)]}+〖slack〗_l^fg=d_l^fg
(11) ∀f∈F.∀p∈G .∀k∈K.∀l∈L x_kl^fg≥∑_(p∈G)▒〖〖X^’〗_kl^fgp 〗
(12) ∀f∈F.∀p∈G.∀l∈L.∀m∈L x_lm^fg≥∑_(p∈G)▒〖〖X^’〗_lm^fgp 〗
(13) ∀f∈F.∀i∈I .∀j∈J.∀k∈K.∀l∈L.∀m∈L t_ij z_ij+t_jk z_jk+t_kl z_kl+t_i y_i+t_j^f+t_k^f y_k+(t_lm+t_l^f )×z_lm≤T^f
(14) ∀f∈F T^f≤T_max^f
(15) ∀g∈G.∀i∈I.∀j∈J
X_ij^g≤〖MY〗_i
(16) ∀f∈F.∀g∈G.∀j∈J.∀k∈K X_jk^fg≤〖MY^’〗_k
(17) ∀f∈F.∀g∈G.∀k∈K.∀l∈L X_kl^fg≤〖MY^’〗_k
(18) ∀g∈G.∀i∈I.∀j∈J
X_ij^g≤〖Mz〗_(ij )
(19) ∀f∈F.∀g∈G.∀j∈j.∀k∈K X_jk^fg≤〖Mz〗_jk
(20) ∀f∈F.∀g∈G.∀k∈K.∀l∈L
X_kl^fg≤〖Mz〗_kl
(21) ∀f∈F.∀g∈G.∀l∈L.∀m∈L X_lm^fg≤〖Mz〗_lm
(22) ∀f∈F.∀g∈G.∀i∈I.∀j∈J.∀k∈K.∀l∈L.∀m∈L X_ij^g .X_jk^fg .X_kl^fg . X_lm^fg ≥0
(23) ∀i∈I.∀j∈J.∀k∈K.∀l∈L.∀m∈L Z_ij .Z_jk .Z_kl .Z_lm .y_i .y_k∈{0.1}
The first objective function is responsible to minimize total costs (variable and fixed) of network utilization. The components of this objective function include establishment costs of donation centers and blood banks, utilization cost of donation centers and blood banks, fixed and variable costs, and variable of using each route and initial costs to establish a bow within the network. The second objective function minimizes the time period that each blood product remains in the network. The third limitation creates a condition in which total blood units taken from donors convert into four blood products in laboratory. Some of these blood units remain intact and the rest convert into three different products including platelet, plasma, and red blood cells. This analysis is conducted based on maximum level of demands for these three products. In this relationship, α represents blood wastes in laboratory that occurs due to security issues such as blood diseases. Limitations (4) and (5) are related to the relationship between laboratories, blood bank centers, and hospitals. Relationships (6) to (9) are related to limitations of donation centers, laboratories, blood bank centers, and hospitals. Limitation (10) guarantees that demands of all hospitals are considered by blood banks or other hospitals. Also, blood group match matrix is considered for each product in these limitations that enables the model to use other blood groups to provide demands of a specific blood group. This can be practical if that specific group is not available when matching condition exists. Limitations (11) and (12) are logical relationships between variables of blood group match. Limitations (13) and (14) are related to corruption characteristics of blood products that guarantee total time that each blood product remains in the system does not exceed the expiration date. Limitations (15) and (16) state that no product should pass those donation centers that are not established. Limitations (18) to (21) guarantee that blood transportation between network layers should be done if there is a link between them. It should be noted that limitations (18) to (21) guarantee that limitation (21-3) is for links between demand zones. Finally, limitations (22) and (23) are related to positive nature and binarity of variables.
Linearization
The suggested model in previous section, due to equation (3), is a nonlinear model. This not only influences time, but the quality of answers is included, as well. This relationship is changed into a linear relationship by substituting X_( ij)^”fg with max⁡{X_jk^fg |f≤3}. Therefore, Limitation (3) is removed and limitations (24) and (25) are added to the model:
(24) ∀f∈F-{4}.∀g∈G.∀j∈J (1-a) ∑_(i∈I)▒X_ij^g ≥∑_(k∈K)▒x_jk^4g +X_( ij)^”fg
(25) ∀f∈F-{4}.∀g∈G.∀j∈J X_( ij)^”fg≥X_jk^fg

4-Robust Match Model

(26) Min Cx
Subject to: c ̃_j x_j+dy≤z ∀c ̃_j∈u_Box^c
a ̃_ij x_j+dy≥b_i ∀a ̃_ij∈u_Box^a
∀f ̃_i∈u_Box^f e_i y≥f ̃_j

Bental et al. [20] proved that in a closed range box, the robust match can change to a solvable model in which u_Box is with a u_ext set that includes final points of u_Box. To present solvable robust mathematical model, the equations should change to solvable mode. For c ̃_j x≤z-dy, we have:
(27) c ̃_j x≤z-dy.∀c ̃_j∈u_Box^c |u_Box^c={c ̃_j∈R^(n_c ):|c ̃_j-c ̅_j |≤ρ_c G_j^c.j=1.2.….n_c }
The equation at the left side includes probabilistic parameters vector while all parameters at the right side are deterministic.
(28) ∑_j▒(c ̅_j x_j+η)
ρ_c G_j^c x_j≤η_j ∀j∈{1.2.….n_c }
ρ_c G_j^c x_j≥-η_j ∀j∈{1.2.….n_c }
For ≥b_i, we need to complete right side parameters, so that they can show the probabilistic nature of parameters. This parameter completion is done as follow:
(29) min┬(a ̃∈u_Box^a )⁡〖{∑_(j=1)^n▒ 〗 a ̃_ij x_j}≥b_i
(30) (min)┬(a ̃_ij:|a ̃_ij-a ̅_ij |≤ρ_a G_ij^a )⁡{ }≥b_i
The presented model in the previous section is based on this method that in Annex (1), its formulation has been proposed.
Epsilon Constraint
This method was proposed by Haimes el al. [21] for the first time and then, its improved versions were developed by Miettnen [22], Ehrgott and Ryan [23]. The basic idea of Epsilon Constraint is that first of all, one of multiple objectives is selected as the main objective of optimization problem and other objectives, in addition to considering upper and lower limit (benefit table constitution), are transferred to problem’s constraints. Therefore, by changing the constraints related to these functions toward the lower limit, all possible Pareto solutions are proposed for multi-objective problem. The general form of Epsilon Constraint is as follow:
Min Z_j (x)
S.t.
Z_k (x)+s_k=ε_k ∀k≠j,x∈X,s_k∈R^+
(31)
Here, sk is an axillary variable that is related to the limitation of the objective function (k). In this study, in order to compare the relative performance of multi-objective imperialist competitive algorithm with the optimized aspect, Epsilon constraint was used.

5-The Proposed Solution

The suggested model is a linear programming mixed by an integer with many limitations and zero and one variables. The total number of variables, integer variables, and limitations of this study increase by the size of the problem. This should be noted that robust matches proposed for the model are nonlinear. According to the above mentioned problem, this is concluded that supply chain design problem is NP-Hard [25]. Also, multi-objective nature of the problem refuses efficient use of the software such as GAMS at large dimensions. Therefore, a meta-heuristic approach is proposed to solve the developed model. The proposed meta-heuristic approach in this study has used non-dominated multi-objective meta-heuristic genetic sorting algorithm and is compared with imperialist competitive algorithm in control issues [26], data clustering issues [27], industrial engineering [28], and fluid engineering [29]. This comparison shows the superiority of non-dominated genetic sorting algorithm.

NSGA-II algorithm
Characteristics of multi-objective genetic algorithm Variables
Software used MATLAB
Population size 1000
Mutation operator Binary race
Intersection operator SCATTER
Generation 100
Period of time 30 minutes
CPU Dual core
Number of replications 50
Cost function Two-objective

This algorithm was proposed by Deb et al. (2002) to solve multi-objective optimization problems. NSGA-II algorithm is an elitist multi-objective evolutionary genetic algorithm. In addition to having a suitable strategy to maintain superior solutions, this algorithm has a defined mechanism to maintain diversity within the population. Studies by Ghavindan et al. (2014) and Hisne and Jone (2014) showed that NSGA-II algorithm that is to solve location problems and designing supply chain, is one of the best meta-heuristic algorithms that presents near-optimal solutions. The genetic algorithm starts working with a number of initial solutions as the initial population. Each population has number of chromosomes that are randomly produced from the solution space. If the quantity is large, a larger part of the solution space is investigated. In this study, the amount of this parameter is considered as 100. In this level, two chromosomes are selected as generator for genocide. Different selection methods lead to different generations. Roulette wheel and Tournament selection method are standard methods for genetic algorithm. Usually for chromosomes of each generation, the cost objectives are estimated and classified in ascending order. Then, the best samples are selected for mating and others are eliminated. is a fraction of . The number of chromosomes that are maintained in each generation equals . = . In fact, each generation includes senior chromosome and – child chromosome that is obtained from parents by mating. In this study, we have . Selection has been done according to the Roulette wheel, so that in this method, better chromosomes have higher chance for being selected that is directly dependent on their fitness. In tournament method, a competition is hold between 2 or 3 chromosomes and the best chromosome is copied to the next population. This is repeated based on the population size. In standard cross operator for the new generation, two generators are selected as father and mother that is the result of mating by two children. Parents must generate – children the complete the new generation. Most of genetic algorithm users, select a random number from normal distribution with the mean of zero. In this method, a value should be chosen for normal distribution. However, those chromosomes are selected for mutation that are not among best chromosomes for next generations. After crossover and mutation operator, the utility function is estimated for each chromosome. Chromosomes are ranked and the best ones are selected again. This cycle is continued to reach optimal solution.
Classification according to crowding distance
Crowding distance approach is proposed to maintain a variety of solutions in Pareto optimum front. To estimate crowding distance attributable to each point on a certain front, points before and after objective functions are selected and according to Fig. 2, a rectangle (two-dimensional) is constituted. It is obvious that if number of functions is more than 2, the points of interest will constitute a cube.

Fig 2. How to estimate crowding distance in a two-objective minimization problem
As shown in Fig 6-2, points on the front of interest are illustrated as bold points and other solutions are shown by hollow circles. To estimate crowding distance of ith solution in this figure, a rectangle has been drawn where solutions are located. The crowding distance for this point equals the average of sides of the rectangle. The shorter crowding distance shows higher density of solutions. It should be noted that in conditions where the problem has more than two objectives, points i-1 and i+1 will not be similar for all objective functions. Estimation steps of crowding distance for solutions on frontier (F) will be as follow:
Calculate number of solutions on frontier (i) and call it (l) (|F|=l). For each i in this set, suppose the initial value of crowding distance (d_i).
Classify the solutions according to each objective function n=1,2,…,M.
For each objective function (m), give long crowding distance to the solutions on frontier border (initial and final points) (d_(I_1^m )=d_(I_l^m )=∞) and to calculate this index for other solutions use relationship (1-2).
(1-2) d_(I_1^m )=d_(I_l^m )+(f_m^((I_(j+1)^m))-f_m^((I_(j-1)^m)))/(f_m^max-f_m^min )
In the above equation, I_j^m shows the ith solution in the classified list according to the objective function (m). The numerator at the right side of the equation (1-2) shows the mth objective function value difference for two neighboring solutions (j). Denominator shows the minimum and maximum differences of mth objective function in the population.
Intersection Operator
In each iteration of algorithm on a part of the existing population (pCrossover), the intersection operator is implemented and new solutions are provided as the number of nPop*pCrossover. Values of each objective function are calculated for all of population members. Here, binary tournament selection was used. In this regard, first, two chromosomes were randomly selected and the chromosome that belonged to lower frontier was used. If both chromosomes belong to a single frontier, the crowding distance is calculated and the chromosome which has larger value is selected.
Mutation Operator
In each iteration of algorithm on a part of the existing population (pMutation), the mutation operator is implemented and new solutions are provided as the number of nPop*pMutation. We calculate the value of each objective function for the provided solutions. To select generator solutions, binary tournament selection was used that was explained in the previous section. Now the current population, population resulted from intersection, and population resulted from mutation are combined together. Non-dominated sorting operation and sorting based on crowding distance were performed on the current population and arrange it. The initial solution enters the iteration after algorithm and the rest of solutions are eliminated. The steps in BSGA-II algorithm are illustrated in Fig 3. The algorithm stops after reaching maximum iteration.

Fig 3. Conducted steps in each iteration of NSGA-II algorithm
5-2. Imperialist Competitive Algorithm
As is obvious, the efficiency of an algorithm and quality of output solutions are dependent on solution representation in solution space. Also, solution representation should allow easy search for solution space. In this study, solution representation is presented as a numerical series. This representation must simultaneously show the locations of blood banks and donation centers, how to allocate customers to blood banks, how to allocate donation centers to the laboratories, and values of transitional products in the network. In this study, solution representation is as Fig 4.

Fig 4. Solution representation
Here, G*I*I shows whole blood transfer values with blood group (G) from the first layer to the second layer; F *G*J*K shows transfer values of product (F) with blood group (G) from the second layer to the third layer, and; F*G*K+L*L shows transfer values of product (F) with the blood group (G) from the third layer to the fourth layer as well as transfers in sub-network of demand zones.
Characteristics of multi-objective imperialist competitive algorithm Variables
Software used MATLAB
Size of initial empire 1000
B 2
A .5
b .2
Duration 30 minutes
CPU Dual core
Number of iterations 50

Generating the Original Empire
Each mode in imperialist competitive algorithm is in a form of an array. Each array includes optimized values. In genetic algorithm technical terms, this array is called chromosome. In this algorithm, the term “country” is used for array. In N-dimensional optimization problem, a country is 1*N array and this array is defined as Country = [p_1,p_3,p_3,…,p_N]. Wherever p_i is an optimized variable, each variable in a country is distinguished as a sociopolitical characteristic such as culture, language, and economic policy. After generating countries, non-domination technique and crowding distance were used to create fronts and sorting the members of frontiers are used. Then, members of frontier (1) are stored in archive. Emperors are selected from the archive of interest according to the predetermined numbers. Then, the rest of answers are allocated to each emperor according to the power that emperor has. To calculate costs or power of each emperor, the objective function was obtained for each emperor. Then, the cost for each objective function was calculated as equation (37).
〖Cost〗_(i,n)=|f_(i,n)^p-f_(i,n)^(p,best) |/(f_(i,total)^(p,max)-f_(i,total)^(p,min) ) (32)
Here, 〖Cost〗_(i,n) is the normalized objective function (i) for emperor (n). Also, f_(i,n)^(p,best), f_(i,total)^(p,max), and f_(i,total)^(p,min) are the best, maximum, and minimum values of objective function in each iteration, respectively. Finally, the level of normalized cost is obtained as equation (38), so that “r” is the objective function.
Toral 〖Cost〗_n=∑_(i=1)^r▒〖Cost〗_(i,n) (33)
The power of each emperor after obtaining normalized cost is calculated as equation (39) and colonies are distributed according to the power of each emperor.
p_n=|(Toral 〖Cost〗_n)/(∑_(i=1)^(N_imp)▒〖Total Cost〗_i )| (34)
Then, the initial number of emperor colonies is determined as equation (40).
〖NC〗_n=round{p_n.N_col} (35)
In this regard, 〖NC〗_n is the initial number of nth emperor colonies, N_col is total number of colonies. 〖NC〗_n is randomly selected and is given to each emperor. The powerful emperor will have more colonies compared to weaker emperors. The condition to stop imperialist is when there is only empire in the world. The general process of imperialist competitive algorithm and its details are presented in Fig. 5, respectively.

Fig 5. General schematic of imperialist competitive algorithm

6-Computational Results

Three proposed problems at small sizes are provided in Table (1). Selection of these sizes has been based on the available problems in literature, because the model has two objectives [30]. The conflicts are clearly illustrated in the table that show the validity of two-objective model. To solve these problems, GAMS has been used and in large size, non-dominated genetic sorting algorithm was used.
Table 1. Comparing the first and second objective functions at absolute mode
Problem size
I/J/K/L First objective function value Second objective function value
3/2/2/3 6.5862E+1 18(h)
5/3/2/4 1.5858E+8 115 (h)
10/4/2/6 3.3025E+8 265 (h)
To investigate the performance of the first probabilistic model, three problems were selected with different numbers of demand zones, blood banks, laboratories, and donation centers and for each size, the model is solved with 5 different levels (ρ=0.2,0.4,0.6,0.60,0.8,1) of probability in the robust model.
Table 2. Results of robust match
problem size
I/J/K/L Objective function Probability level
ρ Objective function values Calculation time (second)
Definite Robust Definite Robust
3/2/2/3 First 0.2 6.5862E+1 3.3100E+8 15 18
0.4 6.0211E+8 19
0.6 8.7323E+8 21
0.8 1.1444E+9 14
1 1.4155E+9 14
Second 0.2 18(h) 93.6 (h) 0.000 0.000
0.4 109.2 (h) 0.000
0.6 124.8 (h) 0.000
0.8 140.4 (h) 0.000
1 156 (h) 0.000
5/3/2/4 First 0.2 1.5858E+8 5.1696E+8 165 4620
0.4 9.3654E+8 4753
0.6 1.6096E+9 4698
0.8 2.1257E+9 4785
1 2.9625E+9 4658
Second 0.2 115 (h) 138 (h) 65 136
0.4 149.1 (h) 139
0.6 162.6 (h) 175
0.8 184 (h) 149
1 207.5 (h) 156
10/4/2/6 First 0.2 3.3025E+8 – 4202 –
0.4 – –
0.6 – –
0.8 – –
1 – –
Second 0.2 265 (h) – 521 –
0.4 – –
0.6 – –
0.8 – –
1 – –
To compare the efficiency of proposed algorithms, the model will be implemented for average-and-large-sized problems using both algorithms and are compared by four indicators mentioned in the previous section. The results of comparing imperialist competitive algorithm and non-dominated genetic sorting algorithm are presented in Table (3). All of these results are solved by probabilistic degree of p=0.8 that is explained in robust model. In Annex (2), criteria calculation is explained.
Table 3. Comparing MOICA and NSGA-II algorithms for the model
MID DM QM SM Comparison index
MOICA NSGA-II MOICA NSGA-II MOICA NSGA-II MOICA NSGA-II Problem size
0.581 0.403 1.209 1.322 0.705 0.875 0.682 0.346 3/2/2/3
0.714 0.415 1.145 1.457 0.554 0.756 0.652 0.441 5/3/2/4
0.638 0.531 1.231 1.568 0.405 0.501 0.644 0.488 10/4/2/6
0.787 0.720 0.732 0.878 0.433 0.647 0.708 0.495 25/7/4/11
0.746 0.633 0.690 0.952 0.318 0.417 0.771 0.505 45/12/7/21
For statistical comparison of algorithms, pairwise t-test was used and a significant difference was observed between algorithms. This test showed confidence level of 90% that can be observed in Table (4).
Table 4. Pairwise t-test to compare two algorithms
Criteria Mean Standard deviation Standard error mean Confidence level of 95% t Degree of freedom Sig
Lower limit Upper limit
SM 0-.23640 0.06792 0.03038 0-.32074 0-.15206 -7.782 4 0.001
QM 0.15620 0.05596 0.02502 0.08672 0.22568 6.241 4 0.003
DM 0.23400 0.09983 0.04464 0.11005 0.35795 5.241 4 0.006
MID 0-.15280 0.09090 0.04065 0-.26567 0-.03993 -3.758 4 0.02

According to the significance of pairwise t-tests and this fact that their levels are smaller than (0.05), it can be concluded that the proposed algorithm at all four levels has a significant difference with NSGA-II algorithm. Therefore, it is possible to use NSGA-II algorithm with more confidence in other problems. In figures (6) to (10), the Pareto points related to both algorithms for problems with different sizes are proposed simultaneously to observe the superiority of NSGA-II algorithm compared to the proposed algorithm. Therefore, this algorithm is more suitable for solution. To observe the quality of the proposed algorithm in providing acceptable solutions, Pareto solution resulted from Epsilon constraint for problems with small sizes is compared with the proposed algorithms.

Fig 6. Comparing the Pareto solutions of proposed algorithm with Epsilon constraint algorithm for size problem 3/2/3/2

Fig 7. Comparing Pareto solutions of the proposed algorithm with Epsilon constraint algorithm for size problem 5/3/2/4

Fig 8. Comparing Pareto solutions of the proposed algorithms with Epsilon constraint algorithm for size problem 10/4/2/6

Fig 9. Comparing Pareto solutions of the proposed algorithm with Epsilon constraint algorithm for size problem 25/7/4/12

Fig 10. Comparing Pareto solutions of the proposed algorithm with Epsilon constraint algorithm for size problem 45/12/7/21

7-Sensitivity Analysis

The mathematical model solution gives a schematic view of problem solving in simulated conditions. However, the effects of basic parameters in model are not considered. In this section, the effects of parameters on objective function are investigated. Those parameters that are considered to analyze sensitivity include the loss rate, demand, supply or blood donation, and probability level. Sensitivity analysis in this section has been done considering probabilistic level (0.8) for the problem with the size of 3/2/2/3.
Pareto solutions
Fig. 11 shows Pareto optimized solutions obtained from Epsilon constraint. This figure has been obtained for the first problem at the uncertainty level of 0.8. As expected, by increased first objective function, the second objective function decreased that show the contradiction between objective functions.

Fig 11. Pareto solution of Epsilon algorithm for the model
Loss rate
The first parameter that is considered to analyze sensitivity is blood loss rate at the laboratory that is shown by α. Blood loss is caused due to various reasons; reasons such as product transfer space, blood diseases, hemolysis, and other unknown causes. According to Fer et al. [24], almost 779% of bloods losses are due to expiration. Blood loss in hospitals ranges between 193% and 307%. In general, it is not possible to measure exact amount of blood loss. In this study, the nominal amount of blood loss has been considered as 10%. Fig. 12 shows changes of cost objective function with changes in loss rate from 10% to 30% at 5% intervals. Time objective function is not dependent on loss rate. As can be seen, increased loss rate in laboratories, increases cost objective function.

Fig 12. Effect of blood loss rate changes on cost objective function
Uncertainty level
Uncertainty level in parameters can have significant effect on objective functions. These changes are shown in Figs (13) and (14) for five uncertainty levels (ρ=0.2,0.4,0.6,0.60,0.8,1) on two objective and cost functions. Obviously, with increased uncertainty, both cost and time functions significantly increase.

Fig 13. The effect of uncertainty level on cost objective function value

Fig 14. The effect of uncertainty level on cost objective function value

Blood donation amount
The third parameter that will be analyzed regarding sensitivity is blood supply or blood donation. This issue that donation rate in collection centers is highly random is proved and as we know, only a small part of potential donors, donate blood. This can unpleasantly change the supply chain. Therefore, blood supply is a very important issue. Fig. 15 shows how to change optimal front with decreased blood provision. In this figure, donation amount in each solution has decreased by 10%. This should be noted that decreased blood donation causes loss in chain. Therefore, cost function value is directly related to penalty rate. This is true for demands, as well.

Fig 15. The effect of blood donation rate on cost objective function
Demand
The last sensitivity analysis is done on demand parameter. Blood demand is dependent on various factors such as population, age, gender, accidents, and unexpected events. Therefore, blood demand prediction requires attempts and it is likely that unreliable results are provided. Problem sensitivity analysis regarding demands such as blood supply has been done, so that in each solution set, demand level increase by 10%. By increased demand, it means increased demand for all model products. With increased demand, the cost objective function increases that can be observed in Fig. 16.

Fig 16. The effect of blood demands changes on cost objective function

8-Conclusion

Blood management is an issue that has been taken into consideration. Although there have been numerous advances in technology, no alternative has been obtained for blood products. Blood is not a normal product. Blood supply by donors is not normal and demands for blood products are probabilistic at the best mode. Indeed, demand and supply match in blood supply chain is not an easy task. Also, blood products are perishable and this has complicated the situation. Also, it should be noted that since blood shortage can increase mortality, it causes additional costs and damages for the society. The reasons mentioned above only constitute a part of reasons that encourage researchers to conduct studies on blood supply context and more specifically, blood supply chain. In this study, mixed integer non-linear programming is proposed for designing blood supply chain. This supply chain is a four-level supply chain consisting of blood donation centers, blood laboratories, blood bank center, and demand centers. The goals of model include: 1. Locating donation centers and blood bank centers 2. Determining the amount of transmission products between both points of supply chain. In this study, most of problem complexities are considered. Complexities such as blood loss in laboratories, blood products decomposition in the laboratory, and transferring blood products among blood products. Since model parameters have uncertain nature in real world, their amounts are not clear. In this study, the model parameters are uncertain. To solve uncertainty problem of parameters, robust optimization method was used in this study. The research model was analyzed in GAMS using certain method. However, because in larger sizes, the problem cannot be solved by exact methods, imperialist competitive meta-heuristic algorithm and multi-objective genetic algorithm were used. After solving the model by these methods, it was identified that non-dominated genetic sorting method has better performance. In other words, the model of this study is two-objective that the first objective minimizes blood supply chain designing costs and the second objective function minimizes the time that a blood product remains in the network. In the proposed model, the considered supply chain has four products including whole blood, plasma, platelet, and red blood cells. In this model, complex problems of blood group match are considered. Therefore, in situations where a blood group is not available sufficiently, those blood groups that can be replaced are used. This model was formulized non-linearly, but became linear by proposing a method. Since this model is two-objective, Epsilon constraint was used to obtain Pareto solution; then, to solve the model, two meta-heuristic algorithms were used.

References

• J. Beliën and H. Forcé, “Supply chain management of blood products: A literature review,” European Journal of Operational Research, vol. 211, pp. 16-1, 2012.
• W. P. Pierskalla, “Supply chain management of blood banks,” in Operations research and health care, ed: Springer, pp. 145-103, 2005.
• K. Katsaliaki, “Cost-effective practices in the blood service sector,” Health policy, vol.86, pp. 281-216, 2008.
• G. B. Schreiber, K. S. Schlumpf, S. A. Glynn, D. J. Wright, Y. Tu, M. R. King, et al., “Convenience, the bane of our existence, and other barriers to donating,” Transfusion, vol. 46, pp. 553-545, 2006.
• J. S. Rytilä and K. M. Spens, “Using simulation to increase efficiency in blood supply chains,” Management Research News, vol. 21, pp. 811-801, 2006.
• G. P. Prastacos, “Blood inventory management: an overview of theory and practice,” Management Science, vol. 30, pp. 800-111, 1984.
• J. Rautonen, “Redesigning supply chain management together with the hospitals,” Transfusion, vol. 41, pp. 111S-200S, 2007.
• J. R. Hess, “Red cell freezing and its impact on the supply chain,” Transfusion Medicine, vol. 14, pp. 8-1, 2004.
• N. Van Dijk, R. Haijema, J. Van Der Wal, and C. S. Sibinga, “Blood platelet production: a novel approach for practical optimization,” Transfusion, vol. 41, pp. 420-411, 2009.
• J. T. Blake, “On the use of Operational Research for managing platelet inventory and ordering,” Transfusion, vol. 41, pp. 401-316, 2009.
• M. S. Daskin, C. R. Coullard, and Z.-J. M. Shen, “An inventory-location model: Formulation, solution algorithm and computational results,” Annals of Operations Research, vol. 110, pp. 106-83, 2002.
• Z.-J. M. Shen, C. Coullard, and M .S. Daskin, “A joint location-inventory model,” Transportation Science, vol. 31, pp. 55-40, 2003.
• E. Cetin and L. S. Sarul, “A blood bank location model: A multiobjective approach,” European Journal of Pure and Applied Mathematics, vol. 2, pp. 124-112, 2009.
• G. Şahin, H. Süral, and S. Meral, “Locational analysis for regionalization of Turkish Red Crescent blood services,” Computers & operations research, vol. 34, pp. 104-612, 2007.
• Y. Sha and J. Huang, “The Multi-period Location-allocation Problem of Engineering Emergency Blood Supply Systems,” Systems Engineering Procedia, vol. 5, pp. 28-21, 2012.
• M. J. Fontaine, Y. T. Chung, W. M. Rogers, H. D. Sussmann, P. Quach, S. A. Galel, et al., “Improving platelet supply chains through collaborations between blood centers and transfusion services,” Transfusion, vol. 41, pp. 2041-2040, 2009.
• A.Jabbarzadeh, B. Fahimnia, and S. Seuring, “Dynamic supply chain network design for the supply of blood in disasters: A robust model with real world application,” Transportation Research Part E: Logistics and Transportation Review, Vol. 70, pp. 225-244, 2014.
• Ben-Tal and A. Nemirovski, “Robust convex optimization,” Mathematics of Operations Research, vol. 23, pp. 805-161, 1998.
• Vira and Y. Y. Haimes, Multiobjective decision making: theory and methodology: North-Holland, 1983.
• K. Miettinen, Nonlinear multiobjective optimization vol. 12: Springer, 1999.
• M. Ehrgott and D. M .Ryan, “Constructing robust crew schedules with bicriteria optimization,” Journal of Multi‐Criteria Decision Analysis, vol. 11, pp. 150-131, 2002.
• R. M. Far, F. S. Rad, Z. Abdolazimi, and M. M. D. Kohan, “Determination of Rate and Causes of Wastage of Blood and Blood Products in Iranian Hospitals,” Turkish Journal of Hematology, vol. 31, p. 161, 2014.
• Zhou, G., Min, H., & Gen, M. (2002). The balanced allocation of customers to multiple distribution centers in the supply chain network: a genetic algorithm approach. Computers & Industrial Engineering, 43(1), 251-261.
• Lucas, C., Nasiri-Gheidari, Z., & Tootoonchian, F. (2010). Application of an imperialist competitive algorithm to the design of a linear induction motor. Energy conversion and management, 51(7), 1407-1411.
• Niknam, T., Fard, E. T., Pourjafarian, N., & Rousta, A. (2011). An efficient hybrid algorithm based on modified imperialist competitive algorithm and K-means for data clustering. Engineering Applications of Artificial Intelligence, 24(2), 306-317.
• Nazari-Shirkouhi, S., Eivazy, H., Ghodsi, R., Rezaie, K., & Atashpaz-Gargari, E. (2010). Solving the integrated product mix-outsourcing problem using the imperialist competitive algorithm. Expert Systems with Applications, 37(12), 7615-7626.
• Yousefi, M., Darus, A. N., & Mohammadi, H. (2012). An imperialist competitive algorithm for optimal design of plate-fin heat exchangers. International Journal of Heat and Mass Transfer, 55(11), 3178-3185.
• Fahimnia, B., Jabbarzadeh, A., Ghavamifar, A., & Bell, M. (2015). Supply chain design for efficient and effective blood supply in disasters. International Journal of Production Economics.
• Nozari, H.; Najafi, E.; Fallah, M.; Hosseinzadeh Lotfi, F. Quantitative Analysis of Key Performance Indicators of Green Supply Chain in FMCG Industries Using Non-Linear Fuzzy Method. Mathematics 2019, 7, 1020. https://doi.org/10.3390/math7111020

Annex 1

Blood supply chain network designing model under uncertainty
According to robust model description, robust match model, and blood supply chain designing are formulized as follow considering blood groups match.
1 Min Z_1
2
Min Z_2
∑_(i∈I)▒g_i y_i+∑_(i∈I)▒〖(f ̅_i/λ〗 y_i+η_i^f)+∑_(k∈K)▒g_k y_k+∑_(k∈K)▒〖(f ̅_k/λ〗 y_k+η_k^f)+∑_(i∈I)▒∑_(j∈J)▒ρ_ij z_ij+∑_(i∈I)▒〖∑_(j∈J)▒〖(h ̅_ij/λ〗 z_ij+η_ij^h)+∑_(j∈J)▒∑_(k∈K)▒ρ_jk z_jk+∑_(j∈J)▒〖∑_(k∈k)▒〖(h ̅_jk/λ〗 z_jk+η_jk^h)〗+∑_(k∈K)▒∑_(l∈L)▒ρ_kl z_kl+∑_(k∈K)▒〖∑_(l∈L)▒〖(h ̅_kl/λ〗 z_kl+η_kl^h)+〗 ∑_(l∈L)▒∑_(m∈L)▒ρ_lm z_lm+∑_(l∈L)▒〖∑_(m∈L)▒〖(h ̅_lm/λ〗 z_lm+η_lm^h)〗+∑_(g∈G)▒∑_(i∈I)▒∑_(j∈J)▒(C ̅_ij X_ij^g η_ij^c ) +∑_(g∈G)▒∑_(f∈F)▒∑_(j∈J)▒∑_(k∈K)▒〖(C ̅_jk X_jk^fg η_jk^c )+∑_(g∈G)▒∑_(f∈F)▒∑_(k∈K)▒∑_(l∈L)▒(C ̅_kl X_kl^fg η_kl^c ) +∑_(g∈G)▒∑_(f∈F)▒∑_(l∈L)▒∑_(m∈L)▒(C ̅_lm X_lm^fg η_lm^c ) +M×∑_f▒∑_g▒∑_l▒〖slack〗_1^fg ≤Z_1 〗 〗 3
∑_f▒∑_i▒∑_j▒∑_k▒∑_l▒∑_m▒〖(t ̅_ij z_ij+η_ij^t+t ̅_jk z_jk+η_jk^t+t ̅_kl z_kl+η_kl^t+t ̅_i y_i+η_i^t+t ̅_j^f+t ̅_k^f y_k+η_k^h 〗+(t ̅_lm+t ̅_l^f )×z_lm+η_lm^h+η_l^h)≤Z_2 4
ρ_f ξ_i^f y_i≤η_i^f λ ∀i∈I 5
ρ_f ξ_i^f y_i≥-η_i^f λ ∀i∈I 6
ρ_f ξ_k^f y_k≤η_k^f λ ∀k∈K 7
ρ_f ξ_k^f y_k≥-η_k^f λ ∀k∈K 8
ρ_h ξ_ij^h z_ij≤η_ij^h λ ∀i∈I . ∀j∈J 9
ρ_h ξ_ij^h z_ij≥-η_ij^h λ ∀i∈I . ∀j∈J 10
ρ_h ξ_jk^h z_jk≤η_jk^h λ ∀j∈J . ∀k∈K 12
ρ_h ξ_jk^h z_jk≥-η_jk^h λ ∀j∈J . ∀k∈K 13
ρ_h ξ_kl^h z_kl≤η_kl^h λ ∀k∈K . ∀l∈L 14
ρ_h ξ_kl^h z_kl≥-η_kl^h λ ∀k∈K . ∀l∈L 15
ρ_h ξ_lm^h z_lm≤η_lm^h λ ∀l∈L . ∀m∈L 16
ρ_h ξ_lm^h z_lm≥-η_lm^h λ ∀l∈L . ∀m∈L 17
ρ_c ξ_ij^c X_ij^g≤η_ij^c ∀i∈I . ∀j∈J . ∀g∈G 18
ρ_c ξ_ij^c X_ij^g≥-η_ij^c ∀i∈I . ∀j∈J . ∀g∈G 19
ρ_c ξ_jk^c X_jk^fg≤η_jk^c ∀j∈J . ∀k∈K .∀f∈F .∀g∈G
20
ρ_c ξ_jk^c X_jk^fg≥-η_jk^c ∀j∈J . ∀k∈K .∀f∈F .∀g∈G
21
ρ_c ξ_kl^c X_kl^fg≤η_kl^c ∀j∈J . ∀k∈K .∀f∈F .∀g∈G
22
ρ_c ξ_kl^c X_kl^fg≥-η_kl^c ∀k∈K . ∀l∈L .∀f∈F .∀g∈G
23
ρ_c ξ_lm^c X_lm^fg≤η_lm^c ∀l∈L . ∀m∈L .∀f∈F .∀g∈G
24
ρ_c ξ_lm^c X_lm^fg≥-η_lm^c ∀l∈L . ∀m∈L .∀f∈F .∀g∈G
25
ρ_t ξ_ij^t z_ij≤η_ij^t ∀i∈I . ∀j∈J 26
ρ_t ξ_ij^t z_ij≥-η_ij^t ∀i∈I . ∀j∈J 27
ρ_t ξ_jk^t z_jk≤η_jk^t ∀j∈J . ∀k∈K 28
ρ_t ξ_jk^t z_jk≥-η_jk^t ∀j∈J . ∀k∈K 29
ρ_t ξ_kl^t z_kl≤η_kl^t ∀k∈K . ∀l∈L 30
ρ_t ξ_kl^t z_kl≥-η_kl^t ∀k∈K . ∀l∈L 31
ρ_t ξ_lm^t z_lm≤η_lm^t ∀l∈L . ∀m∈L 32
ρ_t ξ_lm^t z_lm≥-η_lm^t ∀l∈L . ∀m∈L 33
ρ_t ξ_lm^t z_lm≤η_lm^t ∀l∈L . ∀m∈L 34
ρ_t ξ_lm^t z_lm≥-η_lm^t ∀l∈L . ∀m∈L 35
ρ_t ξ_i^t y_i≤η_i^t ∀i∈I 36
ρ_t ξ_i^t y_i≥-η_i^t ∀i∈I 37
ρ_t ξ_k^t y_k≤η_k^t ∀k∈K 38
ρ_t ξ_k^t y_k≥-η_k^t ∀k∈K 39
(1-a ̅-ρ_a ξ^a ) ∑_(i∈I)▒X_kl^fg ≥∑_(k∈K)▒X_jk^4g +〖X^”〗_ij^fg ∀f∈F-{4}.∀g∈G.∀j∈J
40
〖X^”〗_ij^fg≥X_jk^fg ∀f∈F-{4}.∀g∈G.∀j∈J
41
∑_(j∈J)▒X_jk^fg ≥∑_(l∈L)▒〖X_kl^fg 〗 ∀f∈F.∀g∈G.∀k∈K
42
∑_(k∈K)▒X_kl^fg ≥∑_(m∈L)▒X_lm^fg ∀f∈F.∀g∈G.∀l∈L
43
∑_(j∈J)▒X_ij^g ≤〖dc〗_i^g×y_i ∀i∈I.∀g∈G
44
∑_(g∈G)▒∑_(f∈F)▒∑_(k∈K)▒X_jk^fg ≤〖lc〗_j ∀j∈J 45
∑_(g∈G)▒∑_(f∈F)▒∑_(l∈L)▒X_kl^fg ≤〖cbb〗_k×〖y^’〗_k ∀k∈K 46
∑_(g∈G)▒∑_(f∈F)▒∑_(k∈K)▒X_kl^fg +∑_(g∈G)▒∑_(f∈F)▒∑_(m∈L)▒X_ml^fg ≤〖hc〗_l ∀l∈L 47
∑_(k∈K)▒X_kl^fg +∑_(m∈L)▒X_ml^fg -∑_(m∈L)▒X_lm^fg ≥0 ∀f∈F .∀g∈G.∀l∈L
48
∑_g▒〖{〖BC〗^fgp 〗[∑_(k∈K)▒〖X^’〗_kl^fgp +∑_(m∈L)▒〖〖(X〗^’〗_ml^fgq -〖X^’〗_lm^fgq)]}+〖slack〗_l^fg=d ̅_I^fg+ρ_d ξ_(fgl )^d ∀f∈F .∀p∈G.∀l∈L
49
∑_g▒〖{〖BC〗^fgp 〗[∑_(k∈K)▒〖X^’〗_kl^fgp +∑_(m∈L)▒〖〖(X〗^’〗_ml^fgq -〖X^’〗_lm^fgq)]}+〖slack〗_l^fg=d ̅_I^fg+ρ_d ξ_(fgl )^d ∀f∈F .∀p∈G.∀l∈L
50
X_kl^fg≥∑_(p∈G)▒〖〖X^’〗_kl^fgp 〗 ∀f∈F . ∀p∈G .∀k∈K .∀l∈L
51
X_lm^fg≥∑_(p∈G)▒〖〖X^’〗_lm^fgp 〗 ∀f∈F . ∀p∈G .∀l∈L .∀m∈L
52
(t ̅_ij+ρ_t ξ_ij^t ) z_ij+(t ̅_jk+ρ_t ξ_jk^t ) z_jk+(t ̅_kl+ρ_t ξ_kl^t ) z_kl+(t ̅_i+ρ_t ξ_i^t ) y_i+t ̅_j^f+ρ_t ξ_j^t+(t ̅_k^f+ρ_t ξ_k^t ) y_k+(t ̅_lm+ρ_t ξ_lm^t+t ̅_l^f+ρ_t ξ_l^t )×z_lm≤T^f ∀f∈F . ∀i∈I .∀j∈J .∀k∈K.∀l∈L .∀m∈L
53
T^f≤T_max^f ∀f∈F 54
X_ij^g≤〖MY〗_i ∀g∈G.∀i∈I.∀j∈J
55
X_jk^fg≤〖MY^’〗_k ∀f∈F.∀g∈G.∀j∈J.∀k∈K 56
X_kl^fg≤〖MY^’〗_k ∀f∈F.∀g∈G.∀k∈K.∀l∈L 57
X_ij^g≤〖Mz〗_(ij ) ∀g∈G.∀i∈I.∀j∈J
58
X_jk^fg≤〖Mz〗_jk ∀f∈F.∀g∈G.∀j∈j.∀k∈K 59
X_kl^fg≤〖Mz〗_kl ∀f∈F.∀g∈G.∀k∈K.∀l∈L
60
X_lm^fg≤〖Mz〗_lm ∀f∈F.∀g∈G.∀l∈L.∀m∈L 61
X_ij^g .X_jk^fg .X_kl^fg . X_lm^fg . η_i^f . η_k^f . η_ij^h . η_jk^h . η_kl^h . η_lm^h . η_ij^c . η_jk^c . η_kl^c . η_lm^c . η_ij^t . η_jk^t . η_kl^t . η_lm^t . η_i^t . η_k^t≥0 ∀f∈F.∀g∈G.∀i∈I.∀j∈J.∀k∈K.∀l∈L.∀m∈L 62
Z_ij .Z_jk .Z_kl .Z_lm .y_i .y_k∈{0.1}
∀i∈I.∀j∈J.∀k∈K.∀l∈L.∀m∈L 63

Annex 2

To compare proposed algorithms, various mtrices will be used that in this section, a brief description is provided.
Spacing Metric (SM)
This metric shows uniform distribution of Pareto solutions. The calculation of this metric is according to Equation (1).
SM=(∑_(i=1)^(n-1)▒|d ̅-d_i | )/((n-1)d ̅ ) (1)
Here, d is the Euclidean distance between two Pareto solutions in the present space. Also, d ̅ is the mean of d_i distances. The smaller SM, the better algorithm.
Quality Metric (QM)
The QM considers all obtained Pareto solutions by each of the algorithms. Then, non-dominated operation is done for all solutions. Finally, the quality of each algorithm equals the share of Pareto solutions. High quality shows fitness of algorithm.
Diversification Metric (DM)
This metric shows the extent of Pareto solutions of an algorithm and can be calculated by Equation (2). The larger DM, the better algorithm.
DM=√(((maxf_1i-minf_1i)/(f_(1,total)^max-f_(1,total)^min ))^2+((maxf_2i-minf_2i)/(f_(2,total)^max-f_(2,total)^min ))^2 ) (2)

Mean Ideal Distance (MID)
The amount of this metric equals Pareto points distance with algorithm of interest from ideal point. In this study, the ideal point is considered as minimum level of objective functions in all algorithms, since objective functions consider minimization. MID metric can be calculated by Equation (3).
MID=(∑_(i=1)^n▒√(((f_1i-f_1^best)/(f_(1,total)^max-f_(1,total)^min ))^2+((f_2i-f_2^best)/(f_(2,total)^max-f_(2,total)^min ))^2 ))/n (3)