A robust possibilistic programming approach to ... - Wiley Online Library

0 downloads 0 Views 1MB Size Report
hospital evacuation planning problem under uncertainty ... problem, two well-known metaheuristic algorithms are developed to solve the large-sized problems.
Intl. Trans. in Op. Res. 25 (2018) 157–189 DOI: 10.1111/itor.12331

INTERNATIONAL TRANSACTIONS IN OPERATIONAL RESEARCH

A robust possibilistic programming approach to multiperiod hospital evacuation planning problem under uncertainty Masoud Rabbania , Mohammad Zhalechiana and Amir Farshbaf-Geranmayeha,b a

School of Industrial Engineering, College of Engineering, University of Tehran, Tehran, Iran b Group for Research in Decision Analysis (GERAD), HEC Montr´eal, Montr´eal, Canada E-mail: [email protected] [Rabbani]; [email protected] [Zhalechian]; [email protected] [Farshbaf-Geranmayeh] Received 3 July 2015; received in revised form 22 April 2016; accepted 21 June 2016

Abstract In this paper, a biobjective programming model is developed to address the hospital evacuation problem under uncertainty. It aims to concurrently minimize the total evacuation time and the total weighted number of unevacuated patients in each period. The presented model considers two types of patients and three transportation modes. Moreover, the evacuating hospitals are divided into two groups. In the first group, it is not possible to send vehicles to the evacuating hospitals due to the poor road condition or congestion, whereas there is no such limitation in the second group. A robust possibilistic programming approach is adopted to deal with the inherent uncertainty in the input data. To cope with the computational complexity of the problem, two well-known metaheuristic algorithms are developed to solve the large-sized problems. Finally, several computational experiments and sensitivity analyses are conducted and the results are analyzed. Keywords: hospital evacuation planning; disaster; robust possibilistic programming; metaheuristics

1. Introduction Hospital evacuation is a major issue that hospital’s emergency management is confronting. In the past, the most common reasons for evacuating the hospitals were related to natural disasters such as floods, hurricanes, tornados, and earthquakes. Also, more recently, a variety of potential threats such as loss of power and exposure to hazardous material create the need to evacuate the endangered hospitals (Childers et al., 2014). Several important reasons highlighted the need for designing emergency response procedure, including the large number of health care facilities (i.e., hospitals), their placement and geographical presence in the population, and the huge level of risk to patients associated with evacuation (McGlown, 2001). Planning of evacuations requires a careful strategy and has to be available or generated quickly in the response phase. In this phase, a set of actions are taken to handle the difficult process of hospital evacuation.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research  Published by John Wiley & Sons Ltd, 9600 Garsington Road, Oxford OX4 2DQ, UK and 350 Main St, Malden, MA02148, USA.

158

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Hospital evacuation planning is complex because it involves many factors such as limited resources and specific needs of patients. In reality, patients are categorized according to their treatment level (Najafi et al., 2013). Usually, patients are categorized into two major types, namely, critical and noncritical patients. Critical patients need specific types of transportation vehicles (e.g., ambulances) and they have a higher priority than the other type. Resource allocation is another important issue in the hospital evacuation problem. The risk manager, who is responsible for the evacuation of patients of a city or a country, usually faces with some limitations such as the number of available vehicles that can be used to transfer patients (Tayfur and Taaffe, 2009a). Due to the complexity of evacuation planning, using a mathematical model that considers aforementioned factors can play an essential role in planning and executing a successful evacuation. Research in health care facility evacuation can be divided into two main categories: (1) behavioral and social science, (2) modeling and operations. In the behavioral and social science category, Vogt (1991) studied the organizational characteristics of nursing homes and health care facilities in which a partial or complete evacuation has occurred. There are other researches in this category that are focused on the decision-making process for hospitals and special needs of evacuees (McGlown, 1999, 2001). In the literature, there are a few studies in the second category (i.e., modeling and operations). Tayfur and Taaffe (2009a) presented a deterministic optimization approach to find the scheduling and allocation of resources within a prespecified evacuation time period while minimizing the total cost. In another study, they presented a stochastic model via simulation to address hospital evacuation. They used a simulation-optimization framework to examine nurse and vehicle transport requirements while minimizing the total cost (Tayfur and Taaffe, 2009b). There are other studies focused on patient prioritization. Childers et al. (2009) utilized simulation to investigate different scenarios and obtain optimal evacuation policies while considering the rate of evacuation and survivability. Based on their previous study, Childers et al. (2014) presented a dynamic programming model to find the optimal prioritization strategies while considering different factors including the rate of the evacuating, probability of surviving the transport, and rate of dying while waiting for the evacuation. Notably, resource restriction such as vehicles and available beds in the receiving hospital were not considered in this research. Bish et al. (2014) presented a model to determine patient allocation strategies for transporting patient from evacuating hospitals to a set of potential receiving hospitals. They considered two objective functions, one of them aims to minimize the evacuation risk and the other one aims to minimize the sum of the cumulative threat risk for the patients who are not evacuated during the time period. In their model, they considered available beds in receiving hospitals, vehicles limitation for transporting patients, and different types of patients. There are several studies in the disaster management area that emphasized the importance of air transportation in disaster response phase. Sometimes, emergency affects the availability of transportation resources/routes or there is a need to transport patients to receiving hospitals in a faster way (Childers and Taaffe, 2010). In this respect, other types of transportation such as air transportation can be helpful to evacuate patients. Haghani and Oh (1996) and Oh and Haghani (1997) analyzed the transportation of personnel and different disaster relief commodities, such as food, medicine, and machinery, to maximize the efficiency of the rescue operations. Barbaroso˘glu et al. (2002) focused on the use of helicopters for rescue missions and aid delivery during the initial response phase of disaster management. Ozdamar (2011) coordinated helicopter operations in disaster relief to evacuate injured persons and deliver items such as medicine, vaccines, etc.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

159

Most of the previous studies assumed that vehicles can easily access to the hospitals to pick up patients. However, in many real-world cases, some barriers will challenge this assumption (White et al., 2008). Indeed, some roads that ran through evacuating hospitals may suffer great damage in the aftermath of some disasters as they become blocked by toppled houses, fires, etc. (e.g., Japan’s Great Hanshin-Awaji Earthquake in January 1995, Australia’s Grafton Floods in March 2001). As a result, a great deal of car traffic converged in the limited roads may cause a huge congestion. Such congestion makes it difficult for rescue vehicles to access to the evacuating hospitals (Chang and Winton, 2001). In this situation, pick-up nodes that exist in the vicinity of hospitals can be helpful to transfer patients. In this way, rescue groups can transfer patient from evacuating hospital to pick-up nodes and then patients can be transferred to receiving hospitals via available vehicles at pick-up nodes. In the real world, many parameters such as the demands, respective costs, travel times, and other relevant parameters may change due to the dynamic nature of the environment. Therefore, it is important to make the obtained solution more realistic by addressing the uncertainties. To this end, three types of modeling techniques, namely stochastic programming, fuzzy programming, and robust optimization are introduced in the literature. In this paper, a robust possibilistic programming (RPP) approach that benefits from both advantages of possibilistic and robust programming approaches is used. The literature review demonstrates that there is a gap in the consideration of alternative transportation modes with different transportation time in hospital evacuation planning. Furthermore, most of the existing researches assumed that it is always possible to send vehicles to the evacuating hospitals to pick up patients, while in the real life it is not possible to send vehicles to some of the evacuating hospitals due to the poor road condition or congestion. Additionally, existing literature does not sufficiently address the uncertainties in the hospital evacuation problem. Moreover, there is not any study in the literature that simultaneously minimizes the total evacuation time and the total weighted number of unevacuated patients in each period. In order to fill these gaps, this paper contributes to the literature in the following ways:

r Developing a biobjective mixed integer programming model to minimize the total evacuation time and the total weighted number of unevacuated patients in each period.

r Considering two types of evacuating hospitals. In the first group, it is not possible to send vehicles r r r r

to evacuating hospitals due to the poor road condition or congestion. In the second group of evacuating hospitals vehicles can easily access to the hospitals to pick up patients. Considering two types of transportation, which are land path and aerial path to transport different types of patients. Allocating available vehicles (i.e., ambulances, vans, and helicopters) to pick-up nodes, helicopter stations, and second group hospitals in an efficient way. Applying an efficient RPP approach to cope with the uncertainties. Developing two metaheuristic algorithms to solve the large-sized instances.

