An emergency logistics distribution approach for quick response to urgent relief demand in disasters

(1)

An emergency logistics distribution approach for quick

response to urgent relief demand in disasters

Jiuh-Biing Sheu

*

Institute of Traﬃc and Transportation, National Chiao Tung University, 4F, 114 Chung Hsiao W. Rd., Sec. 1, Taipei 10012, Taiwan, ROC

Received 24 October 2005; received in revised form 18 April 2006; accepted 27 April 2006

Abstract

Quick response to the urgent relief needs right after natural disasters through efficient emergency logistics distribution is vital to the alleviation of disaster impact in the affected areas, which remains challenging in the field of logistics and related study areas. This paper presents a hybrid fuzzy clustering-optimization approach to the operation of emergency logistics co-distribution responding to the urgent relief demands in the crucial rescue period. Based on a proposed three-layer emer-gency logistics co-distribution conceptual framework, the proposed methodology involves two recursive mechanisms: (1) disaster-affected area grouping, and (2) relief co-distribution. Numerical studies with a real large-scale earthquake disaster occurring in Taiwan are conducted, and the corresponding results indicate the applicability of the proposed method and its potential advantages. We hope that this study can not only make the proposed emergency logistics system available with more benefits to the development of emergency logistics systems for the urgent needs of disaster areas around the world but also stimulate more excellent researches concerning emergency logistics management.

Keywords: Logistics; Emergency logistics distribution; Relief demand forecasting; Fuzzy clustering; Multi-objective optimization

1. Introduction

Quick response to the urgent need of relief in affected areas right after disasters is a critical issue for emer-gency logistics, which has aroused growing concerns and research interests recently due to the occurrence of several notorious disasters either natural or man-made (e.g., the recent tsunami disaster in the southeastern Asia and the 911 event in the US). Different from general business logistics, emergency logistics is unique in four aspects that may increase the relative complexity and difficulty in solving the induced relief logistics problems, particularly in terms of emergency logistics distribution. First, the demand-related information, e.g., the severity of disaster-induced effects on the affected areas and casualties, is quite limited in the initial search-and-rescue period, and is intuitively unpredictable using historical data. Second, relative to general

* _{Tel.: +886 2 2349 4963; fax: +886 2 2349 4953.} E-mail address:[email protected]

(2)

business logistics management, the emergency logistics resources and the corresponding requirements may not be fully controllable to decision makers in the supply side, adding more challenging issues to the quick-respon-sive emergency logistics distribution system. Third, the damaged infrastructure aﬀected by disasters may incur more unexpected risks of relief distribution, leading to the induced issue of restructuring the emergency logis-tics networks, which should be resolved in a limited time frame. Fourth, regarding the global relief supply for large-scale natural disasters such as the tremendous tsunami casualty of the southeastern Asia, the resulting international relief and logistics resource management problems can make the entire emergency logistics system more complicated, leading to more serious supply–demand imbalance problems in the process of emergency logistics distribution.

Despite the critical importance of emergency logistics distribution, relative to business logistics and supply chain management, only a limited amount of related research has been carried out. Several pioneering studies are illustrated below for further discussion.

The signiﬁcance of issues on relief supply to areas suﬀering from disasters, e.g., drought and earthquakes, and the resulting logistics problems had been addressed previously (Kembell-Cook and Stephenson, 1984;

Knott, 1987; Ardekani and Hobeika, 1988; Long and Wood, 1995), followed by the emergence of diverse

lin-ear programming based models proposed for emergency logistics planning (Knott, 1988; Rathi et al., 1992;

Brown and Vassiliou, 1993 Fiedrich et al., 2000; Barbarosoglu et al., 2002; Ozdarmar et al., 2004). Therein,

a number of researchers tended to formulate the resulting relief transportation issues as multi-commodity multi-modal ﬂow problems with time windows (Rathi et al., 1992; Haghani and Oh, 1996).

By incorporating knowledge-based rules into a linear programming model, the issue of vehicle scheduling for supplying bulk relief of food to a disaster area has been addressed inKnott (1988). Recently,Brown and

Vassiliou (1993)developed a sophisticated real-time decision support system using optimization approaches,

simulation techniques as well as the decision maker’s judgment for both relief resource allocation and assign-ment following a disaster. Considering the multi-commodity supply problems under emergency conditions, three liner programming formulations are proposed inRathi et al. (1992), where the routes and the supply amount carried on each route are assumed to be known in each of the given origin–destination (O–D) pairs. Their purpose, in reality, is to assign a limited number of vehicles loading multiple types of goods in given pairs of origins and destinations such that the induced multi-commodity flow problem is solved with minimal penalties caused by delivery inefficiency, e.g., early and late delivery as well as shipping on non-preferred vehi-cles. Based on the concept of a proposed time-space network,Haghani and Oh (1996)formulated the large-scale disaster relief transportation problem as a multi-commodity multi-modal network flow model with a single objective function. In their conceptual model, the time-varying status of commodities and vehicles mov-ing in a transportation network is represented by three types of links, includmov-ing routmov-ing, transfer, and supply/ demand carry-over links, to facilitate the analysis of the resulting complicated network flow problem. In

Fied-rich et al. (2000), a dynamic combinatorial optimization model is proposed to ﬁnd the optimal resource rescue

schedule with the goal of minimizing the total number of fatalities during the search and rescue (SAR) period, which refers to the first few days after the disaster. Although the model proposed by Fiedrich et al. aims merely to deal with rescue resource allocation problems, their approach is unique in the estimation of fatality probabilities in various rescue scenarios during the SAR period. This may motivate further research expansion to forecast the time-varying relief demand patterns using the relationships between the survivals and the induced relief demands. In addition, three types of affected areas (termed as the operational areas in their study) coupled with three rescue-related facilities are defined in order to clarify the linkages among the affected areas, rescue facilities, and the corresponding work tasks in the resulting resource allocation process.

More recently, considering the complexity and diﬃculty in solving the emergency logistics distribution problem with a single model, there is a research trend of decomposing the original problem into given mutu-ally correlated sub-problems, then solve them systematicmutu-ally in the same decision scheme. For instance, a bi-level hierarchical decomposition approach is proposed inBarbarosoglu et al. (2002)for helicopter mission planning during a disaster relief operation. Therein, the top-level programming model is formulated to deal with the resulting tactical decision problems, covering the issues of helicopter ﬂeet management, crew assign-ment, and the determination of the tour number undertaken by each helicopter, followed by the base-level programming model aiming to address the corresponding operational decisions, including routing, loading/ unloading, rescue, and re-fueling scheduling problems. Another case studied by Ozdarmar et al. (2004) is

(3)

unique in incorporating the vehicle routing problem into the relief distribution process, in which vehicles are treated as commodities to facilitate decomposing the comprehensive emergency logistics distribution problem into two multi-commodity network sub-problems, and then solved using Lagrangean relaxation.

Although the emergency logistics distribution problem considered here is related to vehicle routing prob-lems (VRP) which have been extensively investigated in previous literature, the nature of the problem of a comprehensive emergency logistics distribution system can be more complicated, and needs to further include the pre-route operational tasks, such as relief demand forecasting and collection as well as efficient relief resource allocation to affected areas. In addition, the typical vehicle routing maneuvers, involving the require-ment of vehicle dispatching and returning to the same depot, do not necessarily hold in the emergency logistics context. However, it is maintained in Ozdarmar et al. (2004) that in some emergency logistics operational cases, any given node receiving relief commodities can be the new depot or the former depot may no longer supply relief, and thus vehicles may stay at their last stop, until the next distribution mission is identified. Surveys and discussions on the existing VRP approaches and their extensions can be readily found in the pre-vious literature (Solomon and Desrosiers, 1988; Dror et al., 1989; Laporte, 1992; Bertsimas et al., 1995; Fisher,

1995; Bramel and Simchi-Levi, 1999), and thus are omitted in this study.