The rest of the paper is organized as follows: in Section 2, problem definition and mathematical model are presented. The proposed RPP is explained in Section 3. In Section 4, solution methods are explained. Section 5 handles the computational experiments to show the efficiency of model. The conclusions are presented in Section 6.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

160

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

2. Problem description and mathematical programming 2.1. Problem description In general, the main objective of the hospital evacuation plan in the response phase is transferring patients from evacuating hospitals to receiving hospitals. In this situation, the available vehicles to transfer patients are usually inadequate. Furthermore, there are many factors that should be considered in the evacuation plan including, capacity of receiving hospitals, conditions and care requirements of patients, etc. In addition to the complexity of an evacuation plan, the available network data on demands (i.e., the number of evacuees in the evacuating hospitals), the number of available vehicles, and the transportation times between each node in the response phase are usually uncertain during the planning horizon. These uncertainties further complicate the evacuation planning. Hence, this paper aims at developing a biobjective, multimode and multiperiod mathematical model for hospital evacuation problem in the disaster response phase while considering uncertainty. Figure 1 shows the schematic view of a simple example of the considered network where there are five evacuation hospitals, five pick-up nodes, two helicopter stations, and four receiving hospitals. The main assumptions of the proposed model are as follows:

r Evacuating/receiving hospital nodes, pick-up nodes, and helicopter stations are known. r The evacuating hospitals are divided into two groups. In the first group, it is not possible to send vehicles to the evacuating hospitals due to the poor road condition or congestion. Hence, rescue groups must transfer patients to pick-up nodes to access available vehicles. However, in the second group, vehicles can easily access the evacuating hospitals and there is no need to use pick-up nodes. r There are two types of patients with different priorities: noncritical and critical patients. r There are different transportation modes including ambulances, helicopters, and vans, each with specific capabilities. Notably, it is assumed that ambulances are only used to transfer critical patients since they are equipped to treat just one critical patient and also the available number of them is significantly less than the available number of vans in the response phase. Moreover,

Fig. 1. Schematic view of a simple example of the considered network.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

r

r r r r

161

it is not possible to transfer critical patients through vans due to their lack of equipment to treat critical patients. There are two choices to transfer patients from pick-up nodes (i.e., for patients of the first group of evacuating hospitals) or evacuating hospitals (i.e., for patients of the second group of evacuating hospitals) to receiving hospitals: (i) Transferring patients directly to receiving hospitals using land path. (ii) Transferring patients to helicopter stations in order to use air path to transfer them to receiving hospitals. It is assumed that each vehicle will be fully loaded before leaving. Capacity limitation is considered for different types of transportation resources. Capacity limitation is considered for receiving hospitals. Since there are limitations in the number of available transportation resources in the response phase, the model has to allocate these resources to the pick-up nodes, helicopter stations, and the second group evacuating hospitals in an efficient way.

In the proposed model, the objective is to evacuate the hospitals in such a way to minimize the total evacuation time and to minimize the total weighted number of unevacuated patients in each period.

2.2. Mathematical programming The following notations are used in the formulation of hospital evacuation problem. It should be noted that the imprecise parameters are denoted with a tilde sign (∼). Indices i i i j v h

index of evacuating hospitals i = {1, . . . , I}, i = (i ∪ i ). index of first group’s evacuating hospitals (i ⊂ i). index of second group’s evacuating hospitals (i ⊂ i). index of receiving hospitals j = {1, . . . , J}. index of pick-up nodes v = {1, . . . , V }. index of helicopter stations h = {1, . . . , H}.  1, noncritical index of patientsk = 2, critical  1, van index of land vehicles f = . 2, ambulance

k f t, s

index of time periods t = {1, . . . , T }, s = {1, . . . , T } (i.e., T: length of the planning horizon).

Parameters teiv j tv

transportation time from evacuating hospital i to vehicle station v.

vj f

transportation time from pick-up node v to receiving hospital j using land vehicle type f.

 tvh vh f  tih

transportation time from pick-up node v to helicopter station h land vehicle type f.

j ti ijf j th

transportation time from evacuating hospital i to receiving hospital j using land vehicle type f.

 NP ki

total number of patients of type k in evacuating hospital i.

ih f

hj

transportation time from evacuating hospital i to helicopter station h using land vehicle type f. transportation time from helicopter station h to receiving hospital j using helicopters.

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

162

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

 T HN

t

total available number of helicopters in response phase at period t.

t

 T VNf Hocap jk V cap Hcap L

total available number of land vehicles of type f in response phase at period t. capacity of receiving hospital j for treating patients of type k. capacity of each van. capacity of each helicopter. number of patients who are transferred together by rescue groups from evacuation hospitals to pick-up nodes. priority for evacuating patients of type k.

Pk

Decision variables xtkivh j number of patients of type k transported from evacuating hospital i to receiving hospital j at time period t through pick-up node v and helicopter station h. number of patients of type k transported from evacuating hospital i to receiving hospital j at time period t ytkiv j through pick-up node v. number of patients of type k transported from evacuating hospital i to receiving hospital j at period t wtkih j through helicopter station h. number of patients of type k transported directly from evacuating hospital i to receiving hospital j at ztki j period t. number of allocated helicopters to helicopter station h at period t. AHht AVvtf

number of allocated land vehicles of type f to pick-up node v at period t.

AIitf SGtiv

number of allocated land vehicles of type f to evacuating hospital i at period t. number of times that rescue groups transfer patients from evacuating hospital i to pick-up node v at period t. number of land vehicles of type f used to transfer patients from pick-up node v to receiving hospital j at time period t. number of land vehicles of type f used to transfer patients from pick-up node v to helicopter station h at period t. number of land vehicles of type f used to transfer patients from evacuating hospital i to helicopter station h at period t. number of helicopters used to transfer patients from helicopter station h to receiving hospital j at period t.

Nv jvt j f Nvhtvh f Nihtih f Nh jht j Ni jit j f

number of land vehicles of type f used to transfer patients from evacuating hospital i to receiving hospital j at period t. number of patients of type k not evacuated till period t.

RCPkt

The biobjective mathematical model is presented as follows: MinOFV1 =

 i

t

+

 t

i

MinOFV2 =

h

t˜ei v SGti v +

v

t˜ihi h f Nihti h f +

t



v

j

i

j

t˜v jv j f Nv jvt j f +

f

 t

f

 t

f

t



t˜i ji j f Ni jit j f +

v

h

 t

h

t˜vhvh f Nvhtvh f

f

t˜h jh j Nh jht j

(1)

j

Pk · RCPkt

k

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

(2)

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

163

s.t.

⎛ ⎞    ⎝ xtki vh j + ytki v j ⎠ ≥ NP˜ ki t

v

v

j

h

j

⎛ ⎞    ⎝ wtki h j + ztki j ⎠ ≥ NP˜ ki t

j

h



j

∀i ⊂ i ∀k

(3)

∀i ⊂ i ∀k

(4)

∀k, t

(5)

∀t

(6)

∀f,t

(7)

⎛ t     ⎝ ˜ ki − NP xski vh j + yski v j + wski h j

i

i

s=1

+

 i

T˜ HN t ≥



v



h

i

j

v

zski j ⎠ = RCPkt

j

i

h

j

j

AHht

h

˜ N tf ≥ TV



AVvtf +

 i

v

 k

xtki vh j +

j

h





∀i ⊂ i

(8)

∀v, t ≤ Hcap × Nh jht j

∀h, j, t

(9)

∀v, h, t

(10)

yt2i v j ≤ Nv jvt j2

∀v, j, t

(11)

wt2i h j ≤ Nihti h2

∀i ⊂ i ∀h, t

(12)

∀i ⊂ i ∀ j, t

(13)

i

k

 i

xtki vh j +

ytki v j ≤ L × SGti v

j

k

 

AIit f

v

wtki h j

i

xt2i vh j ≤ Nvhtvh2

j

 i

 j

zt2i j ≤ Ni jit j2

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

164



M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

∀v, h, t

(14)

yt1i v j ≤ V cap × Nv jvt j1

∀v, j, t

(15)

wt1i h j ≤ V cap × Nihti h1

∀i ⊂ i ∀h, t

(16)

∀i ⊂ i ∀ j, t

(17)

∀k, j

(18)

∀h, t

(19)

Nv jvt j f ≤ AVvtf

∀v, f , t

(20)

Ni jit j f ≤ AIit f

∀i ⊂ i ∀f,t

(21)

i



xt1i vh j ≤ V cap × Nvhtvh1

j

i

 j

zt1i j ≤ V cap × Ni jit j1   t



i

v

Nvhtvh f +

v

i

wtki h j +

h





ztki j

≤ Hocap jk

i

j

h



i

h

ytki v j +



Nh jht j ≤ AHht

j



xtki vh j +



Nihti h f +

h

 j

xtki vh j , ytki v j , wtki h j , ztki j  integer

∀i , i ⊂ i

AHht , AVvtf , AIit f  integer

∀ j, k, v, h

SGtiv , Nh jht j , Nvhtvh f , Nv jvt j f , Nihti h f , Ni jit j f  integer

∀t

RCPkt ≥ 0

∀k, t.

(22)

(23)

The first objective function (1) minimizes the total evacuation time, including travelling time of patients to reach pick-up nodes by rescue groups, travelling time among pick-up nodes to receiving hospitals, travelling time among pick-up nodes to helicopter stations, travelling time among evacuating hospitals to helicopter stations, traveling time among evacuating hospitals to receiving hospitals, and traveling time among helicopter stations to receiving hospitals. The second objective function (2) minimizes the total weighted number of unevacuated patients in each period. Indeed, this objective function aims at evacuating the critical patients sooner than the noncritical ones, since a higher priority is considered for the critical patients. Constraint (3) enforces that all patients (i.e., noncritical and critical) who belong to the first group of evacuating hospitals must be evacuated during the planning horizon. Similarly, constraint (4) enforces that all patients from the second group of evacuating hospitals must be evacuated during the planning horizon. Equation (5) calculates the number of unevacuated patients in each period. Constraint (6) allocates the total number of available helicopters in the response phase to  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

165

the helicopter stations. Constraint (7) allocates the total number of available ambulances and vans in the response phase to the pick-up nodes and the second group of evacuating hospitals. Constraint (8) calculates how many times rescue groups have to transfer patients to pick-up nodes at period t. Constraint (9) calculates the number of helicopters among available helicopters in station h used to transfer patients from helicopter station h to receiving hospital j at period t. Constraints (10)–(13) calculate the number of ambulances used to transfer patients from pick-up node v to helicopter station h, pick-up node v to receiving hospital j, receiving hospital i to helicopter station h, and receiving hospital i to evacuating hospital j, respectively. Similarly, constraints (14)–(17) calculate the number of vans used to transfer patients from pick-up node v to helicopter station h, pick-up node v to receiving hospital j, receiving hospital i to helicopter station h, and receiving hospital i to evacuating hospital j, respectively. Constraint (18) indicates that the total number of patients transferred to receiving hospital j at period t must not exceed the capacity of receiving hospital j at this period. Constraint (19) indicates that the number of used helicopters in station h must not exceed the capacity of total allocated helicopters to station h at time period t. Constraint (20) makes sure that the number of used vans and ambulances in pick-up node v does not exceeded the capacity of total allocated vans and ambulances to pick-up node v at period t. Similarly, constraint (21) makes sure that the number of used vans and ambulances in evacuating hospital i does not exceeded the capacity of total allocated vans and ambulances evacuating hospital i at period t. Constraints (22) and (23) are domain constraints.

3. RPP approach Generally, uncertainty in data can be divided into two main types. The first type is randomness that stems from random nature of parameters. In the literature, when information about distribution functions of the random variables are available, stochastic programming methods are the most applied approaches to deal with randomness (Zhalechian et al., 2016). However, in many real-life situations, there is not enough information about the distribution functions of the random data. Therefore, robust programming methods are used to deal with this sort of uncertainty (Mousazadeh et al., 2015). The second type is epistemic that deals with ill-known and imprecise parameters. Indeed, in this type of uncertainty, there is a lack of knowledge to estimate the exact values of these uncertain parameters. Usually, possibilistic programming approaches are applied to deal with epistemic uncertainty in input data. In our model, the number of patients in evacuating hospitals and the number of available vehicles in the response phase are highly uncertain and depend on many factors such as the nature of disaster and time of impact. Furthermore, transportation times between nodes are another uncertain input data in the model. It is worth noting that there is a lack of knowledge to fit probability distributions for these uncertain parameters. Even, if one could fit probability distributions for these uncertain parameters, it might be possible that these parameters do not behave according to the historical data (Zhalechian et al., 2016). Hence, these parameters are assumed to be imprecise. In order to cope with these possibilistic data in the model, the chance-constraint possibilistic programming approach is usually applied. In the optimization problems, a solution that has both feasibility robustness and optimality robustness is called a robust solution. Feasibility robustness means that the solution should remain feasible for (almost) all possible values of uncertain parameters.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

166

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Fig. 2. Trapezoidal distribution of fuzzy parameters.

Optimality robustness means that the value of objective function for the solution should remain close to optimal value (i.e., there should be a minimum deviation between the optimal value and the value of objective function) for (almost) all possible values of uncertain parameters (Pishvaee et al., 2012). However, the chance-constraint possibilistic programming approach is unable to guarantee the robustness of the obtained solutions. In order to cope with such inefficiency, an RPP approach is used, which is recently proposed by Pishvaee et al. (2012) and Pishvaee and Khalaf (2016). Triangular and trapezoidal distributions are the two most common distributions considered for the imprecise parameters. In this way, we adopt trapezoidal possibility distribution for the imprecise parameters in our model (see Fig. 2). In addition, to cope with the possibilistic constraints, which involve possibilistic data, the possibilistic chance constrained programming (PCCP) approach is usually applied. The possibility measures (Pos) and the necessity measure (Nec) are the basic fuzzy measures introduced in the literature. It should be noted that the Pos and Nec measures indicate the optimistic and pessimistic attitude of decision maker, respectively. Since the considered problem is in the disaster management area, it is more reasonable and conservative to assume that the decision maker has a pessimistic attitude toward satisfying the possibilistic chance constraints. Hence, we use Nec measure to ensure satisfaction of possibilistic chance constraints. Furthermore, the expected value operator is used to model the objective functions. The PCCP counterpart of our compact form of the proposed model can be formulated as follows:  

 Min Z = E A˜ x + E B˜ y (24) s.t.

(25)

˜ ≥ α1 Nec{Cx ≥ D}

(26)

Nec{Ex ≤ F˜ } ≥ α2

(27)

Gx ≤ Hy

(28)

˜ ≥ α3 Nec{Ix = J}

(29)

y ∈ {0, 1} , x ≥ 0 and integer,

(30)

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

167

where A, B,, D, F,, and J are the imprecise parameters and ak represents the minimum confidence level of the kth chance constraint. In this approach, the decision maker determines several initial values for each confidence level subjectively, and then an interactive experiment is done to find a satisfactory level. Based on Dubois and Prade (1987) and Heilpern (1992), if c˜ is a trapezoidal fuzzy number, the expected interval (EI) and expected value (EV) of c˜ can be defined as follows:   

c(1) + c(2) c(3) + c(4) ˜ , E (c) ˜ = ˜ = E (c) , (31) EI [c] 2 2 ˜ = EV [c]

c(1) + c(2) + c(3) + c(4) . 4

In addition, the necessity of c˜ ≤ r can be calculated as follows: ⎧ 1, c(4) ≤ r ⎪ ⎪ ⎪ ⎪ ⎨ r − c (3)   Nec c˜ ≤ r = c − c , c(3) ≤ r ≤ c(4) ⎪ (4) (3) ⎪ ⎪ ⎪ ⎩ 0, c(3) ≥ r It can be shown that (Inuiguchi and Ramık, 2000; Liu and Iwamura, 1998):   Nec c˜ ≤ r ≥ a ⇔ r ≥ (1 − a) c(3) + ac(4) 0 ≤ a ≤ 1.

(32)

(33)