Accordingly, this study presents a quick-responsive emergency logistics distribution approach to coordinate the chain-based relief logistics flows in a specified three-layer relief supply network during the crucial rescue period right after the occurrence of a serious natural disaster. Here the crucial rescue period refers to the initial three days following the onset of a disaster, which is regarded as the most critical period to search and rescue the trapped survivals. Aiming at the scope of the study defined above, the proposed emergency logistics co-distribution approach is unique with the following three distinctive features:

(1) The time-varying relief demand associated with each affected area is predicted using a proposed short-term dynamic relief demand forecast model. The respective time-varying relief demand is predicted which mainly relies on the postulation that the relief demand may correlate highly with the number of survivals trapped in the affected areas without receiving any rescue aids. As such, the aforementioned information concerning the number of trapped survivals after a disaster is quite limited, and may change over time upon updating the disaster conditions. Considering both the uncertain and dynamic features of relief demands, a specific forecast mechanism is proposed to predict the time-varying relief demands in the affected areas.

(2) Prior to vehicle dispatching, a fuzzy-based affected-area clustering procedure is executed for the classi-fication of the affected areas according to the estimated degrees of severity. Here, the affected area groups rated with greater severity degrees indicating relatively urgent need for relief will be scheduled to receive relief demands with higher priority.

(3) Due to the distinctive driven operational feature of emergency logistics, a two-stage demand-chain based dynamic optimization model is formulated to integrate the resulting decisions of multi-source relief supply from the ﬁrst to the second layer, and relief distribution from the second to the third layer in the speciﬁed emergency logistics co-distribution framework.

The remainder of the paper is organized as follows. Section2introduces the methodological framework of the proposed approach, including the embedded multi-source relief supply and co-distribution mechanisms and related models proposed for systematic operations. Section3 depicts a numerical study including a real earthquake case occurred in Taiwan, and the corresponding numerical results generated using the proposed method. Finally, concluding remarks and directions for future research are summarized in Section 4. 2. System speciﬁcation

The hypothetical relief logistics network considered in this study is depicted inFig. 1, which involves three primary chain members: (1) relief suppliers, (2) urgent relief distribution centers (URDC), and (3) relief-demanding areas, forming a specific three-layer relief supply chain. Here, relief suppliers refer to the sources of relief supply from the private or public organizations; the URDC is defined as the relief supply hub which aims to efficiently coordinate the resulting inbound and outbound relief logistics in response to the urgent

(4)

relief demands from the aﬀected areas in the crucial rescue period. Generally, considering the speciﬁc logistics requirements and facility scale, such a distribution center is dominated by the public sector, e.g., the local gov-ernment or the corresponding regional rescue organizations.

Based on the specified relief logistics network, the following hypothesis is postulated for system operations. Given the occurrence of a natural disaster, resulting in different degrees of damage associated with certain affected areas, the functionality of the proposed emergency logistics distribution system is triggered immedi-ately by the system initialization phase, where the members of the specified relief logistics network, including the relief suppliers, URDCs, and the affected areas are identified. Meanwhile, initial relief is collected urgently from all identified suppliers to the corresponding URDCs nearby, followed by the execution of a sequential operational phases presented inFig. 2. Herein, each respective URDC kicks off the emergency logistics distri-bution mechanism through the real-time estimation of time-varying relief demands and grouping the affected areas by immediately delivering the remaining relief to the affected areas according to the estimated distribution priority associated with these areas. The relief supply mechanism is then triggered by requesting the pre-spec-ified multiple sources to supply relief urgently to appropriate URDC based on the estimates of relief demands updated in real time. In addition, uncertainties in the operational environment may occur either from the demand side (e.g., relief estimation errors and incomplete disaster damage information) or from the supply side (e.g., diverse lead times from multiple relief supply sources). Each URDC should continue to update the time-varying relief demand and supply information and mutually share the updated information so as to coordinate the time-varying inbound and outbound logistics flows in the relief distribution process. Such a five-phase sequential operational routine is continued until the urgent relief demands from all affected areas are satisfied, and then taken over by any normal relief distribution strategies to maintain the daily needs of the correspond-ing survival. The correspondcorrespond-ing models executed in these phases are detailed in the followcorrespond-ing section. 3. Model formulation

To facilitate model formulation in the following subsections, ﬁve assumptions are postulated.

(1) The number of aﬀected areas (denoted by I) and the corresponding geographic relationships are known. It is presumed that such information can be readily accessible in real time via advanced disaster detection technology.

Relief supply channels Relief distribution chennals layer-1: relief suppliers layer-2: relief distribution layer-3: affected areas

center

(5)

(2) The relief supply channels including the corresponding supply sources and URDC are given. Practically, such speciﬁc distribution channels can be pre-speciﬁed in the corresponding strategic planning domain, followed by national disaster management programs.

(3) The updated information in terms of disaster-induced damage conditions and casualties associated with each aﬀected area can be obtained during the crucial rescue period.

(4) The time-varying relief demand needed in a given aﬀected area is highly correlated with the number of corresponding local survivals.

(5) Different types of relief are allowed to be loaded in a vehicle to serve the affected areas. Correspondingly, each vehicle is permitted to load with multiple types of relief in any given relief distribution mission. Based on the aforementioned postulations, the models embedded in these five sequential operational phases are formulated, respectively, which are described in the following subsections.

3.1. Time-varying relief demand forecast

Relief demand forecasting is the primary step, which is carried out to forecast the time-varying urgent relief demand associated with each aﬀected area. In this study, two types of urgent relief are mainly considered: the ﬁrst is the daily consuming relief including water and meal boxes, which may vary with the time of day, and the second is the daily-used equipment for refugees, e.g., sleeping bags and camps. Following the

(6)

aforementioned fourth assumption and the concept of safety stock prepared to avoid the phenomenon of relief demand over the corresponding supply in each given affected area, the time-varying relief demand forecast model is formulated as Dl iðtÞ ¼ max al_d iðtÞ L þ z1a STDliðtÞ ffiffiffi L p ;0 n o ; if l2 Rconsm; max al_d iðtÞ þ bli Pt1 e¼1A l iðt eÞ; 0 n o ; if l2 Requip; 8 > < > : ð1Þ where aland bl

i are the two parameters representing the average hourly demand of relief l needed per survival and the corresponding buﬀer demand associated with relief l and aﬀected area i; Al

iðt eÞ represents the time-varying amount of relief l arriving at a given aﬀected area i in a given time interval t e; Dl

iðtÞ represents the time-varying relief demand associated with relief l of an affected area i in a given time interval t; L represents the upper bound preset to regulate the temporal headway between two successive relief distributions to any given affected area without exceeding the corresponding maximum value; Rconsmand Requiprepresent the relief groups in terms of daily consuming relief and daily-used equipment, respectively; z1arepresents the respective statistical value chosen given that the tolerable possibility of time-varying relief demand shortage is set to be a; di(t) represents the estimated number of survivals trapped in an affected area i in a given time interval t. STDl_iðtÞ represents the time-varying standard deviation of relief demand associated with relief l and affected area i, which is given by

STDl_iðtÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pt1 e¼0 D l iðt eÞ DliðtÞ 2 q t 1 ; ð2Þ where Dl

iðtÞ represents the time-varying mean value with respect to the time-varying demand D l iðtÞ, and it is given by Dl_iðtÞ ¼a lPt1 e¼0d l iðt eÞ t : ð3Þ

As expressed in Eq.(1), the estimated demands of daily-consuming relief and daily-used equipment diﬀer mainly in the estimation of the buﬀer demands proposed when the potential supply relief shortage problems are considered. For those daily-consuming relief demands, the amount of potential relief shortage during L is given by z1a STDliðtÞ

ffiffiffi L p

, as shown in the upper term of Eq.(1). Such a formulation is built based on the concept of safety stock, which has been extensively implemented for normal inventory control (Simchi-Levi

et al., 2000). Both the factors of relief shortage probability (a) and the upper bound of the relief distribution