(34)

Notably, fuzzy chance constraints can be converted into their equivalent crisp ones using Equation (34). Hence, the crisp counterpart formulation of the PCCP model can be stated as:     B(1) + B(2) + B(3) + B(4) A(1) + A(2) + A(3) + A(4) x+ y (35) Min Z = 4 4 s.t.   Cx ≥ 1 − α1 D(3) + α1 D(4)

(36)

  Ex ≤ 1 − α2 F(2) + α2 F(1)

(37)

Gx ≤ Hy α   3 Ix ≤ J(3) + 1 − 2  α  3 Ix ≥ J(2) + 1 − 2

(38) α3  J(4) 2 α3  J(1) 2

y ∈ {0, 1} , x ≥ 0 and integer.

(39) (40) (41)

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

168

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Accordingly, the crisp counterpart formulation of PCCP model for the studied problem can be formulated as follows:       teiv(1) + teiv(2) + teiv(3) + teiv(4)

SGti v MinE OFV1 = 4 t v i      tv jv j f (1) + tv jv j f (2) + tv jv j f (3) + tv jv j f (4)  Nv jvt j f + 4 t v j f

+

     tvhvh f (1) + tvhvh f (2) + tvhvh f (3) + tvhvh f (4)  t

+

v

i

h

4

f

     ti ji j f (1) + ti ji j f (2) + ti ji j f (3) + ti ji j f (4)  i

t

+

f

     tihi h f (1) + tihi h f (2) + tihi h f (3) + tihi h f (4)  t

+

h

Nvhtvh f

4

j

f

 t

4

th jh j(1) + th jh j(2) + th jh j(3) + th jh j(4)

Ni jit j f



4

j

h

(42) Nihti h f

Nh jht j

MinOFV2

(43)

s.t. ⎛ ⎞      ⎝ xtki vh j + ytki v j ⎠ ≥ 1 − α1 NPki (3) + α1 NPki (4) t

v

v

j

h

∀i ⊂ i

j

(44)

∀k ⎛ ⎞      ⎝ wtki h j + ztki j ⎠ ≥ 1 − α2 NPki (3) + α2 NPki (4) t

h

j

∀i ⊂ i

j

(45)

∀k ⎛ ⎞ t      ⎝ xskivh j + yskiv j + wski h j + zski j ⎠ s=1



i

v

h

  α3  i

2

j

i

v

j

i

h

j

i

∀k, t

j

  α  NPki(3) + 1 − 3 NPki(4) − RCPkt 2

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

(46)

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189



t      ⎝ xskivh j + yskiv j + wski h j + zski j ⎠ i

s=1



v

  α3  2

i



h

j

i

v

i

j

h

  α  NPki(2) + 1 − 3 NPki(1) − RCPkt 2

j

i

169

⎞ ∀k, t

j

(47)

  t t AHht ≤ 1 − α4 T HN(2) + α4 T HN(1)

∀t

(48)

∀ f , t.

(49)

h



AVvtf +

 i

v

  AIit f ≤ 1 − α5 TV N tf (2) + α5 TV N tf (1)

Several versions of RPP are introduced by Pishvaee et al. (2012). Among them, we use RPP-II, in which the decision maker is only sensitive to the deviation of value of objective functions over the expected optimal value. Finally, the RPP formulation for the studied problem is stated as follows: 

 

 

 MinE Zr + γ Zmax,r − E Zr + δ1r NPki (4) − 1 − α1 NPki (3) − α1 NPki (4) 

  + δ2r NPki (4) − 1 − α2 NPki (3) − α2 NPki (4)     α3  α  + δ3r ∀r NPki(3) + 1 − 3 NPki(4) − NPki(3) 2 2 i α     α  (50) NPki(2) − 3 NPki(2) − 1 − 3 NPki(1) + δ3r 2 2 i    t t t + α4 T HN(1) − T HN(1) + δ4r 1 − α4 T HN(2)    + δ5r 1 − α5 TV N tf (2) + α5 TV N tf (1) − TV N tf (1) s.t. and constraints (8)–(23), (44)–(49). The first term of RPP objective function r aims to minimize its mean value. The second term controls the optimality robustness of the solution. Indeed, it minimizes the maximum possible value of the objective function r and its expected value that can be controlled by constant coefficient γ . The other terms, control the feasibility robustness of the solution. Indeed, they aim to minimize the difference between the worst case value of the imprecise parameters and their chance-constraints value. It should be noted that, in this model αi is a decision variable for all i and δir is a penalty parameter related to term i and OFVr . Also, Zmax,1 is obtained as follows:    Zmax,1 = teiv(4) SGti v + tv jv j f (4) Nv jvt j f + tvhvh f (4) Nvhtvh f t

+

i

v

 t

i

h

f

t

tihi h f (4) Nihti h f +

v

j

f

 t

i

j

f

t

ti ji j f (4) Ni jit j f +

v

h

f

 t

h

th jh j(4) Nh jht j .

j

(51)  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

170

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Notably, the term γ (Zmax − E[Z]) is not considered in OFV2 , since OFV2 does not include imprecise parameters. The proposed RPP approach has two main advantages. First, it benefits from both advantages of possibilistic and robust programming approaches. Second, it reduces the time-consuming experiments to find the appropriate value of constraint’s satisfaction level and assures the optimality of the specified value of satisfaction level (Pishvaee et al., 2012).

3.1. Handling the multiple objective functions In the literature, various approaches have been developed to solve multiobjective programming models, including a priori approach, a posteriori approach, and interactive approach. Fuzzy interactive methods are attractive because of the ability to measure the satisfaction degree of each objective function. This ability can help the decision makers to base the decision on the satisfaction and preference of objective functions (Torabi and Hassini, 2008). In order to solve the biobjective model, this paper proposes a new hybrid approach by combining the presented RPP with TH method, which is proposed by Torabi and Hassini (2008). The steps of the proposed hybrid method are summarized as follows: Step 1: Formulate the biobjective robust possibilistic mathematical model. Step 2: Convert the imprecise objective functions to their crisp forms using the expected value operator. Step 3: Convert the chance constraints into their crisp forms using the necessity measure (Nec). Step 4: Determine the positive ideal solution (PIS) and negative ideal solution (NIS) for each objective function. Step 5: Determine a linear membership function for each objective function as follows: ⎧ 1 if Zr < ZrPIS ⎪ ⎪ ⎪ ⎪ ⎨ ZNIS − Z r r if ZrNIS ≤ Zr ≤ ZrPIS . (52) μr (x) = NIS − Z PIS Z ⎪ r r ⎪ ⎪ ⎪ ⎩ 0 if Zr > ZrNIS Step 6: Convert the biobjective model to an equivalent single-objective one using the TH aggregation function (Torabi and Hassini, 2008), which is calculated by  max λ (x) = ϑλ0 + (1 − ϑ ) ϕr μr (x) (53) r

s.t. λ0 ≤ μr (x)

∀r ∈ {1, 2, 3}

x ∈ F (x) , λ0 , ϑ ∈ [0, 1] ,

(54) (55)

where μr (x) is the satisfaction degree of rth objective function, F (x) indicates the feasible region involving constraints of the equivalent model. Furthermore, ϕr denotes the relative importance of  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189



171

the rth objective function ( r ϕr = 1; ϕr > 0), ϑ is the coefficient of compensation, and λ0 is the minimum satisfaction degree of objectives. Step 7: Specify the importance weight of the objectives (ϕr ) and the coefficient of compensation (ϑ ). Then, solve the single-objective model in order to obtain objective functions. If the decision maker is not satisfied with the current solution, manipulate the parameters to provide another solution that satisfies the decision maker.

4. Solution methods The model presented in Section 2 is a mixed-integer linear programming (MILP) optimization problem. The presented model was coded in GAMS 24.1 software. The reliability and validity of the proposed model was tested on randomly generated data sets on a Pentium 4 computer with 2.6 GHz CPU and 4 GB RAM. However, the required computation time to solve the larger instances is increased considerably for even several days. To overcome this limitation, in this paper we developed two metaheuristic algorithms, namely genetic algorithm (GA) and hybrid imperialistic competitive algorithm (HICA), to find near-optimal solutions in a relatively small amount of time. The GA and HICA are coded in MATLAB 2013b. In the following subsections, a brief description of GA and HICA is provided.