headway (L) associated with any given aﬀected area are taken into account such that the following relief demand condition holds:

Prob relief demand during L 6 al_d

iðtÞ L þ z1a STDliðtÞ ffiffiffi L p ¼ 1 a: ð4Þ

In contrast, the refugees may not be consuming the daily used equipment daily, and thus a constant buﬀer demand bl

i is preset to deal with the potential relief deﬁciency problems.

It is noteworthy that the above simplified treatments proposed to address the concern of buffer demands of daily-consuming relief appear reasonable in the operational cases of emergency relief supply. It is because that the time for data-processing and allocating logistics resources during the crucial rescue period are quite lim-ited, resulting in our simplification incorporating the upper bound of the relief distribution headway (L) rather than using variable lead times for the estimation of time-varying relief demand.

3.2. Grouping the aﬀected areas

This phase aims to efficiently group the affected areas using fuzzy-clustering techniques in which each affected area is characterized with multiple attributes to indicate the degree of relief demand. Considering

(7)

the existence of qualitative and quantitative attributes exhibited in characterizing the urgency of relief demands in the affected areas, a hybrid fuzzy-hierarchical clustering method is proposed to perform the corresponding affected-area classification mechanism before dispatching the vehicles. Such a pre-route demand data processing measure can provide the proposed emergency logistical distribution method with ben-efits, not only for allocating the available relief to appropriate resources efficiently, but also to rapidly respond to a variety of relief demands in these affected areas according to their degrees of urgency. The output from this phase, i.e., the identified affected area groups, is used as the input of the following phase to further deter-mine the emergency distribution priority.

In this phase, four urgency attributes used as the determinants of the grouping aﬀected areas are speciﬁed. (1) u1

iðtÞ represents the time difference between the present time interval t and the time of the previous relief arrival to a given affected area i. In emergency logistics distribution operations, the longer the aforemen-tioned time lag, the more urgent the associated affected area that needs relief.

(2) u2

iðtÞ represents the ratio of the number of casualties observed relative to the total number of population trapped in a given affected area i in a given time interval t. In general, a higher casualty ratio associated with a given affected area i may indicate higher urgency of the corresponding affected area in demanding relief to avoid the continuity of casualties in that area.

(3) u3

iðtÞ represents the proportion of the helpless group including children and the elders relative to the total number of population trapped in a given aﬀected area i observed in a given time interval t. According to

Shiono and Krimgold (1989), the survival probability of trapped victims in the disasters may decrease

with time, depending on their physical conditions and severity of injuries. Therein, children and the elders can be regarded as the two most powerless groups who may need relief and rescue most urgently. Therefore, the corresponding proportion is taken into account in determining the emergency distribution priority.

(4) u4

iðtÞ refers to the significance of building damage conditions, such as serious or complete destruction measured in a given time interval t. In general, the damage conditions of the building may reflect the severity of the disaster effects on the survival probabilities of the trapped people. Accordingly, a rela-tively greater degree of building damage conditions may indicate a higher urgency of relief distribution to the corresponding affected area.

Accordingly, we then have a 4· 1 urgency-attribute vector associated with each given aﬀected area i (Ui(t)) as: UiðtÞ ¼ u1iðtÞ; u 2 iðtÞ; u 3 iðtÞ; u 4 iðtÞ T : ð5Þ

Correspondingly, up to this stage, each given aﬀected area can be represented by a multi-attribute datum char-acterized with a 4· 1 urgency-attribute vector.

Using the speciﬁed urgency-attribute vector, a fuzzy clustering-based algorithm is then employed to per-form the function of aﬀected-area grouping. The proposed algorithm is primarily composed of three sequen-tial procedures including: (1) binary transformation, (2) generation of fuzzy correlation matrix, and (3) clustering. The fundamentals of these procedures are described in the following.

3.2.1. Binary transformation

This procedure aims to transform the measurements of each respective urgency-attribute vector into binary data for further fuzzy-clustering. It should be noted that all these specified urgency attributes are not quan-titative. In addition, there still exist certain errors in the corresponding data collected during a disaster. Con-sidering the inefficiency and uncertainty in measuring these attributes using the corresponding raw data, a specific data-transformation procedure is needed.

Here, two transitional steps are involved. First, for each given aﬀected area i, all the corresponding urgency attributes are measured with ﬁve linguistic terms, including ‘‘very high’’, ‘‘high’’, ‘‘medium’’, ‘‘low’’, and ‘‘very low’’ (i.e., VH, H, M, L, VL for short), using the following logic rules:

(8)

w½uk iðtÞ ¼ VH ; if uk iðtÞ P 0:8ukmax; H ; if 0:6uk max6u k iðtÞ < 0:8u k max; M ; if 0:4uk

max6ukiðtÞ < 0:6ukmax; L; if 0:2uk max6ukiðtÞ < 0:4u k max; VL; otherwise; 8 > > > > > > < > > > > > > : 8i; k; t; ð6Þ where uk

iðtÞ represents the raw data in terms of the kth urgency attribute associated with a given aﬀected area i collected in a given time interval t; uk

maxrefers to the expected maximum value of u k

iðtÞ; and w½u k

iðtÞ represents the linguistic measurement associated with uk

iðtÞ. Note that Eq. (4) applies to those quantitative attributes including u1 iðtÞu 2 iðtÞ, and u 3 iðtÞ; whereas u 4

iðtÞ itself is a qualitative attribute, and does not need the speciﬁc lin-guistic measurement process mentioned above. Second, these measured urgency attributes are then trans-formed into binary codes, where each linguistic criterion is represented by a 4-bit binary code, e.g., ‘‘0000’’ for the linguistic term ‘‘very low’’ and ‘‘1111’’ for ‘‘very high’’, as illustrated inTable 1.

Accordingly, any given measured urgency attribute uk

iðtÞ can then be transformed into a binary code (w k iðtÞ) with four bits (i.e., rk

i;jðtÞ for j = 1, 2, 3, and 4), and is given by wk_iðtÞ ¼ rk i;1ðtÞ; r k i;2ðtÞ; r k i;3ðtÞ; r k i;4ðtÞ h i : ð7Þ

In addition, to facilitate data processing in real-world applications, the following standardization procedure with respect to rk

i;jðtÞ is proposed, and the corresponding standardized value of r k i;jðtÞ (~r k i;jðtÞ) is given by ~ rk i;jðtÞ ¼ rk i;jðtÞ rkjðtÞ STD½^rk jðtÞ ; ð8Þ where rk jðtÞ and STD½^r k

jðtÞ correspond to the values of mean and standard deviation with respect to r k i;jðtÞ, respectively. Finally, we have the standardized binary urgency attribute (~wk

iðtÞ), which is given by ~

w_ikðtÞ ¼ ~hrk_i;1ðtÞ; ~rk_i;2ðtÞ; ~rk_i;3ðtÞ; ~rk_i;4ðtÞi: ð9Þ

3.2.2. Generation of fuzzy correlation matrix

Given the number of affected areas I, as postulated in the first assumption, a time-varying I· I fuzzy cor-relation matrix (W(t)) is estimated in this procedure employing the standardized binary urgency attributes measured in the previous procedure, where each element (wpq(t)) of W(t) represents the correlation between a given pair of affected areas p and q. The mathematical forms of W(t) and wpq(t) are given by

WðtÞ ¼ w11ðtÞ w12ðtÞ w13ðtÞ w1IðtÞ w21ðtÞ w22ðtÞ .. . w31ðtÞ . . . .. . .. . .. . ... wI1ðtÞ wNNðtÞ 2 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 5 II ; ð10Þ Table 1

Deﬁnitions of 4-bit binary codes for linguistic criteria

Linguistic criterion Binary code

rk

i;1ðtÞ rki;2ðtÞ rki;3ðtÞ rki;4ðtÞ

Very high (VH) 1 1 1 1

High (H) 1 1 1 0

Medium (M) 1 1 0 0

Low (L) 1 0 0 0

(9)

wpqðtÞ ¼ 1 1 b ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X4 k¼1 X4 j¼1 ~ rk pjðtÞ ~rkqjðtÞ 2 v u u t ; ð11Þ

where b represents a parameter pre-determined for the upper and lower boundaries of wpq(t), i.e., 1 and 0, respectively; j represents a given bit code. It is also worth noting that according to Eqs. (10) and (11), W(t) turns out to be a symmetric matrix.

In addition, according to the fundamentals of fuzzy clustering technologies, the estimated fuzzy correlation matrix W(t) should be processed through the composition operation, such that the following condition (i.e., Eq. (12)) holds:

f

WðtÞ fWðtÞ ¼ fWðtÞ; ð12Þ

where fWðtÞ represents the composite fuzzy correlation matrix of W(t). To obtain fWðtÞ, a routine of the max-min composition operation with respect to each given element of W(t) (e.g., wpq(t)) is proposed as

wpq¼ max I

i¼1 min½wpiðtÞ; wiqðtÞ

: ð13Þ

Such a routine should be conducted until the condition shown in Eq.(12)is satisﬁed. 3.2.3. Clustering

This procedure clusters the identified affected areas into several groups, where the affected areas with rel-atively similar urgency characteristics are assigned to the same group, and relrel-atively, their urgency attributes can be significantly different from those of any other affected area groups. Such a cluster function is executed in each time interval during the crucial rescue period. To execute this mechanism, five major computational steps are involved in the proposed algorithm, as summarized in the following:

Step 0: Initialize the computational iteration. Let the crucial rescue period start with the first time interval (i.e., t = 1); set the column-search index p = 1; input the composite fuzzy correlation matrix (fWðtÞ) measured in the previous procedure; start the iteration from the first column of the processed fuzzy correlation matrix (~w1ðtÞ), and let p = p. Herein, we can target the affected area associated with the greatest number of casualties as the first one to trigger the affected-area grouping procedure. It is also noteworthy that any targeted affected area is regarded as the representative of a given affected-area group in the following cluster process.

Step 1: Given an aﬀected area (denoted by p), remove the row of fWðtÞ associated with p (~wpðtÞ T

). Note that once an aﬀected area is selected as the target customer denoted by ~wpðtÞ, it is not necessary to re-clus-ter the target area in a given time inre-clus-terval t. Accordingly, the corresponding row ~wpðtÞ

T

can be removed to facilitate the following clustering process.

Step 2: Find the maximum element in ~wpðtÞ, denoted by ^wpqðtÞ, and then conduct the following cluster pro-cedures in sequence:

• If the condition ^wpqðtÞ > k1 holds, then assign the aﬀected area q to the same group as the target aﬀected area p, and remove the row of fWðtÞ associated with q (~wqðtÞ

T

). This represents that the given aﬀected area q has the urgency conditions that are very much similar to that of the target aﬀected area p, and thus both of them are assigned to the same group.

• Go back to Step 2 to continue checking other elements of ~wpðtÞ until there does not exist any ele-ment that meets the aforeele-mentioned clustering condition. If so, remove ~wpðtÞ from fWðtÞ, represent-ing that all the elements of ~wpðtÞ have been considered, and thus the corresponding clustering process based on the urgency attributes of the target aﬀected area p can be ended in the given time interval.

• If there are aﬀected areas assigned at this stage, then let all the assigned areas be the target areas (i.e., let p = q), and go back to Step 1 to process the elements of fWðtÞ associated with these target areas. • Otherwise, let p = p, and go back to Step 1 to continue the cluster process for the next aﬀected-area

(10)

Step 3: Conduct the following termination rules to stop the mechanism of aﬀected-area grouping: • If no column remains, then stop the cluster procedure;

• Else, go back to Step 1 for the next iteration.

Herein, k1 is a pre-determined threshold for identifying the relative similarity between a given pair of aﬀected areas, and in practice it can be speciﬁed based on the conditions of available logistics resources and relief supplies.

Suppose that all the given aﬀected areas planned to be served in a given time interval t are clustered into ‘‘G’’ groups through the aforementioned clustering procedure. Thus, up to this stage, we have the group-based standardized binary urgency-attribute matrix associated with a given aﬀected-area group g ( eWgðtÞ) given by

e

WgðtÞ ¼ WigðtÞ; ig¼ 1; 2; . . . ; ‘g

16‘g; ð14Þ

where eWgðtÞ is a (16 · ‘g) group-based standardized binary urgency-attribute matrix which is composed of the standardized binary urgency-attribute vectors (WigðtÞ) of all the aﬀected areas (ig) belonging to group g, and ‘g

represents the number of aﬀected areas involved in a given group g. Here, WigðtÞ is a (16 · 1) vector, which can

be further expressed as WigðtÞ ¼ ~w k igjðtÞ; j ¼ 1; . . . ; 4; k ¼ 1; . . . ; 4 h i 161: ð15Þ

3.3. Determination of distribution priority

In this phase, a respective evaluation measure (Xg(t)) is proposed to determine the relief distribution prior-ity associated with each clustered aﬀected-area group g. Xg(t) is given by

XgðtÞ ¼ P4 k¼1-k P 8ig P4 j¼1w~kigjðtÞ h i ‘g ; ð16Þ

where -k represents the weight associated with a given urgency attribute k, which is specified to determine the relative significance of these attributes in influencing the corresponding degree of urgency.

Using the estimates of Xg(t), the clustered G affected-area groups can then be ranked in order by the cor-responding relief distribution priority, where the affected areas in the same group will be given with the same priority in the relief distribution. Note that in some cases, it is possible that some affected-area groups may have consistent values of Xg(t), in which the mean value of the urgency attribute (i.e.,

P 8ig P4 j¼1w k ig jðtÞ ‘g ) associated

with the highest weight (-k) is suggested to be employed for further distinction. 3.4. Group-based relief distribution

In this phase, a composite weighted multi-objective optimization model is formulated to deal with the prob-lem of distributing the clustered multi-type urgent relief from multiple urgent relief distribution centers to mul-tiple affected-area groups. Each of the given affected-area group g, planned to be served in a given time interval t, is associated with the two respective objective functions with the goals of maximizing the time-varying relief demand fill rate (F1

gðtÞ) and minimizing the time-varying distribution costs (F 2 gðtÞ). Here, F 1 gðtÞ and F 2 gðtÞ are given, respectively, by max F1 gðtÞ ¼ P 8l PM m¼1 P 8ig2gX l m;igðtÞ P 8l P 8ig2gD l igðtÞ ; 8ðg; tÞ; ð17Þ min F2 gðtÞ ¼ X 8l XM m¼1 X 8ig2g CSl m ðu 1 igðtÞÞ þ CT l m;ig h i Xl m;igðtÞ; 8ðg; tÞ; ð18Þ

(11)

where CSl

mand CT

l

m;ig represent, respectively, the unit setup and transportation costs associated with a given

type of relief l and distribution center m for relief distribution to a given aﬀected area ig; M represents the num-ber of relief distribution centers; and Xl

m;igðtÞ is a decision variable denoting the time-varying quantities of a

given type of relief l distributed from a given relief distribution center m to a given aﬀected area igin a given time interval t.

Considering the differing distribution priority determined above, two respective weights (xr and xc) are introduced, subject to the condition xr+ xc= 1. In addition, the difference in measurement scales associated with the relief demand fill rate and distribution costs may also influence the determination of optimal solu-tions. Therefore, the proposed multi-objective functions associated with each given affected-area group g are rewritten as a composite normalized form (max Fg(t)) which is given by

max FgðtÞ ¼ xr F1 gðtÞ F 1 g F1 g F1g " # xc F2 gðtÞ F 2 g F2 g F2g " # ; 8ðg; tÞ; ð19Þ where F1 gand F 1

gare the expected maximum and minimum values associated with F 1 gðtÞ, whereas F 2 gand F 2 gare

the expected maximum and minimum values associated with F2_gðtÞ.