4.1. Genetic algorithm Genetic algorithm (GA) is a nature-inspired metaheuristic that belongs to the class of evolutionary algorithms (EAs). Indeed, this algorithm optimizes problems using techniques inspired by natural evolution including inheritance, crossover, and mutation. GA starts from a population of randomly generated individuals and then this population evolves by iteratively applying selection, crossover, and mutation operators. During selection, individual solutions are selected based on a fitness-based process where the qualities of the represented solutions are measured. In each iteration, except the first iteration, a new population is generated by applying a combination of crossover and mutation ¨ on a number of solutions (Sorensen, 2015). In the crossover operation, a pair of solutions, called “parent,” is selected to produce new solutions, called “offspring.” Unlike the crossover operator, the mutation operator applies on only one solution. This operator prevents the algorithm from trapping in local optima by exploring the solution space. Figure 3 shows the steps of GA.

4.2. Imperialist competitive algorithm Imperialist competitive algorithm (ICA) is an EA, which uses sociopolitical evolution of human as a source of inspiration for developing a strong optimization strategy (Atashpaz-Gargari and Lucas, 2007). This algorithm starts with an initial population that is called country. The countries include two categories: imperialists and colonies. Imperialists are the best countries and colonies are the countries that belong to the imperialists respective to their power. Thus, imperialists and their colonies together create the initial empires. In order to distribute the colonies among empires and  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

172

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Fig. 3. Pseudo code of developed GA.

Fig. 4. Moving a colony toward its relevant imperialist with random angle.

begin the competition between imperialists, the total power of an empire that includes the power of the imperialists and the power of its colonies needs to be calculated. Hence, the total power of an empire is calculated as the sum of power of imperialist country and percentage of mean power of its colonies. In the next step, colonies start moving toward their relevant imperialist (Assimilation). Figure 4 depicts the assimilation operator in detail. Afterwards, a competition begins among all the imperialists. Any empire that is not able to succeed in this competition and increase its power will be eliminated from the competition. Therefore, weaker empires will gradually lose their power and finally they will collapse. Finally, during the competition among empires, all the countries will  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

173

Fig. 5. Pseudo code of developed HICA.

converge to a stage where just one empire exists and all the other countries become colonies of that empire. As can be seen in Fig. 4, in order to move colonies toward imperialists in other direction, a θ angle is applied with uniform distribution. It should be noted that d is the distance between the colony and the imperialist.

4.2.1. Hybrid imperialist competitive algorithm In the proposed hybrid imperialist competitive algorithm (HICA), two extra operators are added to the ICA. These operators include: crossover and revolution. Figure 5 shows the steps of our proposed HICA.

Crossover. In order to improve the colonies, they share their information using a crossover operator. In this operator, we adopt arithmetic crossover for the continuous part and binary crossover for the binary part to share the information between colonies.

Revolution. Revolution operator is used to prevent the algorithm from getting stuck at a local minimum. In the proposed HICA, a number of colonies are randomly selected to undergo the revolution operator. Inversion method is adopted for the continuous part and reversion method is adopted for the binary part.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

174

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

4.3. Tuning the parameters In order to find an optimal combination of the parameters in the metaheuristic algorithms, parameter adjustment must be done. There are three major methods in the literature to adjust the parameters of metaheuristic algorithms: (i) fractional factorial design, (ii) response surface methodology (RSM), and (iii) Taguchi design of experiments, which is more commonly employed among other methods. In this paper, the Taguchi method is applied to tune the effective parameters of the proposed GA and HICA. 4.3.1. Taguchi method Procedure of Taguchi method, which is a special case of fractional factorial design, can be summarized by the following steps (Cheng and Chang, 2007): Step 1: Determining the factors that affect the response variable. Step 2: Determining the levels of each effective parameter during a trial-and-error procedure. Step 3: Selecting the appropriate orthogonal array. Step 4: Conducting experiments based on the orthogonal array distribution. Step 5: Analyzing the results to determine the optimal level of each factor based on the signalto-noise ratio (S/N) ratio method or analysis of variance (ANOVA). As mentioned, there are two approaches to analyze the results of Taguchi method, the ANOVA and the S/N ratio method. The first method is usually used for the experiments that are repeated once and the second method is used for the experiments that have multiple runs (Sadeghi et al., 2014). In this paper, analyzing the results of Taguchi method is done by the S/N ratio method. S/N ratio method uses the S/N as the response variable, where S represents controllable factors and N stands for noise factors that affect the response (Roy, 1990).

5. Computational experiments In this section, at first the robustness and desirability of the obtained solutions using RPP model in comparison to PCCP model are assessed. Then, a sensitivity analysis has been done on parameter of TH method to obtain solutions that satisfy the decision maker. Afterwards, to validate the proposed metaheuristic algorithms for a number of small- and medium-sized instances, exact solutions of GAMS software and near-optimal solutions of metaheuristic algorithms are compared. Finally, the performance of HICA is evaluated by the performance of GA over a number of large-sized instances.

5.1. RPP model evaluation To evaluate the robustness and desirability of the obtained solutions using RPP model, we have considered an instance and generated 10 random values for imprecise parameters (i.e., distribution  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

175

Table 1 Parameters of test problem Parameters teiv tv jv j f tvhvh f tihih f ti ji j f th jh j L Pk V cap

Uniform∼ Uniform∼ Uniform∼ Uniform∼ Uniform∼ Uniform∼ 8 [0.3, 0.7] 6

(5, 10) (30, 40) (5, 10) (5, 10) (35, 45) (5, 10)

Hcap NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

5 [23, 24, 16, 25, 21] [10, 11, 13, 10, 8] [15, 15] [25, 30] [5, 5] [45, 50, 30] [30, 30, 35]

functions of parameters follow Table 1). Afterward, the performance of the obtained solutions is tested by replacing the obtained solutions (x∗ , y∗ ) in the compact form of the model presented below: 

Min Z = Arand x∗ + Brand y∗ + δ1 R + δ2 R

(56)

s.t. Cx ≤ D

(57)

Ex ≤ F y

(58)

Gx ≤ Hrand + R

(59)



Gx + R ≥ Hrand 

R, R ≥ 0,

(60) (61)

where R and R are the only decision variables in this model that determine the violation of chance constraints. Moreover, the parameters δ1 and δ2 are the penalty values that can be estimated similar to the penalty values of the RPP model. Notably, the average and standard deviation of objective function values (OFVs) are used as performance measures to evaluate the RPP model (for more detailed information, see Pishvaee et al., 2012). The results related to the evaluation of our RPP model in comparison to PCCP model are represented in Tables 2 and 3. The results show that the RPP model outperforms the PCCP model under different values for ai (i.e., PCCP 0.7 and PCCP 0.9) both in terms of optimality (see Fig. 6) and robustness (see Fig. 7). Since the obtained solutions of RPP model are less than the obtained solutions of PCCP models both in average value and standard deviation, the results justify the application of RPP model to cope with uncertainty.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

176

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Table 2 Performance of PCCP model under αi = 0.7 and αi = 0.9 PCCP ai = 0.7

Problem size

PCCP ai = 0.9

|I| × |V | × |H| × |J| × |T |

Number

OFV1

OFV2

OFV1

OFV2

5×6×2×3×2

1 2 3 4 5

3153.60 3144.96 3291.15 3313.84 3090.87

217 166 194 205 194

3259.80 3186.82 3376.25 3490.48 3308.64

247 172 236 281 211

Average St. Dev.

3198.88 79.95

195.20 15.42

3324.40 94.53

229.40 33.29

Table 3 Performance of RPP model Problem size

RPP

|I| × |V | × |H| × |J| × |T |

Number

OFV1

OFV2

5×6×2×3×2

1 2 3 4 5

3076.20 3063.06 3183.85 3021.28 2966.43

194 153 172 187 179

Average St. Dev.

3062.16 65.63

177 12.87

Fig. 6. Average values of objective function values.

5.2. Sensitivity analyses In this section, several sensitivity analyses are performed in order to validate the presented model and investigate the effect of data parameters on the test cases. The results are discussed in depth to provide some managerial implications.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

177