Expanded from the above formulation, the aggregate objective function (eFðtÞ) involving the distribution priority associated with each aﬀected-area group determined in the previous phase is formulated for urgent relief distribution from multiple relief distribution centers to multiple aﬀected-area groups. eFðtÞ is given by

max eFðtÞ ¼X G

g¼1

c_gðtÞ FgðtÞ; 8t; ð20Þ

where cg(t) represents the time-varying weight associated with a given aﬀected-area group g, and is given by c_gðtÞ ¼_P_GXgðtÞ

g¼1XgðtÞ

; 8t: ð21Þ

As deﬁned previously, Xg(t) estimated in the previous phase is used to determine the relief distribution priority associated with each clustered aﬀected-area group g.

In addition, the limitations in terms of relief supply and demand conditions, and the corresponding resource availability should also be considered, and thus the corresponding constraints (Eqs. (22)–(25)) are proposed: XM m¼1 Xl_m;i gðtÞ 6 D l igðtÞ; 8ðg; ig; l; tÞ; ð22Þ XG g¼1 X 8ig2g Xl_m;i_gðtÞ 6 Ql mðt 1jt 1Þ; 8ðl; m; tÞ; ð23Þ XG g¼1 X 8ig2g X 8l UNVl Xlm;igðtÞ 6 X 8t UNLt NVmtðtÞ; 8ðm; tÞ; ð24Þ Xl m;igðtÞ P 0; 8ðg; ig; l; m; tÞ; ð25Þ where Ql

mðt 1jt 1Þ represents the time-varying amount of relief l remaining at a given relief distribution center m observed at the end of a given time interval t 1; UNVlrepresents the unit volume associated with relief l; UNLtrepresents the unit loading capacity associated with a given type of vehicle t; NVmtðtÞ represents the number of vehicle trips associated with a given type of vehicle t available at a given relief distribution cen-ter m in a given time incen-terval t.

Among the above constraints, Eq.(22)is speciﬁed to ensure that the aggregate amount of relief l distributed to a given aﬀected area ig in a given time interval t should not exceed the corresponding estimated relief demand Dl

igðtÞ which is forecasted earlier in the 1

st

operational phase (see Eq. (1)), and similarly, Eq. (23) regulates the time-varying relief amount distributed from any given distribution center without exceeding the corresponding relief amount available in any given time interval. In addition, Eq. (24) indicates that

(12)

the time-varying relief amount distributed from any given relief distribution center should not exceed the cor-responding aggregate transportation capacity, which is determined by the respective unit loading capacity UNLt associated with vehicle t and the corresponding number of vehicle trips NVmtðtÞ. Furthermore, Eq. (25)characterizes a feasible numerical domain associated with each decision variable Xl

m;igðtÞ.

3.5. Dynamic relief supply

This phase aims to transport the optimal relief supply amounts efficiently in multiple relief supply channels (i.e., from multiple relief suppliers to distribution centers) in each given time interval. Differing from the pre-vious relief distribution phase which serves to distribute mixed types of relief in each given relief distribution channel (i.e., a given pair of relief distribution center and affected area), the phase of relief supply aims to transport homogeneous supplied relief in each relief supply channel with the goal of minimizing the transpor-tation costs. Similarly, considering the degree of urgency shown in each given relief distribution center, the objective function is formulated as a composite weighted form with respect to the corresponding transporta-tion costs associated with these relief supply channels. Accordingly, we have

minX 8l XM m¼1 gmðtÞ XNl nl¼1 CSl_n_lþ CTlnl;m Xl nl;mðtÞ " # ; 8t; ð26Þ

where CSl_n_l and CTl_n_l_;mrepresent, respectively, the unit setup and transportation costs associated with a given type of relief l for the relief supply from a given supplier nlto a given distribution center m; nlrepresents a given supply source which speciﬁcally supplies a given type of relief l; Nlis the number of relief suppliers asso-ciated with relief l; Xl_n_l_;mðtÞ is a decision variable representing the time-varying quantities of a given type of relief l supplied from a given relief supply source nlto a given relief distribution center m in a given time inter-val t; gm(t) represents the time-varying weight indicating the urgency degree associated with a given relief dis-tribution center m, and is given by

gmðtÞ ¼ PG g¼1cgðtÞ P 8ig2gXe l m;igðtÞ PM m¼1 PG g¼1cgðtÞ P 8ig2gXe l m;igðtÞ ; 8ðm; tÞ: ð27Þ

In Eq.(27), cg(t) represents the respective time-varying weight associated with a given aﬀected-area group g and can be estimated by Eq. (21); eXl

m;igðtÞ refers to the optimal solution of X

l

m;igðtÞ which is determined in

the previous operational phase.

Similarly, the limitations with respect to the relief supply and demand conditions in the relief supply chan-nels should be considered, and the corresponding constraints are formulated, as shown in Eqs.(28)–(30):

XG g¼1 X 8ig2g e Xl m;igðtÞ 6 XNl nl Xl nl;mðtÞ 6 CAP l m; 8ðl; m; tÞ; ð28Þ XM m¼1 Xl_n l;mðtÞ 6 Q l nlðt 1jt 1Þ; 8ðl; nl; tÞ; ð29Þ Xl_n l;mðtÞ P 0; 8ðl; m; nl; tÞ: ð30Þ

Here, CAPl_mrepresents the storage capacity associated with a given relief distribution center m with respect to a given type of relief l; Ql_n_lðt 1jt 1Þ represents the time-varying amount of relief l remaining at a given relief supply source nlobserved at the end of a given time interval t 1. Eq.(28)is proposed to ensure that the aggregate amount of relief l supplied to a given relief distribution center m in a given time interval t should neither be less than the corresponding outbound distribution amount nor be greater than the corresponding storage capacity. Similarly, the time-varying aggregate relief amount supplied from any given supply source nl in a given time interval t should not exceed the amount of relief remaining at the end of the previous time interval t 1, thus leading to the formulation of Eq.(29), followed by Eq.(30)specifying the feasible numer-ical domain associated with each decision variable Xl

(13)