Fig. 7. Standard deviation values of objective function values. Table 4 Sensitivity analysis results for the test problem Problem size

Satisfaction degree

Objective function

|I| × |V | × |H| × |J| × |T |

ϑ

μ1

μ2

OFV1

OFV2

5×6×2×3×2

0.9 0.8 0.7 0.6 0.5–0.1

0.610 0.607 0.603 0.603 0.594

0.631 0.654 0.687 0.718 0.749

3428.532 3402.126 3388.923 3388.923 3362.517

319 293 266 257 243

In converting a biobjective problem to a single-objective one using TH method, two parameters: including the relative importance of the kth objective function, and the coefficient of compensation, have to be specified. Relative importance of objective functions (i.e., weight vector) is determined based on the preferences of the decision maker. In our problem, based on the expert’s opinion, we set the weight vector of objectives as ϕ = (0.6, 0.4). Coefficient of compensation determines the balancing amount of a compromise solution. Higher values for coefficient of compensation lead to more balanced compromise solutions. On the other hand, the lower values for it lead to more unbalanced compromise solutions. Therefore, in order to obtain a solution that satisfies the decision maker, a sensitivity analysis is conducted on the coefficient of compensation. Table 4 shows the results of the sensitivity analysis on the considered test problem. As shown in Table 4, increasing the value of ϑ resulted in higher value for the minimum of membership functions and lower value for the weighted average of membership functions. Figures 8 and 9 show a better explanation of the effect of ϑ on the minimum membership functions and the weighted average membership functions, respectively. In order to investigate the effect of the total available number of vehicles and priority for evacuating patients, we considered a test problem. Table 5 shows the parameters of the considered the test problem; other parameters follow Table 1. Figure 10 shows the effect of the total available number of vehicles (i.e., helicopter, van, and ambulance) on the total evacuation time. Generally, increasing the number of vehicles will decrease the total evacuation time. However, the slope of the blue line, which is related to an increase  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

178

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Fig. 8. Min μ_r (x) versus ϑ.

Fig. 9. Weighted average versus ϑ. Table 5 Parameters of the considered test problem No. |I| × |V | × |H| × |J| × |T |

Parameters

1

2

3

4

5

6

7

7×3×3×4×2

NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

40 9 35 35 8 90 35

35 10 35 35 8 105 40

48 8 − − − 85 40

50 5 − − − 120 50

60 10 − − − − −

37 10 − − − − −

53 12 − − − − −

in T HN t , is significantly greater than the slope of the red and gray lines. Since helicopters are faster than other vehicles (i.e., vans and ambulances), the number of available helicopters has a more significant role in reducing the total evacuation time. Indeed, when the number of available helicopters is increased, more critical and noncritical patients can be transferred by helicopters that decrease the total evacuation time.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

179

Fig. 10. Number of vehicles versus total evacuation time.

Fig. 11. Alteration of Pk versus number of unevacuated patients.

Figure 11 illustrates the sensitivity analysis on different values of Pk (i.e., priority for evacuating patients). In this sensitivity analysis, we considered six scenarios, which are reported in Table 6. As expected, increasing the priority of critical-patients forces the model to evacuate more critical patients during the first period. However, considering a very high priority for critical patients is not always the best policy in the hospital evacuation because it may postpone the evacuation of a large number of noncritical patients. As a consequence, decision makers should pay special attention to select proper priorities for the patients.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

180

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Table 6 Different scenarios Scenario

1

2

3

4

5

6

Priority of noncritical patients Priority of critical patients

0.2 0.8

0.3 0.7

0.4 0.6

0.5 0.5

0.6 0.4

0.7 0.3

Table 7 Levels of GA parameters Level Parameter

1

2

3

Npop It Pc Pm Ts

100 200 0.5 0.3 7

200 350 0.65 0.45 9

300 500 0.8 0.6 11

Table 8 Levels of HICA parameters Level Parameter

1

2

3

Npop It Nemp Rp

100 300 5 0.5

200 450 7 0.65

300 600 9 0.75

5.3. Tuning the parameters In this section, parameters of the two proposed metaheuristics are tuned using Taguchi method. In Tables 5 and 6, effective parameters in GA and HICA are reported, respectively. There are four effective factors in GA, including population size(Npop ), iteration number (It), crossover percentage(Pc ), mutation percentage (Pm ), and tournament size (Ts). Effective factors in the HICA include population size(Npop ), iteration number (It), empires number (Nemp ), and revolution probability(R p ). By a trial-and-error procedure, three levels for each factor are considered for both algorithms in Tables 7 and 8. Based on the parameters of the test problem shown in Table 1, Taguchi design of experiments along with the results of experiments for GA and HICA are shown in Tables 9 and 10, respectively. Due to the stochastic nature of the algorithms, five replications are performed in each experiment to achieve more reliable results. In order to analyze the results of Taguchi design of experiment, it is better to convert the values by one of the existing methods in the literature such as relative deviation index (RDI) and relative  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

181

Table 9 Taguchi design of experiments and summary of results for GA Factor levels

Replication

Experiment number

Npop

It

Pc

Pm

Ts

1

2

3

4

5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

100 100 100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 200 300 300 300 300 300 300 300 300 300

200 200 200 350 350 350 500 500 500 200 200 200 350 350 350 500 500 500 200 200 200 350 350 350 500 500 500

0.50 0.50 0.50 0.65 0.65 0.65 0.80 0.80 0.80 0.65 0.65 0.65 0.80 0.80 0.80 0.50 0.50 0.50 0.80 0.80 0.80 0.50 0.50 0.50 0.65 0.65 0.65

0.30 0.30 0.30 0.45 0.45 0.45 0.60 0.60 0.60 0.60 0.60 0.60 0.30 0.30 0.30 0.45 0.45 0.45 0.45 0.45 0.45 0.60 0.60 0.60 0.30 0.30 0.30

7 9 11 7 9 11 7 9 11 7 9 11 7 9 11 7 9 11 7 9 11 7 9 11 7 9 11

0.5988 0.6068 0.5996 0.6153 0.6146 0.6092 0.6190 0.6148 0.6137 0.6090 0.6135 0.6172 0.6220 0.6144 0.6223 0.6163 0.6159 0.6277 0.6243 0.6172 0.6131 0.6177 0.6300 0.6303 0.6173 0.6266 0.6142

0.6032 0.6128 0.6041 0.6166 0.6192 0.6207 0.6156 0.6228 0.6170 0.6176 0.6087 0.6080 0.6137 0.6189 0.6095 0.6146 0.6160 0.6144 0.6148 0.6159 0.6095 0.6165 0.6302 0.6176 0.6273 0.6216 0.6301

0.603 0.6082 0.6016 0.6144 0.6131 0.6147 0.6144 0.613 0.6165 0.6094 0.6107 0.6304 0.6194 0.613 0.6081 0.6173 0.6208 0.618 0.6106 0.6124 0.6101 0.615 0.6171 0.6256 0.6207 0.6131 0.6157

0.5975 0.6015 0.6052 0.6233 0.6196 0.6097 0.6172 0.6173 0.6141 0.6188 0.6173 0.6042 0.6174 0.6221 0.614 0.6266 0.629 0.6176 0.6248 0.6159 0.6131 0.6155 0.6247 0.6111 0.6122 0.617 0.627

0.5963 0.6029 0.6105 0.5467 0.6139 0.6143 0.6161 0.6245 0.6248 0.6104 0.6267 0.6125 0.6139 0.6181 0.6103 0.6254 0.6124 0.6282 0.6145 0.6101 0.6071 0.626 0.6334 0.6294 0.6209 0.6122 0.628

Table 10 Taguchi design of experiments and summary of results for HICA Factor levels

Replication

Experiment number

Npop

It

Pc

Pm

1

2

3

4

5

1 2 3 4 5 6 7 8 9

100 100 100 200 200 200 300 300 300

300 450 600 300 450 600 300 450 600

5 7 9 7 9 5 9 5 7

0.50 0.65 0.75 0.75 0.50 0.65 0.65 0.75 0.50

0.5806 0.5965 0.5962 0.5931 0.5966 0.5955 0.5978 0.6063 0.5960

0.5989 0.6066 0.5914 0.5922 0.5957 0.5951 0.6026 0.6020 0.6103

0.5951 0.5947 0.5896 0.5889 0.6059 0.6057 0.5979 0.6031 0.6007

0.5898 0.5909 0.5929 0.5960 0.5983 0.5980 0.5959 0.5995 0.6174

0.5946 0.5794 0.5914 0.5982 0.5905 0.6093 0.5810 0.6115 0.6115

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

182

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Fig. 12. The mean S/N plot for GA parameters.

Fig. 13. The mean S/N plot for HICA parameters.

percentage deviation (RPD). In this section, we utilized RPD index to convert the objective values. The RPD index formula is as follows:   Maxsol − Algsol × 100, (62) RDI = Maxsol where Algsol is the obtained value of the metaheuristics for each experiment. Also, Maxsol denotes the best solution. The mean S/N plots for each level of GA and HICA parameters are shown in Figs. 12 and 13, respectively. Table 11 presents the tuned parameters of two proposed algorithms.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

183

Table 11 Tuned parameters of factors in GA and HICA Parameter GA

Value

Parameters ICA

Value

Npop It Pc Pm Ts

300 500 0.5 0.6 9

Npop It Nemp Rp

300 600 7 0.5

5.4. Comparative study In this section, the performance of the developed metaheuristic algorithms is evaluated in a comparative study over a number of randomly generated data sets. Distribution functions of parameters follow Table 1. Other parameters of the data sets are represented in Table 12. Altogether, we considered 18 test problems ranging from small- to large-sized problems. In Tables 13 and 14, optimal solutions obtained from GAMS software and the developed metaheuristics (GA and HICA) are compared for small- and medium-sized problems. The percentage of relative gaps for the five replications of each test problem and computation time are averaged and recorded as the average percentage of relative gap. Since large-sized test problems cannot be solved by GAMS software, only the comparison between GA and HICA is reported in Table 13. Hence, for small- and medium-sized test problems, the GA and HICA are compared with the optimal solutions obtained from GAMS as [100 × (GGAMS − GMeta )/GMeta ], where GGAMS and GMeta are the TH objective function value obtained from GAMS software and solution obtained from metaheuristic algorithm, respectively. Similarly, for large-sized test problems, a gap between GA and HICA is calculated as [100 × (SHICA − SGA )/SGA ], where SHICA and SGA are the solutions obtained from HICA and GA metaheuristics, respectively. It should be noted that in all the test problems HICA performs better than GA as shown in Table 15. Figure 14 shows the average percentage of relative gap of GA and HICA. Furthermore, the average percentage of relative gap of GA in comparison to HICA is shown in Fig. 15. Figure 16 shows the computation time of GA and HICA. The results show that the developed HICA and GA can find near-optimal solutions, since there is a small average percentage of relative gap (i.e., less the 4%) between the optimal solution and the obtained solutions from HICA and GA. As expected, the computation time of developed metaheuristic algorithms are increased by increasing the size of test problems. Moreover, it can be seen that the number of evacuating hospitals and the length of the planning horizon are two important parameters that significantly make the model more complex and increase the computation time. Furthermore, increasing the value of these parameters will increase the average percentage of relative gap between HICA and GA algorithms. Based on Table 15 and Fig. 15, it can be concluded that HICA can perform better than GA at solving large-sized test problems.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

184

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

Table 12 Small- and medium-sized data sets 1–6 No. Data set

|I| × |V | × |H| × |J| × |T |

Parameters

1

3×3×2×2×1

NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

8 2 11 20 8 10 10

4 8 − − − 15 10

4 4 − − − − −

− − − − − − −

− − − − − − −

− − − − − − −

− − − − − − −

2

5×4×2×3×1

NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

23 10 20 45 10 45 30

24 5 − − − 50 30

16 8 − − − 30 35

25 10 − − − − −

21 6 − − − − −

− − − − − − −

− − − − − − −

3

5×6×3×4×2

NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

23 10 8 20 4 35 20

24 9 12 30 7 35 25

16 8 − − − 45 25

25 10 − − − 40 35

21 10 − − − − −

− − − − − − −

− − − − − − −

4

7×5×3×4×2

NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

40 9 20 30 8 90 35

35 10 20 35 8 105 40

48 8 − − − 85 40

50 5 − − − 120 50

60 10 − − − − −

37 10 − − − − −

53 12 − − − − −

5

7×6×3×4×2

NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

40 12 30 35 10 80 30

35 5 35 35 10 90 35

48 7 − − − 120 30

50 10 − − − 100 45

60 12 − − − − −

37 10 − − − − −

53 8 − − − − −

6

7×7×3×5×2

NP1,i NP2,i TV N1t TV N2t T HN t Hocap j1 Hocap j2

90 10 75 50 15 210 100

105 12 80 50 20 150 70

87 10 − − − 150 80

140 15 − − − 250 65

102 12 − − − 200 95

110 12 − − − − −

130 15 − − − − −

1

2

3

4

5

6

7

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

185

Table 13 Average percentage of relative gaps and CPU time in GA for small- and medium-sized test problems GAP Replications Data set

I

V

H

J

T

1

2

3

4

5

Average gap

Average time

1 2 3 4 5 6

3 5 5 7 7 7

3 4 6 5 6 7

2 2 3 3 3 3

2 3 4 4 4 5

1 1 2 2 2 2

0.055 0.360 0.805 1.140 1.650 3.051

0.066 0.168 0.690 1.330 1.540 3.164

0.055 0.300 0.104 1.045 1.650 3.051

0.077 0.240 0.805 0.950 1.540 3.164

0.066 0.264 0.104 1.235 1.870 3.277

0.064 0.266 0.501 1.140 1.650 3.141

63 65 65 90 92 98

Table 14 Average percentage of relative gaps and CPU time in HICA for small- and medium-sized test problems GAP Replications Data set

I

V

H

J

T

1

2

3

4

5

Average gap

Average time

1 2 3 4 5 6

3 5 5 7 7 7

3 4 6 5 6 7

2 2 3 3 3 3

2 3 4 4 4 5

1 1 2 2 2 2

0.044 0.060 0.081 0.751 1.210 2.147

0.033 0.036 0.058 0.086 1.430 2.373

0.055 0.048 0.069 0.950 1.320 2.260

0.044 0.048 0.069 0.475 1.100 2.486

0.033 0.060 0.081 0.618 1.320 2.396

0.042 0.050 0.071 0.576 1.276 2.332

75 77 83 98 103 114

Table 15 Average percentage of relative gaps and CPU time for large-sized test problems HICA

GA Replications

Data set

I

V

H

J

T

Time

1

2

3

4

5

Gap

Time

7 8 9 10 11 12 13 14 15 16 17 18

15 15 15 15 25 25 25 25 40 40 40 40

10 15 16 18 15 20 25 28 35 40 45 50

5 5 5 4 10 10 8 8 15 15 10 9

9 8 6 6 20 15 15 12 30 25 25 20

3 3 3 3 5 5 5 5 7 7 7 7

158 160 161 158 488 499 491 502 954 953 956 972

3.630 4.200 4.140 3.515 4.510 4.859 5.170 4.935 8.645 10.340 9.888 10.094

3.740 4.200 4.025 3.420 4.290 5.085 5.060 5.145 8.645 10.230 9.991 10.094

3.740 3.960 4.255 3.420 4.620 4.972 5.170 5.145 8.835 10.450 9.785 9.991

3.520 4.080 3.910 3.420 4.620 5.198 4.950 5.040 8.930 10.450 9.888 9.991

3.850 3.960 4.140 3.515 4.180 4.972 5.060 4.935 8.835 10.340 9.991 10.094

3.696 4.080 4.094 3.458 4.444 5.017 5.082 5.040 8.778 10.362 9.909 10.053

190 193 187 197 594 597 590 609 1089 1118 1110 1123

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

186

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189 3.5 GA

GAP (%)

3

HICA

2.5 2 1.5 1 0.5 0

1

2

3

4

5

Test problem No.

6

Fig. 14. The average percentage of relative gap of GA and HICA.

12

GAP (%)

10 8 6 4 2 0

7

8

9

10

11

12

13

14

15

16

17

18

Test problem No.

Fig. 15. The average percentage of relative gap of GA in comparison to HICA.