Through the execution of the aforementioned ﬁve sequential operational phases, the optimal solutions of the time-varying relief supply and distribution ﬂows (i.e., eXl

nl;mðtÞ and eX

l

m;igðtÞÞ are obtained to determine the

corresponding emergency logistics ﬂows in the speciﬁed emergency logistics network in each of the given time intervals. In addition, for the operations in the next time interval, the relief inventory conditions in both sup-ply sources (Ql

nlðtjtÞ) and distribution centers (Q

l mðtjtÞ) should be updated by Ql_n lðtjtÞ ¼ Q l nlðt 1jt 1Þ þ DS l nl XM m¼1 e Xl_n l;mðtÞ; 8ðl; nl; tÞ; ð31Þ Ql_mðtjtÞ ¼ Ql mðt 1jt 1Þ þ XNl nl e Xl nl;mðtÞ XG g¼1 X 8ig2g e Xl m;igðtÞ; 8ðl; m; tÞ; ð32Þ where DSl

nl represents the constant production increment of relief l generated at a given relief supply source nl

in each given time interval. 4. Numerical results

The main purpose of this numerical study is to demonstrate the efficiency of the proposed method that can be used in practical operations of quick-responsive urgent relief distribution during large-scale natural disas-ters. The case study illustrates a massive earthquake measured 7.3 on the Richter scale, which occurred in cen-tral Taiwan, on September 21, 1999 (termed as the 921 Chichi earthquake). Given the disaster-related information including the number, locations, and severities of the affected areas, as well as the relief supply network, the proposed method is employed to determine the time-varying relief logistics flows in a three-day crucial rescue period. Here the unit length of a time interval is given by 4 hours, thrice a three-day. There are a total of 12 time intervals involved in the test scenario. The resulting performance is then evaluated by comparing with the historical data which was yielded using the existing emergency logistics management strategies.

The problem background of the case is summarized below. Reportedly, the earthquake and its aftershocks caused 2455 deaths in total, more than 8000 injuries, with 38,935 homes destroyed (Executive Yuan of Taiwan, 1999). The affected areas, Taichung and Nantou Counties, are located in central Taiwan. Here, the study is aimed at 24 most severely affected (in terms of the dead casualties) areas of Taichung. To centralize the actions of rescue and relief distribution to the affected areas, three tentative refuge centers (termed as relief distribu-tion centers 1, 2, and 3) to collect and supply relief were erected at the three towns (Dongshi, Shinkang, and Wufan) in Taichung Country. Meanwhile, the support of trapped-casualty rescue and relief supply from 6 unaffected counties, including Taipei, Taoyuan, Hsinchiu, Chengghua, Tainan, and Kaohsiung (termed as supply sources 1 to 6), was requested immediately by the corresponding refuge centers. Nevertheless, due to the lack of coordination between the refuge centers and relief supply sources as well as the overestimation of relief demands from affected areas, there arose serious relief supply–demand imbalance problems. One typ-ical example is the excessive accumulation of supplied relief at these refuge centers without proper modulation during the rescue period, as illustrated in Fig. 3. In addition, the allocation of relief distribution resources, such as vehicles and volunteers, and the corresponding vehicle dispatching strategies implemented at these ref-uge centers became disordered, resulting in significant delays in transporting relief to certain affected areas.

Based on the background of the above problem, a simpliﬁed 6· 3 · 24 relief supply network is formed in the test scenario, as shown inFig. 4, where the geographical relationships among these relief demand and sup-ply units are speciﬁed.

The next step is to generate the time-varying disaster-induced casualty and damage data used for relief demand forecasting and affected-area grouping. Due to the lack of the corresponding time-varying casualty and damage data (e.g., the number of casualties observed in each of the given affected areas in each time inter-val), the simulated data processed from the official statistics in the 921 earthquake special report (Executive

Yuan of Taiwan, 1999) were used in the test scenario. Table 2 summarizes the aggregate statistics in terms

of the disaster effects as well as the corresponding population data associated with the affected areas of the study site. Conveniently, the time-varying number of casualties associated with each affected area observed in each given time interval was assumed to follow a Poisson process with the real mean value obtained from

(14)

the corresponding statistics, and thus, estimated via simulation. Accordingly, the simulated time-varying casu-alties and damage data used in the three-day test scenario were generated.

Fig. 3. Illustration of relief over-supply conditions at refuge centers.

(15)

In addition, the determination of model parameters must precede the execution of relief distribution. Here, the corresponding parameters were set mainly for three diﬀerent scenarios: (1) relief demand forecast, (2) aﬀected area grouping, and (3) relief distribution, which are summarized inTables 3 and 4.

Regarding relief demand forecast, the key parameters estimated include al, bl_i, a and L, where alwas deter-mined based on the basic requirement of personal survival under disaster conditions; and the others were pre-set considering the efficiency and the corresponding resource limitations of the relief distribution operations given the size of the existing distribution network. For instance, considering the emergency of relief supply to the affected areas, the corresponding relief shortage probability (a) and the upper bound of the relief dis-tribution headway (L) were set, respectively, to be 0.01 and 6 hours to ensure that 99% of the time-varying relief demand associated with each of the given affected areas can be satisfied in 6 hours under the proposed relief distribution operations.

The parameters used for grouping those affected areas involve b, k1, and -k, where b and k1are clustering-related parameters which were calibrated with the time-varying disaster-induced effect data generated previously; -k, as mentioned in Section3.2, represents the defining the relative significance of a given urgency attribute k in determining the distribution priority of the affected areas, and is tentatively set to be mutually equal (i.e., -k = 0.25, k = 1, 2, 3, 4) in the case study.

The parameters present in the phase of relief distribution are classiﬁed into three groups: (1) cost-related parameters, (2) weights (i.e., xrand xc), and (3) the capacities of the corresponding facilities and resources. Here, the capacity-related parameters such as the capacities of facility storage (CAPl_m) and vehicle loading (UNLt) were determined mainly based on the operational requirements of general logistics distribution cases. In contrast, the cost-related parameters may need to consider the factors of relief availability and the corresponding accessibility in each given distribution channel under disaster conditions. As a result, they were set to be greater than the normal unit costs in this scenario. Considering the tremendous number of the unit transportation costs which need to be set in this scenario, they are summarized speciﬁcally in Table 4 for clarity.

Table 2

Statistics of disaster eﬀects and population of the aﬀected areas

Aﬀected area Population Number of helpless people Casualty Building damage condition

area-1 1572 230 100 VS area-2 5377 1106 42 S area-3 1764 327 33 VS area-4 9330 2369 251 S area-5 11 149 3110 125 S area-6 26 581 8864 172 VS area-7 3874 961 136 VS area-8 2398 552 44 S area-9 1084 283 69 S area-10 870 199 28 S area-11 1672 347 14 M area-12 1304 265 47 M area-13 2899 330 71 VS area-14 961 314 19 S area-15 4385 1081 106 S area-16 2890 477 38 M area-17 1746 293 62 M area-18 9881 2019 17 S area-19 27 759 8910 237 S area-20 5633 1992 89 VS area-21 3701 656 25 S area-22 11 933 3712 83 S area-23 2908 983 27 M area-24 1856 507 48 S

(16)

Based on the generated time-varying disaster effect data and predetermined parameters, the proposed relief distribution method was then used to determine the optimal time-varying relief logistics flows in the specified relief supply network. Since the method proposed is a hybrid fuzzy-based optimization approach, the computational procedures which embed the fundamentals of the fuzzy clustering and multi-objective optimization models are needed to search for the optimal solutions. Conveniently, two computational stages were carried out to determine the optimal relief logistics flows in the numerical study. First, a fuzzy cluster-ing-based algorithm written in the Turbo C++ language was executed to forecast the time-varying relief demands and cluster affected areas into several groups with specific distribution priority, as described in the operational phases 1, 2, and 3 of the proposed methodology. Using the existing LINGO 8.0 package as an optimal-solution tool, the multi-objective models embedded in phases 4 and 5 were carried out to search the optimal solutions in terms of the time-varying relief logistics flows in the specified relief supply network.

The numerical results of the aﬀected-area grouping are summarized inTable 5to demonstrate the perfor-mance of the proposed method in determining the aﬀected-area groups and the corresponding distribution Table 3

(17)

Estimated unit transportation costs

m ig

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

CTl

m;ig (US$; l = 1: water)

1 0.8 0.8 0.8 0.8 0.8 0.8 0.8 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2

2 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1

3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7

CTl_m;i_g (US$; l = 2: meal box)

1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9

2 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

3 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

CTl

m;ig (US$; l = 3: sleeping bag)

1 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 6 6 6 6 6 6 6 6 6

2 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 6 6 6 6 6 6 6 6 6

3 6 6 6 6 6 6 6 4 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3

CTl_m;i_g (US$; l = 4: camp)

1 10 10 10 10 10 10 10 12 12 12 12 12 12 12 12 13 13 13 13 13 13 13 13 13 2 12 12 12 12 12 12 12 10 10 10 10 10 10 10 10 13 13 13 13 13 13 13 13 13 3 13 13 13 13 13 13 13 12 12 12 12 12 12 12 12 10 10 10 10 10 10 10 10 10 CT1 n1;m CT 2 n2;m CT 3 n3;m CT 4 n4;m nl m nl m nl m nl m 1 2 3 1 2 3 1 2 3 1 2 3 CTl nl;m(US$) 1 0.3 0.2 0.3 1 0.5 0.5 0.6 1 5 5 6 1 8 6 7 2 0.2 0.1 0.2 2 0.3 0.3 0.4 2 4 3 4 2 6 5 6 3 0.1 0.1 0.1 3 0.3 0.2 0.2 3 3 2 3 3 5 5 4 4 0.1 0.1 0.1 4 0.2 0.1 0.1 4 3 2 2 4 4 4 3 5 0.2 0.2 0.2 5 0.2 0.2 0.2 5 3 3 3 5 5 5 4 6 0.2 0.2 0.1 6 0.3 0.2 0.2 6 4 3 4 6 7 6 6 J.-B. Sheu / Transport ation Research Part E 4 3 (2007) 687–709 703

(18)

priority. As can be seen inTable 5, it appears that the clustered affected-area groups and the associated dis-tribution priority remain to be the same with time. This is due to the corresponding disaster-induced effects, e.g., the casualties and damage conditions reported in the affected areas, which tend to be stabilized in a cer-tain time after the disaster, thus contributing to the cluster results converging to a cercer-tain extent two days after the disaster. In addition, it is noteworthy that the original relief distribution strategies implemented in these three distribution centers might not consider the corresponding distribution urgency associated with the affected areas. Correspondingly, the supplied relief was distributed, either once or twice a day, to these affected areas merely based on their geographical adjacency with respect to the relief distribution centers. Conse-quently, some areas with relatively higher urgency such as demand nodes 4 and 7 presented inFig. 3might suffer from more disaster impacts due to longer delay in relief distribution.

After the aforementioned affected-area grouping and distribution priority determination phases, the urgent relief supply and distribution operations were implemented following the optimal solutions determined by the models embedded in the proposed 4 and 5 phases. The problems defined in both the phases of relief distribu-tion and supply are dynamic linear programming (DLP) models formulated with a total of 3240 (i.e., 4· (6 · 3 + 3 · 24)) decision variables, involving 999 and 234 constraints, respectively. To determine the time-varying relief logistics flows in each given time interval, the proposed models were coded following the format of the LINGO package. Then, the corresponding local optimal solutions of a given time interval were searched by employing the algorithms embedded in LINGO, and used for the recursive estimation in the next time interval. Integrated with the previous operational phases, such a routine is continued until the end of the test period.

Table 5

Summary of the clustered aﬀected-area groups

Time interval t = 1 Time interval t = 2 Time interval t = 3

Group g with distribution priority Group member ig (demand node) Group g with distribution priority Group member ig (demand node) Group g with distribution priority Group member ig (demand node) Crucial rescue period: Day-1

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

Crucial rescue period: Day-2

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

Crucial rescue period: Day-3

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

(19)

It is worth mentioning that the computational inefficiency may not be a significant issue existing in these two phases because the process of searching optimal solutions is recursive in the proposed DLP models. That is, the entire three-day test period is divided into 9 intervals; for each given time interval, the time-varying opti-mal solutions of decision variables are determined mainly based on the present logistics resources and inven-tories updated by the corresponding optimal solutions of the previous time interval. Therefore, only a subset of decision variables and constraints are involved in each given time interval. In addition, according to our previous research experience in fuzzy clustering (Hu and Sheu, 2003;Sheu, 2006), the affected-area grouping procedures conducted in the previous phases may also facilitate searching the local optimal solutions. Note that through the fuzzy clustering process, all the affected areas are classified into certain groups, where the affected areas with similar urgency attributes are bounded together so as to respond efficiently to their relief demands based on the respective distribution priority. Accordingly, the time-consuming computational prob-lems were not found in the study case.

To quantitatively assess the eﬃciency of the proposed method, particularly in quickly responding to relief demands in the aﬀected areas and coordinating multiple relief supply sources in diverse disaster severity sce-narios, three criteria are proposed:

(1) AT , which represents the average time diﬀerence between two successive relief arrivals to a given aﬀected area in a day.

(2) AF , which represents the average relief demand ﬁll rate during the three-day test period.

(3) TC, which represents the total emergency logistics costs, mainly including the inbound and outbound transportation costs and the corresponding increase in inventory costs at the relief distribution centers in the three-day test period.

In addition, to evaluate the relative performance of the proposed method, we compared the numerical results with those obtained under the condition that all the relief distribution centers and supply sources are isolated in the emergency logistics distribution operations. Correspondingly, by implementing the existing strategy, each relief distribution center may merely serve certain aﬀected areas located in the same town with-out any coordination with the other distribution centers and relief supply sources, mimicking the existing strategy implemented at that time. The results are summarized inTable 6.

Overall, the comparison results ofTable 6indicate the potential advantages of the proposed approach in urgent relief distribution. Using the proposed method, the aggregate relative improvement reaches as high as 30.6%, outperforming the existing strategy either in the eﬃciency of relief distribution in response to the urgent relief demands of aﬀected areas (measured by AT and AFÞ or in saving the operational costs of emer-gency logistics (measured by TCÞ. This may also imply that the proposed relief distribution methodology cou-pled with the operational strategies is promising for practical uses in emergency logistics management under large-scale disaster conditions.

In addition, three supportive generalizations are made as follows for further discussion.

First, as can be seen inTable 6, the aggregate relative improvement in the proposed system performance results mainly from the time saving in continuously distributing relief to the affected areas during the crucial rescue period. In the study case, the average time headway of relief supply to a given affected area is 4.6 hours, which is improved significantly by 38.7%, upon employing the proposed method. Such a numerical result is meaningful particularly for the applications in emergency logistics management. It should be noted that dif-fering from general business logistics and supply chain management, the efficiency of relief supply to affected

Table 6

Comparison of system performance

Criteria strategy AT (h) AF (%) TC(US$: million)

The proposed approach 4.6 74 15.8

The existing strategy 7.5 58 21.2

Relative improvement (%) 38.7 27.6 25.5

(20)

areas determines not only the operational performance of the emergency logistics system in the supply side but also the survival of trapped people in the affected areas. From a psychological point of view, shortening the time headway of the supplied relief arrivals to affected areas may create the image of the governmental attempt in rescue, and also firm up the will power of trapped people, thus stabilizing the disaster effects.

Second, through appropriate affected-area clustering coupled with the identification of relief distribution priority, the affected areas in the same group can be served more efficiently based on the similarity of their urgency attributes. As observed inTable 5, the 24 targeted affected areas are classified into 6 or 7 groups in most time intervals, and served in the order of urgency using the proposed method. Consequently, 74% of the relief demands associated with the affected area in each group can be satisfied within 5 hours in a day on average, which is superior to the performance of the existing strategy. To a certain extent, this implies that higher emergency logistics service performance can be achieved using the proposed relief distribution methodology.

Third, the coordination between the layers of relief supply and demand through the relief distribution cen-ters is vital in relief logistics control. As mentioned previously, the problem of relief supply–demand imbalance is a common critical issue existing in emergency logistics management. However, through the integration of relief demand forecast and demand-driven pull-based relief supply strategies, the aggregate emergency logistics costs can be reduced by 25.5% in the study case. An induced beneﬁcial eﬀect is that the relief over-supply con-ditions exhibited in the previous relief supply cases (see Fig. 3) no longer exist upon using the proposed method.

Considering the diversity of disaster eﬀects and relief supply conditions in real operational cases, in the fol-lowing tests, nine experimental scenarios by adjusting the corresponding parameters and input data were designed to investigate the applicability of the proposed method under various extreme operational condi-tions. The description of the designed experimental scenarios and the corresponding parameters/input-data adjusted are summarized inTable 7. Using the same evaluation measures, i.e., AT , AF , and TC, the corre-sponding numerical results associated with these scenarios are summarized inTable 8, in which the original performance of the proposed method shown in the previous scenario is included to indicate the relative per-formance of these scenarios.

The numerical results inTable 8are summarized as in the following:

(1) The relief demand fill rate and the corresponding logistics costs appear to have a significant trade-off relationship in influencing the system performance of the relief distribution. As shown in the results of experimental scenario 1, to increase the average relief demand fill rate from 0.74 to 0.87, the aggregate Table 7

Description of the experimental scenarios

Scenario Scenario description Model adjustment

1 The relief demand ﬁll rate is fully considered in relief distribution. Correspondingly, the relief distribution costs are not taken into account during the crucial rescue period

xris set to be 1 in Eq.(1)

2 The relief over-supply to aﬀected areas is allowed during the crucial rescue period

The corresponding constraint (Eq.(22)) is removed

3–5 The number of vehicles available at each relief distribution center in a given time interval is quite limited

The mean value of vehicle trips at each relief distribution center (E½NVm

tðtÞ) is reduced by 25%, 50%, and 75% 6 The storage capacity of each relief distribution center is not

considered. Correspondingly, the continuous relief supply to the distribution center is supposed to be allowed during the crucial rescue period, and thus, diverse relief accumulation conditions, including over-accumulation, may occur at the relief distribution centers

Remove the upper bound (CAPl_m) from Eq.(28)and re-set it to be inﬁnity

7–9 The deﬁciency in supplied relief may occur either at supply sources or at the relief distribution centers in the initial time interval. This may occur practically at the onset of a disaster

Reset the corresponding initial inventories Ql

nlð0j0Þ and

Ql

mð0j0Þ to be under certain lower levels (e.g., 20%, 40%, and 60%)

(21)

relief logistics costs can be raised about 3.8 million US dollars. This implies that the cost of improving the average relief demand ﬁll rate by 0.1 during the crucial rescue period is about 3 million US dollars in this case study. Nevertheless, it should be noted that the upper bound of the average relief demand ﬁll rate is 0.87 subject to the relief supply conditions and the availability of supply resources, e.g., the num-ber of vehicles available in each given time interval.

(2) The constraint in terms of the upper bound of the relief distributed to a given affected area may not have a significant effect on the improvement of the relief distribution performance. As can be seen from the experimental scenario 2, removing the corresponding constraint from the proposed model does not seem to be significantly beneficial for the improvement of either the average relief demand fill rate or the dis-tribution time saving. It is inferred that the corresponding performance may also depend on both the time-varying relief storage and distribution resource allocation conditions.

(3) The number of vehicles available at each relief distribution center appears to be a critical factor in deter-mining the system performance of the relief distribution. As observed from the results of experimental scenarios 3 to 5, the reduction in the number of vehicles associated with each relief distribution center has caused a significant negative effect on the entire system performance, particularly in both the average time headway of relief arrivals to a given affected area (AT ) and the corresponding demand fill rate (AF ). This may infer that the sufficient number of vehicles serving relief distribution to affected areas coupled with appropriate vehicle dispatching strategies is a prerequisite for emergency logistics management. (4) Due to the limitation of logistics distribution resources, e.g., vehicles and servers, the continuous relief

supply strategy does not help to improve the average relief demand ﬁll rate, but raises the costs of relief supply and the diﬃculty in managing the accumulated relief instead. Consequently, the corresponding logistics cost is raised by 35.4%, as observed in scenario 6.

(5) The deficiency of relief supply in the initial time interval does not seem to be a significant factor influ-encing the system performance of the proposed method. As can be seen in the results of experimental scenarios 7 to 9, such relief deficiency conditions at the onset of a disaster may cause a delay in distrib-uting relief to affected areas, and influence the corresponding demand fill rate in the first time interval. However, through the coordination of relief distribution centers and appropriate relief distribution strat-egies using the proposed method, the system performance can still be under control to a certain extent. Overall, the measures of relative improvement shown above have indicated that the proposed method per-mits serving as a decision support tool of emergency logistics distribution, particularly for quick response to the urgent need of relief in the large-scale disaster affected areas.

Table 8

Numerical results of the experimental scenarios

Criteria experimental scenarios AT (h) AF (%) TC(US$: million)

The previous scenario 4.6 74 15.8

Scenario-1 4.4 87 19.6 Relative performance (%) 4.3 17.6 24.1 Scenario-2 4.5 77 16.5 Relative performance (%) 2.2 4.1 4.4 Scenario-3 5.3 72 14.7 Relative performance (%) 15.2 2.7 7.0 Scenario-4 6.2 63 16.6 Relative performance (%) 34.8 14.9 5.1 Scenario-5 7.7 51 18.9 Relative performance (%) 67.4 31.1 19.6 Scenario-6 4.6 74 21.4 Relative performance (%) 0.0 0.0 35.4 Scenario-7 4.7 72 15.6 Relative performance (%) 2.2 2.7 1.3 Scenario-8 4.9 69 15.1 Relative performance (%) 6.5 6.8 4.4 Scenario-9 5.3 68 14.7 Relative performance (%) 15.2 8.1 7.0

(22)

5. Concluding remarks

This paper has presented a novel emergency logistics distribution approach for quickly responding to the urgent relief demands of the affected areas during a three-day crucial rescue period. The proposed approach involves five major mechanisms: (1) time-varying relief demand forecasting, (2) affected-area grouping, (3) determination of distribution priority, (4) group-based relief distribution, and (5) dynamic relief supply. The methodologies mainly used in this study include fuzzy clustering and multi-objective dynamic program-ming models. Based on a proposed three-layer relief supply network, the time-varying relief demand associ-ated with each given affected area is forecasted using the above mentioned model, followed by a fuzzy clustering-based approach formulated to group these affected areas and to associate the respective distribution priority with them. A two-stage demand-driven relief supply strategy is conducted, based on the optimal solu-tions obtained from the proposed relief supply and distribution models.

A numerical study with a real large-scale earthquake disaster occurred in Taiwan was conducted to illus-trate the applicability and potential advantages of the proposed method. By comparing the performance of the proposed logistics distribution method with that of the existing distribution strategy executed by the refuge centers, the numerical results revealed that the overall emergency logistics system performance could be improved by up to 30.6%, resulting mainly from the signiﬁcant improvement in the average time headway of relief arrivals to the aﬀected areas. In addition, nine experimental scenarios with parameter adjustments were conducted.

Nevertheless, there is still a great potential for improving the operations performance of emergency logis-tical distribution. One example is integrating more elaborate emergency logistics resource allocation method-ologies and vehicle dispatching strategies. Furthermore, the provision of real-time and accurate information in terms of survivals and disaster-induced damage conditions in traﬃc networks is needed so as to update the time-varying relief demands and adjust the relief distribution priority in the operations of quick-responsive emergency logistics distribution. In addition, expanding the present study scope and methodology to the inter-national emergency logistics domain is also of our interest, and warrants further research. Such model exten-sion is important particularly for the use of transnational emergency logistics and rescue in a large-scale disaster event, e.g., the earthquake and tsunami event in Southeast Asia.

Finally, it is expected that the proposed emergency logistics distribution approach can make beneﬁts avail-able not only for improving the performance of emergency logistics management, but also for clarifying the importance of the coordination among the relief supply members, e.g., the relief supply sources and refuge centers, in quick response to the real needs of the aﬀected areas in a relief supply network.

Acknowledgement

This research was supported by the National Science Council grants, Taiwan, NSC 95-2416-H-009-004. The author wishes to thank the anonymous referees for their helpful comments.

References

Ardekani, S.A., Hobeika, A., 1988. Logistics problems in the aftermath of the 1985 Mexico City earthquake. Transportation Quarterly 42, 107–124.

Barbarosoglu, G., Ozdamar, L., Ahmet, C., 2002. An interative approach for hierarchical analysis of helicopter logistics in disaster relief operations. European Journal of Operational Research 140, 118–133.

Bertsimas, D.J., Chervi, P., Peterson, M., 1995. Computational approaches to stochastic vehicle routing problems. Transportation Science 29, 342–352.

Bramel, J., Simchi-Levi, D., 1999. The capacitated VRP with unequal demands. In: Bramel, J., Simchi-Levi, D. (Eds.), The Logic of Logistics. Springer, New York, pp. 81–106.

Brown, G.G., Vassiliou, A.L., 1993. Optimizing disaster relief: real-time operational and tactical decision support. Naval Research Logistics 40, 1–23.

Dror, M., Laporte, G., Trudeau, P., 1989. Vehicle routing with stochastic demands: properties and solution frameworks. Transportation Science 23, 166–176.