1200

CPU time

1000

GA HICA

800 600 400 200 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Test problem No.

Fig. 16. The computational time of GA and HICA.

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

187

Table 16 Results of the paired t-test Paired t-test 95% confidence interval of the difference Samples

Pooled St. Dev.

Lower

Upper

t-Value

p-Value

DF

GA-HICA

378.84

−190

295

0.44

0.666

38

In order to see whether a significant difference exists between the CPU times of GA and HICA, a paired t-test is conducted. Paired t-test statistics are listed in Equation (63): ¯ D t= √ , SD / n

 ¯ = where D

n

¯ D

! ! ¯ 2 (Di − D) " and SD = . n−1

(63)

The paired t-test is conducted in Minitab 16. Table 16 shows the results of the test. According to Table 16, p-value is greater than 0.05, therefore we would fail to reject the null hypothesis and it means that the CPU times of GA and HICA are not significantly different from each other. Hence, this similarity in CPU time and outperformance of HICA comparing to GA in terms of objective function values prove that HICA performs better than GA at solving the presented hospital evacuation mathematical model. 6. Conclusions This paper proposed a new biobjective model for hospital evacuation problem under uncertainty in order to minimize total evacuation time and the number of unevacuated patients in each period. In the proposed model, patients are clustered into two groups, including critical patients and noncritical patients. In addition, the evacuating hospitals were divided into two groups. In the first group, it is not possible to send vehicles to the evacuating hospitals due to the poor road condition or congestion. However, in the second group, vehicles can easily access the evacuating hospitals. Furthermore, different transportation modes including ambulances, helicopters, and vans, each with specific capabilities were considered. In order to cope with the uncertainties and to solve the biobjective programming model, an efficient hybrid approach was proposed. In addition, two metaheuristic algorithms (i.e., GA and HICA) were developed to solve large-sized problems. The presented model can provide a decision framework for top managers who can help them handle the difficult process of hospital evacuation. The computational experiments indicate that using aerial vehicles (i.e., helicopters) to transfer patients can significantly reduce the total evacuation time. Hence, it is of crucial importance to provide enough helicopters in the response phase of an evacuation. Moreover, the experiments showed the importance of selecting proper priorities for patients. The analysis further assessed the robustness and desirability of the obtained solutions using RPP model, and the results justified the application of RPP approach to cope with uncertain parameters in the presented hospital evacuation mathematical model. Finally, the obtained results showed that the developed algorithms are capable of solving real-world (i.e., large-sized) problems.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

188

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

References Atashpaz-Gargari, E., Lucas, C., 2007. Imperialist competitive algorithm: an algorithm for optimization inspired by imperialistic competition. Paper presented at the IEEE Congress on Evolutionary Computation, 2007 (CEC 2007), Singapore. ¨ Barbaroso˘glu, G., Ozdamar, L., Cevik, A., 2002. An interactive approach for hierarchical analysis of helicopter logistics in disaster relief operations. European Journal of Operational Research 140, 1, 118–133. Bish, D.R., Agca, E., Glick, R., 2014. Decision support for hospital evacuation and emergency response. Annals of Operations Research 221, 1, 89–106. Chang, B.P., Winton, K., 2001. Analysis of traffic situation in disaster area after the Great Hanshin-Awaji earthquake by using aerial photos. Proceedings of Civil Engineering Project Study 18, 849–858. Cheng, B.-W., Chang, C.-L., 2007. A study on flowshop scheduling problem combining Taguchi experimental design and genetic algorithm. Expert Systems with Applications 32, 2, 415–421. Childers, A.K., Mayorga, M.E., Taaffe, K.M., 2014. Prioritization strategies for patient evacuations. Health Care Management Science 17, 1, 77–87. Childers, A., Taaffe, K., 2010. Healthcare facility evacuations: lessons learned, research activity, and the need for engineering contributions. Journal of Healthcare Engineering 1, 1, 125–140. Childers, A.K., Visagamurthy, G., Taaffe, K., 2009. Prioritizing patients for evacuation from a health-care facility. Transportation Research Record: Journal of the Transportation Research Board 2137, 1, 38–45. Dubois, D., Prade, H., 1987. The mean value of a fuzzy number. Fuzzy Sets and Systems 24, 3, 279–300. Haghani, A., Oh, S.-C., 1996. Formulation and solution of a multi-commodity, multi-modal network flow model for disaster relief operations. Transportation Research Part A: Policy and Practice 30, 3, 231–250. Heilpern, S., 1992. The expected value of a fuzzy number. Fuzzy Sets and Systems 47, 1, 81–86. Inuiguchi, M., Ramık, J., 2000. Possibilistic linear programming: a brief review of fuzzy mathematical programming and a comparison with stochastic programming in portfolio selection problem. Fuzzy Sets and Systems 111, 1, 3–28. Liu, B., Iwamura, K., 1998. Chance constrained programming with fuzzy parameters. Fuzzy Sets and Systems 94, 2, 227–237. McGlown, K.J., 1999. Determinants of the evacuation of health care facilities. PhD dissertation, University of Alabama at Birmingham. McGlown, K.J., 2001. Evacuation of health care facilities: a new twist to a classic model. Natural Hazards Review 2, 2, 90–99. Mousazadeh, M., Torabi, S.A, Zahiri, B., 2015. A robust possibilistic programming approach for pharmaceutical supply chain network design. Computers & Chemical Engineering 82, 115–128. Najafi, M., Eshghi, K., Dullaert, W., 2013. A multi-objective robust optimization model for logistics planning in the earthquake response phase. Transportation Research Part E: Logistics and Transportation Review 49, 1, 217–249. Ozdamar, L., 2011. Planning helicopter logistics in disaster relief. OR Spectrum 33, 3, 655–672. Oh, S.-C., Haghani, S., 1997. Testing and evaluation of a multi-commodity multi-modal network flow model for disaster relief management. Journal of Advanced Transportation 31, 3, 249–282. Pishvaee, M.S., Khalaf, M.F., 2016. Novel robust fuzzy mathematical programming methods. Applied Mathematical Modelling 40, 1, 407–418. Pishvaee, M., Razmi, J., Torabi, S.A., 2012. Robust possibilistic programming for socially responsible supply chain network design: a new approach. Fuzzy Sets and Systems 206, 1–20. Roy, R.K., 1990. A Primer on the Taguchi Method. Van Nostrand Reinhold, New York. Sadeghi, J., Sadeghi, S., Niaki, S.T.A., 2014. Optimizing a hybrid vendor-managed inventory and transportation problem with fuzzy demand: an improved particle swarm optimization algorithm. Information Sciences 272, 126–144. ¨ Sorensen, K., 2015. Metaheuristics—the metaphor exposed. International Transactions in Operational Research 22, 1, 3–18. Tayfur, E., Taaffe, K., 2009a. A model for allocating resources during hospital evacuations. Computers & Industrial Engineering 57, 4, 1313–1323. Tayfur, E., Taaffe, K., 2009b. Simulating hospital evacuation—the influence of traffic and evacuation time windows. Journal of Simulation 3, 4, 220–234.  C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research 

M. Rabbani et al. / Intl. Trans. in Op. Res. 25 (2018) 157–189

189

Torabi, S.A., Hassini, E., 2008. An interactive possibilistic programming approach for multiple objective supply chain master planning. Fuzzy Sets and Systems 159, 2, 193–214. Vogt, B., 1991. Issues in nursing home evacuations. International Journal of Mass Emergencies and Disasters 9, 2, 247–265. White, R.A., Blumenberg, E., Brown, K.A., Contestabile, J.M., Haghani, A., Howitt, A.M., Lambert, T.C., Morrow, B.H., Setzer, M.H., Stanley, E.M., Vel´asquez III, A., 2008. The Role of Transit in Emergency Evacuation. Transportation Research Board, Washington, DC. Zhalechian, M., Tavakkoli-Moghaddam, R., Zahiri, B., Mohammadi, M., 2016. Sustainable design of a closed-loop location-routing-inventory supply chain network under mixed uncertainty. Transportation Research Part E: Logistics and Transportation Review 89, 182–214.

 C 2016 The Authors. C 2016 International Federation of Operational Research Societies International Transactions in Operational Research